Bug Summary

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