Bug Summary

File:programs/Simulation/HDGeant/hitCDC.c
Location:line 305, column 7
Description:Value stored to 'status' is never read

Annotated Source Code

1/*
2 * hitCDC - registers hits for Central Drift Chamber
3 *
4 * This is a part of the hits package for the
5 * HDGeant simulation program for Hall D.
6 *
7 * version 1.0 -Richard Jones July 16, 2001
8 *
9 * changes: Wed Jun 20 13:19:56 EDT 2007 B. Zihlmann
10 * add ipart to the function call hitCentralDC
11 */
12
13#include <stdlib.h>
14#include <stdio.h>
15#include <math.h>
16
17#include <HDDM/hddm_s.h>
18#include <geant3.h>
19#include <bintree.h>
20#include <gid_map.h>
21
22#include "calibDB.h"
23extern s_HDDM_t* thisInputEvent;
24
25typedef struct {
26 int writeenohits;
27 int showersincol;
28 int driftclusters;
29} controlparams_t;
30
31extern controlparams_t controlparams_;
32
33void gpoiss_(float*,int*,const int*); // avoid solaris compiler warnings
34
35// Drift speed 2.2cm/us is appropriate for a 90/10 Argon/Methane mixture
36static float DRIFT_SPEED = 0.0055;
37static float TWO_HIT_RESOL = 25.;
38static int MAX_HITS = 1000;
39static float THRESH_KEV = 1.;
40static float THRESH_MV = 1.;
41static float STRAW_RADIUS = 0.776;
42static float CDC_TIME_WINDOW = 1000.0; //time window for accepting CDC hits, ns
43static float ELECTRON_CHARGE =1.6022e-4; /* fC */
44static float GAS_GAIN = 1e5;
45
46binTree_t* centralDCTree = 0;
47static int strawCount = 0;
48static int pointCount = 0;
49static int stripCount = 0;
50static int initialized = 0;
51static float cdc_drift_time[78];
52static float cdc_drift_distance[78];
53static float BSCALE_PAR1=0.;
54static float BSCALE_PAR2=0.;
55
56/* void GetDOCA(int ipart, float x[3], float p[5], float doca[3]); disabled 6/24/2009 */
57
58typedef int (*compfn)(const void*, const void*);
59extern void polint(float *xa, float *ya,int n,float x, float *y,float *dy);
60
61// Sort function for sorting clusters
62int cdc_cluster_sort(const void *a,const void *b) {
63 const s_CdcStrawTruthHit_t *ca=a;
64 const s_CdcStrawTruthHit_t *cb=b;
65 if (ca->t < cb->t)
66 return -1;
67 else if (ca->t > cb->t)
68 return 1;
69 else
70 return 0;
71}
72
73// Simulation of the ASIC response to a pulse due to a cluster
74double asic_response(double t) {
75 double func=0;
76 double par[11]={-0.01986,0.01802,-0.001097,10.3,11.72,-0.03701,35.84,
77 15.93,0.006141,80.95,24.77};
78 if (t < par[3]) {
79 func=par[0]*t+par[1]*t*t+par[2]*t*t*t;
80 }
81 else {
82 func+=(par[0]*par[3]+par[1]*par[3]*par[3]+par[2]*par[3]*par[3]*par[3])
83 *exp(-(t-par[3])*(t-par[3])/(par[4]*par[4]));
84 func+=par[5]*exp(-(t-par[6])*(t-par[6])/(par[7]*par[7]));
85 func+=par[8]*exp(-(t-par[9])*(t-par[9])/(par[10]*par[10]));
86 }
87 return func;
88}
89
90// Simulation of signal on a wire
91double cdc_wire_signal(double t,s_CdcStrawTruthHits_t* chits) {
92 int m;
93 double asic_gain=0.5; // mV/fC
94 double func=0;
95 for (m=0; m < chits->mult; m++) {
96 if (t > chits->in[m].t) {
97 double my_time=t-chits->in[m].t;
98 func+=asic_gain*chits->in[m].q*asic_response(my_time);
99 }
100 }
101 return func;
102}
103
104void AddCDCCluster(s_CdcStrawTruthHits_t* hits, int ipart, int track, int n_p,
105 float t, float xyzcluster[3])
106{
107 // measured charge
108 float q=0.;
109
110 // drift radius
111 float dradius=sqrt(xyzcluster[0]*xyzcluster[0]+xyzcluster[1]*xyzcluster[1]);
112
113 // Find the drift time for this cluster. Drift time depends on B:
114 // (dependence derived from Garfield calculations)
115 float B[3],Bmag,x[3];
116 transformCoord(xyzcluster,"local",x,"global")transformcoord_(xyzcluster,"local",x,"global",strlen("local")
,strlen("global"));
117 gufld_db_(x,B);
118 Bmag=sqrt(B[0]*B[0]+B[1]*B[1]+B[2]*B[2]);
119 float d2=dradius*dradius;
120 float d3=dradius*d2;
121 int i=(int)(dradius/0.01);
122 float my_t,my_t_err;
123 // Check for closeness to boundaries of the drift table
124 if (i>=75){
125 // Do a crude linear extrapolation
126 my_t=cdc_drift_time[75]+((cdc_drift_time[77]-cdc_drift_time[75])/0.02)
127 *(dradius-cdc_drift_distance[75]);
128 }
129 else{
130 int index;
131 if (i<1) index=0;
132 else index=i-1;
133 // Interpolate over the drift table to find an approximation for the drift
134 // time
135 polint(&cdc_drift_distance[index],&cdc_drift_time[index],4,dradius,&my_t,
136 &my_t_err);
137 }
138 float tdrift=my_t/(1.-BSCALE_PAR1-BSCALE_PAR2*Bmag);
139
140 //Longitudinal diffusion
141 int two=2;
142 float rndno[2];
143 grndm_(rndno,&two);
144 float rho = sqrt(-2*log(rndno[0]));
145 float phi = rndno[1]*2*M_PI3.14159265358979323846;
146 float dt=(7.515*dradius-2.139*d2+12.63*d3)*rho*cos(phi);
147 tdrift+=dt;
148
149 // Prevent unphysical times (drift electrons arriving at wire before particle
150 // passes the doca to the wire)
151 double v_max=0.08; // guess for now based on Garfield, near wire
152 double tmin=dradius/v_max;
153 if (tdrift < tmin) {
154 tdrift=tmin;
155 }
156 float total_time=t+tdrift;
157
158 // Skip cluster if the time would go beyond readout window
159 if (total_time > CDC_TIME_WINDOW)
160 return;
161
162 // Average number of secondary ion pairs for 50/50 Ar/CO2 mixture
163 float n_s_per_p=1.94;
164 if (controlparams_.driftclusters == 0) {
165 /* Total number of ion pairs. On average for each primary ion
166 pair produced there are n_s secondary ion pairs produced. The
167 probability distribution is a compound poisson distribution
168 that requires generating two Poisson variables.
169 */
170 int n_s,one=1;
171 float n_s_mean = ((float)n_p)*n_s_per_p;
172 gpoiss_(&n_s_mean,&n_s,&one);
173 int n_t = n_s+n_p;
174 q = ((float)n_t)*GAS_GAIN*ELECTRON_CHARGE;
175 }
176 else {
177 // Distribute the number of secondary ionizations for this primary
178 // ionization according to a Poisson distribution with mean n_s_over_p.
179 // For simplicity we assume these secondary electrons and the primary
180 // electron stay together as a cluster.
181 int n_s, one=1;
182 gpoiss_(&n_s_per_p,&n_s,&one);
183 // Energy deposition, equivalent to anode charge, in units of fC
184 q = GAS_GAIN*ELECTRON_CHARGE*(float)(1+n_s);
185 }
186
187 // Add the hit info
188 int nhit;
189 for (nhit = 0; nhit < hits->mult; nhit++) {
190 if (fabs(hits->in[nhit].t - total_time) < TWO_HIT_RESOL) {
191 break;
192 }
193 }
194 if (nhit < hits->mult) { /* merge with former hit */
195 /* Use the time from the earlier hit but add the charge*/
196 hits->in[nhit].q += q;
197 if (hits->in[nhit].t > total_time) {
198 hits->in[nhit].t = total_time;
199 hits->in[nhit].d = dradius;
200 hits->in[nhit].itrack = gidGetId(track);
201 hits->in[nhit].ptype = ipart;
202 }
203
204 /* hits->in[nhit].t =
205 (hits->in[nhit].t * hits->in[nhit].q + tdrift * dEsum) /
206 (hits->in[nhit].q += dEsum);
207 */
208 }
209 else if (nhit < MAX_HITS) { /* create new hit */
210 hits->in[nhit].t = total_time;
211 hits->in[nhit].q = q;
212 hits->in[nhit].d = dradius;
213 hits->in[nhit].itrack = gidGetId(track);
214 hits->in[nhit].ptype = ipart;
215
216 hits->mult++;
217 }
218 else {
219 fprintf(stderrstderr,"HDGeant error in hitCentralDC: ");
220 fprintf(stderrstderr,"max hit count %d exceeded, truncating!\n",MAX_HITS);
221 }
222}
223
224/* register hits during tracking (from gustep) */
225
226void hitCentralDC (float xin[4], float xout[4],
227 float pin[5], float pout[5], float dEsum,
228 int track, int stack, int history, int ipart )
229{
230 float x[3], t;
231 float dx[3], dr;
232 float dEdx;
233 float xlocal[3];
234 float xinlocal[3];
235 float xoutlocal[3];
236 float dradius,drin,drout;
237 float trackdir[3];
238 float alpha;
239
240 if (!initialized) {
241 mystr_t strings[50];
242 float values[50];
243 int nvalues = 50;
244 int status = GetConstants("CDC/cdc_parms", &nvalues, values, strings);
245
246 if (!status) {
247 int ncounter = 0;
248 int i;
249 for ( i=0;i<(int)nvalues;i++) {
250 //printf("%d %s \n",i,strings[i].str);
251 if (!strcmp(strings[i].str,"CDC_DRIFT_SPEED")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[i].str) && __builtin_constant_p ("CDC_DRIFT_SPEED"
) && (__s1_len = strlen (strings[i].str), __s2_len = strlen
("CDC_DRIFT_SPEED"), (!((size_t)(const void *)((strings[i].str
) + 1) - (size_t)(const void *)(strings[i].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("CDC_DRIFT_SPEED"
) + 1) - (size_t)(const void *)("CDC_DRIFT_SPEED") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (strings[i].str, "CDC_DRIFT_SPEED"
) : (__builtin_constant_p (strings[i].str) && ((size_t
)(const void *)((strings[i].str) + 1) - (size_t)(const void *
)(strings[i].str) == 1) && (__s1_len = strlen (strings
[i].str), __s1_len < 4) ? (__builtin_constant_p ("CDC_DRIFT_SPEED"
) && ((size_t)(const void *)(("CDC_DRIFT_SPEED") + 1)
- (size_t)(const void *)("CDC_DRIFT_SPEED") == 1) ? __builtin_strcmp
(strings[i].str, "CDC_DRIFT_SPEED") : (__extension__ ({ __const
unsigned char *__s2 = (__const unsigned char *) (__const char
*) ("CDC_DRIFT_SPEED"); register int __result = (((__const unsigned
char *) (__const char *) (strings[i].str))[0] - __s2[0]); if
(__s1_len > 0 && __result == 0) { __result = (((__const
unsigned char *) (__const char *) (strings[i].str))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((__const unsigned char *) (__const char *) (strings[i].str
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((__const unsigned char *) (__const char *) (strings
[i].str))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("CDC_DRIFT_SPEED") && ((size_t)(const void *)(("CDC_DRIFT_SPEED"
) + 1) - (size_t)(const void *)("CDC_DRIFT_SPEED") == 1) &&
(__s2_len = strlen ("CDC_DRIFT_SPEED"), __s2_len < 4) ? (
__builtin_constant_p (strings[i].str) && ((size_t)(const
void *)((strings[i].str) + 1) - (size_t)(const void *)(strings
[i].str) == 1) ? __builtin_strcmp (strings[i].str, "CDC_DRIFT_SPEED"
) : (__extension__ ({ __const unsigned char *__s1 = (__const unsigned
char *) (__const char *) (strings[i].str); register int __result
= __s1[0] - ((__const unsigned char *) (__const char *) ("CDC_DRIFT_SPEED"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("CDC_DRIFT_SPEED"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("CDC_DRIFT_SPEED"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("CDC_DRIFT_SPEED"
))[3]); } } __result; }))) : __builtin_strcmp (strings[i].str
, "CDC_DRIFT_SPEED")))); })) {
252 DRIFT_SPEED = values[i];
253 ncounter++;
254 }
255 if (!strcmp(strings[i].str,"CDC_TWO_HIT_RESOL")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[i].str) && __builtin_constant_p ("CDC_TWO_HIT_RESOL"
) && (__s1_len = strlen (strings[i].str), __s2_len = strlen
("CDC_TWO_HIT_RESOL"), (!((size_t)(const void *)((strings[i]
.str) + 1) - (size_t)(const void *)(strings[i].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("CDC_TWO_HIT_RESOL"
) + 1) - (size_t)(const void *)("CDC_TWO_HIT_RESOL") == 1) ||
__s2_len >= 4)) ? __builtin_strcmp (strings[i].str, "CDC_TWO_HIT_RESOL"
) : (__builtin_constant_p (strings[i].str) && ((size_t
)(const void *)((strings[i].str) + 1) - (size_t)(const void *
)(strings[i].str) == 1) && (__s1_len = strlen (strings
[i].str), __s1_len < 4) ? (__builtin_constant_p ("CDC_TWO_HIT_RESOL"
) && ((size_t)(const void *)(("CDC_TWO_HIT_RESOL") + 1
) - (size_t)(const void *)("CDC_TWO_HIT_RESOL") == 1) ? __builtin_strcmp
(strings[i].str, "CDC_TWO_HIT_RESOL") : (__extension__ ({ __const
unsigned char *__s2 = (__const unsigned char *) (__const char
*) ("CDC_TWO_HIT_RESOL"); register int __result = (((__const
unsigned char *) (__const char *) (strings[i].str))[0] - __s2
[0]); if (__s1_len > 0 && __result == 0) { __result
= (((__const unsigned char *) (__const char *) (strings[i].str
))[1] - __s2[1]); if (__s1_len > 1 && __result == 0
) { __result = (((__const unsigned char *) (__const char *) (
strings[i].str))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((__const unsigned char *) (__const
char *) (strings[i].str))[3] - __s2[3]); } } __result; }))) :
(__builtin_constant_p ("CDC_TWO_HIT_RESOL") && ((size_t
)(const void *)(("CDC_TWO_HIT_RESOL") + 1) - (size_t)(const void
*)("CDC_TWO_HIT_RESOL") == 1) && (__s2_len = strlen (
"CDC_TWO_HIT_RESOL"), __s2_len < 4) ? (__builtin_constant_p
(strings[i].str) && ((size_t)(const void *)((strings
[i].str) + 1) - (size_t)(const void *)(strings[i].str) == 1) ?
__builtin_strcmp (strings[i].str, "CDC_TWO_HIT_RESOL") : (__extension__
({ __const unsigned char *__s1 = (__const unsigned char *) (
__const char *) (strings[i].str); register int __result = __s1
[0] - ((__const unsigned char *) (__const char *) ("CDC_TWO_HIT_RESOL"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("CDC_TWO_HIT_RESOL"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("CDC_TWO_HIT_RESOL"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("CDC_TWO_HIT_RESOL"
))[3]); } } __result; }))) : __builtin_strcmp (strings[i].str
, "CDC_TWO_HIT_RESOL")))); })) {
256 TWO_HIT_RESOL = values[i];
257 ncounter++;
258 }
259 if (!strcmp(strings[i].str,"CDC_MAX_HITS")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[i].str) && __builtin_constant_p ("CDC_MAX_HITS"
) && (__s1_len = strlen (strings[i].str), __s2_len = strlen
("CDC_MAX_HITS"), (!((size_t)(const void *)((strings[i].str)
+ 1) - (size_t)(const void *)(strings[i].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("CDC_MAX_HITS"
) + 1) - (size_t)(const void *)("CDC_MAX_HITS") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (strings[i].str, "CDC_MAX_HITS"
) : (__builtin_constant_p (strings[i].str) && ((size_t
)(const void *)((strings[i].str) + 1) - (size_t)(const void *
)(strings[i].str) == 1) && (__s1_len = strlen (strings
[i].str), __s1_len < 4) ? (__builtin_constant_p ("CDC_MAX_HITS"
) && ((size_t)(const void *)(("CDC_MAX_HITS") + 1) - (
size_t)(const void *)("CDC_MAX_HITS") == 1) ? __builtin_strcmp
(strings[i].str, "CDC_MAX_HITS") : (__extension__ ({ __const
unsigned char *__s2 = (__const unsigned char *) (__const char
*) ("CDC_MAX_HITS"); register int __result = (((__const unsigned
char *) (__const char *) (strings[i].str))[0] - __s2[0]); if
(__s1_len > 0 && __result == 0) { __result = (((__const
unsigned char *) (__const char *) (strings[i].str))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((__const unsigned char *) (__const char *) (strings[i].str
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((__const unsigned char *) (__const char *) (strings
[i].str))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("CDC_MAX_HITS") && ((size_t)(const void *)(("CDC_MAX_HITS"
) + 1) - (size_t)(const void *)("CDC_MAX_HITS") == 1) &&
(__s2_len = strlen ("CDC_MAX_HITS"), __s2_len < 4) ? (__builtin_constant_p
(strings[i].str) && ((size_t)(const void *)((strings
[i].str) + 1) - (size_t)(const void *)(strings[i].str) == 1) ?
__builtin_strcmp (strings[i].str, "CDC_MAX_HITS") : (__extension__
({ __const unsigned char *__s1 = (__const unsigned char *) (
__const char *) (strings[i].str); register int __result = __s1
[0] - ((__const unsigned char *) (__const char *) ("CDC_MAX_HITS"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("CDC_MAX_HITS"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("CDC_MAX_HITS"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("CDC_MAX_HITS"
))[3]); } } __result; }))) : __builtin_strcmp (strings[i].str
, "CDC_MAX_HITS")))); })) {
260 MAX_HITS = (int)values[i];
261 ncounter++;
262 }
263 if (!strcmp(strings[i].str,"CDC_THRESH_KEV")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[i].str) && __builtin_constant_p ("CDC_THRESH_KEV"
) && (__s1_len = strlen (strings[i].str), __s2_len = strlen
("CDC_THRESH_KEV"), (!((size_t)(const void *)((strings[i].str
) + 1) - (size_t)(const void *)(strings[i].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("CDC_THRESH_KEV"
) + 1) - (size_t)(const void *)("CDC_THRESH_KEV") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (strings[i].str, "CDC_THRESH_KEV"
) : (__builtin_constant_p (strings[i].str) && ((size_t
)(const void *)((strings[i].str) + 1) - (size_t)(const void *
)(strings[i].str) == 1) && (__s1_len = strlen (strings
[i].str), __s1_len < 4) ? (__builtin_constant_p ("CDC_THRESH_KEV"
) && ((size_t)(const void *)(("CDC_THRESH_KEV") + 1) -
(size_t)(const void *)("CDC_THRESH_KEV") == 1) ? __builtin_strcmp
(strings[i].str, "CDC_THRESH_KEV") : (__extension__ ({ __const
unsigned char *__s2 = (__const unsigned char *) (__const char
*) ("CDC_THRESH_KEV"); register int __result = (((__const unsigned
char *) (__const char *) (strings[i].str))[0] - __s2[0]); if
(__s1_len > 0 && __result == 0) { __result = (((__const
unsigned char *) (__const char *) (strings[i].str))[1] - __s2
[1]); if (__s1_len > 1 && __result == 0) { __result
= (((__const unsigned char *) (__const char *) (strings[i].str
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((__const unsigned char *) (__const char *) (strings
[i].str))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("CDC_THRESH_KEV") && ((size_t)(const void *)(("CDC_THRESH_KEV"
) + 1) - (size_t)(const void *)("CDC_THRESH_KEV") == 1) &&
(__s2_len = strlen ("CDC_THRESH_KEV"), __s2_len < 4) ? (__builtin_constant_p
(strings[i].str) && ((size_t)(const void *)((strings
[i].str) + 1) - (size_t)(const void *)(strings[i].str) == 1) ?
__builtin_strcmp (strings[i].str, "CDC_THRESH_KEV") : (__extension__
({ __const unsigned char *__s1 = (__const unsigned char *) (
__const char *) (strings[i].str); register int __result = __s1
[0] - ((__const unsigned char *) (__const char *) ("CDC_THRESH_KEV"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("CDC_THRESH_KEV"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("CDC_THRESH_KEV"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("CDC_THRESH_KEV"
))[3]); } } __result; }))) : __builtin_strcmp (strings[i].str
, "CDC_THRESH_KEV")))); })) {
264 THRESH_KEV = values[i];
265 ncounter++;
266 }
267 }
268 if (ncounter==4) {
269 printf("CDC: ALL parameters loaded from Data Base\n");
270 } else if (ncounter<5) {
271 printf("CDC: NOT ALL necessary parameters found in Data Base %d out of 5\n",ncounter);
272 } else {
273 printf("CDC: SOME parameters found more than once in Data Base\n");
274 }
275 }
276 //
277 // Get drift table and scale factors from the database
278 //
279 // First check for non-zero field in the magnet bore
280 float x[3]={0.,0.,65};
281 float B[3];
282 gufld_db_(x,B);
283 if (fabs(B[2])>1e-3){
284 nvalues=78;
285 status=GetColumn("CDC/cdc_drift_table",&nvalues,cdc_drift_time,"t");
286 int k;
287 for (k=0;k<nvalues;k++){
288 cdc_drift_time[k]*=1000.; // Scale fron micro-secons to ns
289 cdc_drift_distance[k]=0.01*(float)k; // 100 micron increments;
290 }
291
292 nvalues=2;
293 status = GetConstants("CDC/cdc_drift_parms", &nvalues, values, strings);
294 for ( k=0;k<(int)nvalues;k++) {
295 if (!strcmp(strings[k].str,"bscale_par1")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[k].str) && __builtin_constant_p ("bscale_par1"
) && (__s1_len = strlen (strings[k].str), __s2_len = strlen
("bscale_par1"), (!((size_t)(const void *)((strings[k].str) +
1) - (size_t)(const void *)(strings[k].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("bscale_par1"
) + 1) - (size_t)(const void *)("bscale_par1") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (strings[k].str, "bscale_par1")
: (__builtin_constant_p (strings[k].str) && ((size_t
)(const void *)((strings[k].str) + 1) - (size_t)(const void *
)(strings[k].str) == 1) && (__s1_len = strlen (strings
[k].str), __s1_len < 4) ? (__builtin_constant_p ("bscale_par1"
) && ((size_t)(const void *)(("bscale_par1") + 1) - (
size_t)(const void *)("bscale_par1") == 1) ? __builtin_strcmp
(strings[k].str, "bscale_par1") : (__extension__ ({ __const unsigned
char *__s2 = (__const unsigned char *) (__const char *) ("bscale_par1"
); register int __result = (((__const unsigned char *) (__const
char *) (strings[k].str))[0] - __s2[0]); if (__s1_len > 0
&& __result == 0) { __result = (((__const unsigned char
*) (__const char *) (strings[k].str))[1] - __s2[1]); if (__s1_len
> 1 && __result == 0) { __result = (((__const unsigned
char *) (__const char *) (strings[k].str))[2] - __s2[2]); if
(__s1_len > 2 && __result == 0) __result = (((__const
unsigned char *) (__const char *) (strings[k].str))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p ("bscale_par1"
) && ((size_t)(const void *)(("bscale_par1") + 1) - (
size_t)(const void *)("bscale_par1") == 1) && (__s2_len
= strlen ("bscale_par1"), __s2_len < 4) ? (__builtin_constant_p
(strings[k].str) && ((size_t)(const void *)((strings
[k].str) + 1) - (size_t)(const void *)(strings[k].str) == 1) ?
__builtin_strcmp (strings[k].str, "bscale_par1") : (__extension__
({ __const unsigned char *__s1 = (__const unsigned char *) (
__const char *) (strings[k].str); register int __result = __s1
[0] - ((__const unsigned char *) (__const char *) ("bscale_par1"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("bscale_par1"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("bscale_par1"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("bscale_par1"
))[3]); } } __result; }))) : __builtin_strcmp (strings[k].str
, "bscale_par1")))); })) {
296 BSCALE_PAR1 = values[k];
297 }
298 if (!strcmp(strings[k].str,"bscale_par2")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(strings[k].str) && __builtin_constant_p ("bscale_par2"
) && (__s1_len = strlen (strings[k].str), __s2_len = strlen
("bscale_par2"), (!((size_t)(const void *)((strings[k].str) +
1) - (size_t)(const void *)(strings[k].str) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)(("bscale_par2"
) + 1) - (size_t)(const void *)("bscale_par2") == 1) || __s2_len
>= 4)) ? __builtin_strcmp (strings[k].str, "bscale_par2")
: (__builtin_constant_p (strings[k].str) && ((size_t
)(const void *)((strings[k].str) + 1) - (size_t)(const void *
)(strings[k].str) == 1) && (__s1_len = strlen (strings
[k].str), __s1_len < 4) ? (__builtin_constant_p ("bscale_par2"
) && ((size_t)(const void *)(("bscale_par2") + 1) - (
size_t)(const void *)("bscale_par2") == 1) ? __builtin_strcmp
(strings[k].str, "bscale_par2") : (__extension__ ({ __const unsigned
char *__s2 = (__const unsigned char *) (__const char *) ("bscale_par2"
); register int __result = (((__const unsigned char *) (__const
char *) (strings[k].str))[0] - __s2[0]); if (__s1_len > 0
&& __result == 0) { __result = (((__const unsigned char
*) (__const char *) (strings[k].str))[1] - __s2[1]); if (__s1_len
> 1 && __result == 0) { __result = (((__const unsigned
char *) (__const char *) (strings[k].str))[2] - __s2[2]); if
(__s1_len > 2 && __result == 0) __result = (((__const
unsigned char *) (__const char *) (strings[k].str))[3] - __s2
[3]); } } __result; }))) : (__builtin_constant_p ("bscale_par2"
) && ((size_t)(const void *)(("bscale_par2") + 1) - (
size_t)(const void *)("bscale_par2") == 1) && (__s2_len
= strlen ("bscale_par2"), __s2_len < 4) ? (__builtin_constant_p
(strings[k].str) && ((size_t)(const void *)((strings
[k].str) + 1) - (size_t)(const void *)(strings[k].str) == 1) ?
__builtin_strcmp (strings[k].str, "bscale_par2") : (__extension__
({ __const unsigned char *__s1 = (__const unsigned char *) (
__const char *) (strings[k].str); register int __result = __s1
[0] - ((__const unsigned char *) (__const char *) ("bscale_par2"
))[0]; if (__s2_len > 0 && __result == 0) { __result
= (__s1[1] - ((__const unsigned char *) (__const char *) ("bscale_par2"
))[1]); if (__s2_len > 1 && __result == 0) { __result
= (__s1[2] - ((__const unsigned char *) (__const char *) ("bscale_par2"
))[2]); if (__s2_len > 2 && __result == 0) __result
= (__s1[3] - ((__const unsigned char *) (__const char *) ("bscale_par2"
))[3]); } } __result; }))) : __builtin_strcmp (strings[k].str
, "bscale_par2")))); })) {
299 BSCALE_PAR2 = values[k];
300 }
301 }
302 }
303 else{
304 nvalues=78;
305 status=GetColumn("CDC/cdc_drift_table::NoBField",&nvalues,cdc_drift_time,"t");
Value stored to 'status' is never read
306 int k;
307 for (k=0;k<nvalues;k++){
308 cdc_drift_time[k]*=1000.; // Scale fron micro-secons to ns
309 cdc_drift_distance[k]=0.01*(float)k; // 100 micron increments;
310 }
311 }
312
313 initialized = 1;
314 }
315
316 x[0] = (xin[0] + xout[0])/2;
317 x[1] = (xin[1] + xout[1])/2;
318 x[2] = (xin[2] + xout[2])/2;
319 t = (xin[3] + xout[3])/2 * 1e9;
320 dx[0] = xin[0] - xout[0];
321 dx[1] = xin[1] - xout[1];
322 dx[2] = xin[2] - xout[2];
323 transformCoord(xin,"global",xinlocal,"local")transformcoord_(xin,"global",xinlocal,"local",strlen("global"
),strlen("local"));
324 transformCoord(xout,"global",xoutlocal,"local")transformcoord_(xout,"global",xoutlocal,"local",strlen("global"
),strlen("local"));
325
326 /*
327 xlocal[0] = (xinlocal[0] + xoutlocal[0])/2;
328 xlocal[1] = (xinlocal[1] + xoutlocal[1])/2;
329 xlocal[2] = (xinlocal[2] + xoutlocal[2])/2;
330 */
331
332 /* For particles that range out inside the active volume, the
333 * "out" time seems to be set to something enormously high.
334 * This screws up the hit. Check for this case here by looking
335 * at xout[3] and making sure it is less than 1 second. If it's
336 * not, then just use xin[3] for "t".
337 */
338 if (xout[3] > 1.0)
339 t = xin[3] * 1e9;
340
341 drin = sqrt(xinlocal[0]*xinlocal[0] + xinlocal[1]*xinlocal[1]);
342 drout = sqrt(xoutlocal[0]*xoutlocal[0] + xoutlocal[1]*xoutlocal[1]);
343
344 trackdir[0] =-xinlocal[0] + xoutlocal[0];
345 trackdir[1] =-xinlocal[1] + xoutlocal[1];
346 trackdir[2] =-xinlocal[2] + xoutlocal[2];
347 alpha=-(xinlocal[0]*trackdir[0]+xinlocal[1]*trackdir[1])
348 /(trackdir[0]*trackdir[0]+trackdir[1]*trackdir[1]);
349 xlocal[0]=xinlocal[0]+trackdir[0]*alpha;
350 xlocal[1]=xinlocal[1]+trackdir[1]*alpha;
351 xlocal[2]=xinlocal[2]+trackdir[2]*alpha;
352
353 // Deal with tracks exiting the ends of the straws
354 if (fabs(xlocal[2]) >= 75.45) {
355 float sign = (xoutlocal[2] > 0)? 1. : -1.;
356 int ring = getring_wrapper_();
357 if (ring <= 4 || (ring >= 13 && ring <= 16) || ring >= 25) {
358 alpha=(sign*75.45-xinlocal[2])/trackdir[2];
359 xlocal[0]=xinlocal[0]+trackdir[0]*alpha;
360 xlocal[1]=xinlocal[1]+trackdir[1]*alpha;
361 xlocal[2]=sign*75.45;
362 }
363 else if (fabs(xlocal[2]) >= 75.575) {
364 alpha=(sign*75.575-xinlocal[2])/trackdir[2];
365 xlocal[0]=xinlocal[0]+trackdir[0]*alpha;
366 xlocal[1]=xinlocal[1]+trackdir[1]*alpha;
367 xlocal[2]=sign*75.575;
368 }
369 }
370
371 /* This will get called when the particle actually passes through
372 * the wire volume itself. For these cases, we should set the
373 * location of the hit to be the point on the wire itself. Do
374 * determine if this is what is happening, we check drout to
375 * see if it is very close to the wire and drin to see if it is
376 * close to the tube.
377 *
378 * For the other case, when drin is close to the wire, we assume
379 * it is because it is emerging from the wire volume and
380 * automatically ignore those hits by returning immediately.
381 */
382 if (drin < 0.0050)
383 return; /* entering straw within 50 microns of wire. ignore */
384
385 if ((drin > (STRAW_RADIUS-0.0200) && drout<0.0050) ||
386 (drin < 0.274 && drin > 0.234 && drout<0.0050))
387 {
388 /* Either we entered within 200 microns of the straw tube and left
389 * within 50 microns of the wire or we entered the stub region near the
390 * donuts at either end of the straw (the inner radius of the feedthrough
391 * region is 0.254 cm) and passed near the wire. Assume the track passed
392 * through the wire volume.
393 */
394
395 x[0] = xout[0];
396 x[1] = xout[1];
397 x[2] = xout[2];
398 t = xout[3] * 1e9;
399 xlocal[0] = xoutlocal[0];
400 xlocal[1] = xoutlocal[1];
401 xlocal[2] = xoutlocal[2];
402
403 /* For dx, we will just assume it is twice the distance from
404 * the straw to wire.
405 */
406 dx[0] *= 2.0;
407 dx[1] *= 2.0;
408 dx[2] *= 2.0;
409
410 /* We will approximate the energy loss in the straw to be twice the
411 energy loss in the first half of the straw */
412 dEsum *= 2.0;
413 }
414
415 /* Distance of hit from center of wire */
416 dradius = sqrt(xlocal[0]*xlocal[0] + xlocal[1]*xlocal[1]);
417
418 /* Calculate dE/dx */
419
420 dr = sqrt(dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]);
421 if (dr > 1e-3)
422 {
423 dEdx = dEsum/dr;
424 }
425 else
426 {
427 dEdx = 0;
428 }
429
430 /* post the hit to the truth tree */
431
432 if (history == 0)
433 {
434 int mark = (1<<30) + pointCount;
435 void** twig = getTwig(&centralDCTree, mark);
436 if (*twig == 0)
437 {
438 s_CentralDC_t* cdc = *twig = make_s_CentralDC();
439 s_CdcTruthPoints_t* points = make_s_CdcTruthPoints(1);
440 int a = thisInputEvent->physicsEvents->in[0].reactions->in[0].vertices->in[0].products->mult;
441 points->in[0].primary = (stack <= a);
442 points->in[0].track = track;
443 points->in[0].t = t;
444 points->in[0].z = x[2];
445 points->in[0].r = sqrt(x[0]*x[0] + x[1]*x[1]);
446 points->in[0].phi = atan2(x[1],x[0]);
447 points->in[0].dradius = dradius;
448 points->in[0].px = pin[0]*pin[4];
449 points->in[0].py = pin[1]*pin[4];
450 points->in[0].pz = pin[2]*pin[4];
451 points->in[0].dEdx = dEdx;
452 points->in[0].ptype = ipart;
453 points->in[0].trackID = make_s_TrackID();
454 points->in[0].trackID->itrack = gidGetId(track);
455 points->mult = 1;
456 cdc->cdcTruthPoints = points;
457 pointCount++;
458 }
459 }
460
461 /* post the hit to the hits tree, mark sector as hit */
462
463 if (dEsum > 0)
464 {
465 s_CdcStrawTruthHits_t* hits;
466
467 int layer = getlayer_wrapper_();
468 int ring = getring_wrapper_();
469 int sector = getsector_wrapper_();
470
471 if (layer == 0) /* in a straw */
472 {
473 int mark = (ring<<20) + sector;
474 void** twig = getTwig(&centralDCTree, mark);
475
476 if (*twig == 0)
477 {
478 s_CentralDC_t* cdc = *twig = make_s_CentralDC();
479 s_CdcStraws_t* straws = make_s_CdcStraws(1);
480 straws->mult = 1;
481 straws->in[0].ring = ring;
482 straws->in[0].straw = sector;
483 straws->in[0].cdcStrawTruthHits = hits = make_s_CdcStrawTruthHits(MAX_HITS);
484 cdc->cdcStraws = straws;
485 strawCount++;
486 }
487 else
488 {
489 s_CentralDC_t* cdc = (s_CentralDC_t*) *twig;
490 hits = cdc->cdcStraws->in[0].cdcStrawTruthHits;
491 }
492
493
494 /* Simulate number of primary ion pairs*/
495 /* The total number of ion pairs depends on the energy deposition
496 and the effective average energy to produce a pair, w_eff.
497 On average for each primary ion pair produced there are n_s_per_p
498 secondary ion pairs produced.
499 */
500 int one=1;
501 // Average number of secondary ion pairs for 50/50 Ar/CO2 mixture
502 float n_s_per_p=1.94;
503 //Average energy needed to produce an ion pair for 50/50 mixture
504 float w_eff=29.5e-9; // GeV
505 // Average number of primary ion pairs
506 float n_p_mean = dEsum/w_eff/(1.+n_s_per_p);
507 int n_p; // number of primary ion pairs
508 gpoiss_(&n_p_mean,&n_p,&one);
509
510 if (controlparams_.driftclusters==0) {
511 AddCDCCluster(hits,ipart,track,n_p,t,xlocal);
512 }
513 else{
514 // Loop over the number of primary ion pairs
515 int n;
516 for (n=0; n < n_p; n++) {
517 // Generate a cluster at a random position along the path within the
518 // straw
519 int one=2;
520 float rndno[1];
521 grndm_(rndno,&one);
522 xlocal[0]=xinlocal[0]+trackdir[0]*rndno[0];
523 xlocal[1]=xinlocal[1]+trackdir[1]*rndno[0];
524 xlocal[2]=xinlocal[2]+trackdir[2]*rndno[0];
525
526 AddCDCCluster(hits,ipart,track,n_p,t,xlocal);
527 }
528 }
529 }
530 }
531}
532
533/* entry points from fortran */
534
535void hitcentraldc_(float* xin, float* xout,
536 float* pin, float* pout, float* dEsum,
537 int* track, int* stack, int* history, int* ipart)
538{
539 hitCentralDC(xin,xout,pin,pout,*dEsum,*track,*stack,*history, *ipart);
540}
541
542
543/* pick and package the hits for shipping */
544
545s_CentralDC_t* pickCentralDC ()
546{
547 s_CentralDC_t* box;
548 s_CentralDC_t* item;
549
550 if ((strawCount == 0) && (stripCount == 0) && (pointCount == 0))
551 {
552 return HDDM_NULL(void*)&hddm_s_nullTarget;
553 }
554
555 box = make_s_CentralDC();
556 box->cdcStraws = make_s_CdcStraws(strawCount);
557 box->cdcTruthPoints = make_s_CdcTruthPoints(pointCount);
558
559 while ((item = (s_CentralDC_t*) pickTwig(&centralDCTree)))
560 {
561 s_CdcStraws_t* straws = item->cdcStraws;
562 int straw;
563
564 s_CdcTruthPoints_t* points = item->cdcTruthPoints;
565 int point;
566 for (straw=0; straw < straws->mult; ++straw)
567 {
568 int m = box->cdcStraws->mult;
569
570 s_CdcStrawTruthHits_t* hits = straws->in[straw].cdcStrawTruthHits;
571
572 // Sort the clusters by time
573 qsort(hits->in,hits->mult,sizeof(s_CdcStrawTruthHit_t),(compfn)cdc_cluster_sort);
574
575 /* compress out the hits below threshold */
576 int i,iok=0;
577
578 if (controlparams_.driftclusters == 0) {
579 for (iok=i=0; i < hits->mult; i++)
580 {
581 if (hits->in[i].q >0.)
582 {
583 if (iok < i)
584 {
585 hits->in[iok] = hits->in[i];
586 }
587 ++iok;
588 }
589 }
590 }
591 else{
592
593 // Temporary histogram in 1 ns bins to store waveform data
594 int num_samples=(int)CDC_TIME_WINDOW;
595 float *samples=(float *)malloc(num_samples*sizeof(float));
596 for (i=0;i<num_samples;i++) {
597 samples[i]=cdc_wire_signal((float)i,hits);
598 //printf("%f %f\n",(float)i,samples[i]);
599 }
600
601 int returned_to_baseline=0;
602 float q=0.;
603 int FADC_BIN_SIZE=1;
604 for (i=0; i<num_samples; i+=FADC_BIN_SIZE) {
605 if (samples[i] >= THRESH_MV) {
606 if (returned_to_baseline == 0) {
607 hits->in[iok].itrack = hits->in[0].itrack;
608 hits->in[iok].ptype = hits->in[0].ptype;
609 hits->in[iok].t=(float) i;
610 returned_to_baseline = 1;
611 iok++;
612 }
613 q += (float)FADC_BIN_SIZE*samples[i];
614 }
615 if (returned_to_baseline && (samples[i] < THRESH_MV)) {
616 returned_to_baseline = 0;
617 if (iok > 0 && q > 0.) {
618 hits->in[iok-1].q=q;
619 q=0.;
620 }
621 //break;
622 }
623 }
624 if (q > 0) {
625 hits->in[iok-1].q = q;
626 }
627 free(samples);
628 }
629
630 if (iok)
631 {
632 hits->mult = iok;
633 box->cdcStraws->in[m] = straws->in[straw];
634 box->cdcStraws->mult++;
635 }
636 else if (hits != HDDM_NULL(void*)&hddm_s_nullTarget)
637 {
638 FREE(hits)free(hits);
639 }
640 }
641 if (straws != HDDM_NULL(void*)&hddm_s_nullTarget)
642 {
643 FREE(straws)free(straws);
644 }
645
646 for (point=0; point < points->mult; ++point)
647 {
648 int m = box->cdcTruthPoints->mult++;
649 box->cdcTruthPoints->in[m] = points->in[point];
650 }
651 if (points != HDDM_NULL(void*)&hddm_s_nullTarget)
652 {
653 FREE(points)free(points);
654 }
655 FREE(item)free(item);
656 }
657
658 strawCount = stripCount = pointCount = 0;
659
660 if ((box->cdcStraws != HDDM_NULL(void*)&hddm_s_nullTarget) &&
661 (box->cdcStraws->mult == 0))
662 {
663 FREE(box->cdcStraws)free(box->cdcStraws);
664 box->cdcStraws = HDDM_NULL(void*)&hddm_s_nullTarget;
665 }
666 if ((box->cdcTruthPoints != HDDM_NULL(void*)&hddm_s_nullTarget) &&
667 (box->cdcTruthPoints->mult == 0))
668 {
669 FREE(box->cdcTruthPoints)free(box->cdcTruthPoints);
670 box->cdcTruthPoints = HDDM_NULL(void*)&hddm_s_nullTarget;
671 }
672 if ((box->cdcStraws->mult == 0) &&
673 (box->cdcTruthPoints->mult == 0))
674 {
675 FREE(box)free(box);
676 box = HDDM_NULL(void*)&hddm_s_nullTarget;
677
678 }
679 return box;
680}