File: | programs/Simulation/HDGeant/hitCDC.c |
Location: | line 122, column 10 |
Description: | Value stored to 'd3' during its initialization is never read |
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" |
23 | extern s_HDDM_t* thisInputEvent; |
24 | |
25 | typedef struct { |
26 | int writeenohits; |
27 | int showersincol; |
28 | int driftclusters; |
29 | } controlparams_t; |
30 | |
31 | extern controlparams_t controlparams_; |
32 | |
33 | void 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 |
36 | static float DRIFT_SPEED = 0.0055; |
37 | static float TWO_HIT_RESOL = 25.; |
38 | static int MAX_HITS = 1000; |
39 | static float THRESH_KEV = 1.; |
40 | static float THRESH_MV = 1.; |
41 | static float STRAW_RADIUS = 0.776; |
42 | static float CDC_TIME_WINDOW = 1000.0; //time window for accepting CDC hits, ns |
43 | static float ELECTRON_CHARGE = 1.6022e-4; /* fC */ |
44 | static float GAS_GAIN = 1e5; |
45 | |
46 | binTree_t* centralDCTree = 0; |
47 | static int strawCount = 0; |
48 | static int pointCount = 0; |
49 | static int stripCount = 0; |
50 | static int initialized = 0; |
51 | static float cdc_drift_time[78]; |
52 | static float cdc_drift_distance[78]; |
53 | static float BSCALE_PAR1=0.; |
54 | static float BSCALE_PAR2=0.; |
55 | |
56 | static int itrack; |
57 | |
58 | /* void GetDOCA(int ipart, float x[3], float p[5], float doca[3]); disabled 6/24/2009 */ |
59 | |
60 | typedef int (*compfn)(const void*, const void*); |
61 | extern void polint(float *xa, float *ya,int n,float x, float *y,float *dy); |
62 | |
63 | // Sort function for sorting clusters |
64 | int 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 |
76 | double 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 |
93 | double 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 | |
106 | void 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 | |
219 | void 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(¢ralDCTree, 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(¢ralDCTree, 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 | |
540 | void 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 | |
550 | s_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(¢ralDCTree))) |
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 | } |