File: | programs/Simulation/HDGeant/hitStart.c |
Location: | line 133, column 7 |
Description: | Value stored to 'dEdx' is never read |
1 | /* |
2 | * hitStart - registers hits for Start counter |
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 hitStartCntr |
11 | * |
12 | * Programmer's Notes: |
13 | * ------------------- |
14 | * 1) In applying the attenuation to light propagating down to the end |
15 | * of the counters, there has to be some point where the attenuation |
16 | * factor is 1. I chose it to be the midplane, so that in the middle |
17 | * of the counters the attenuation factor is 1. |
18 | * 2) In applying the propagation delay to light propagating down to the |
19 | * end of the counters, there has to be some point where the timing |
20 | * offset is 0. I chose it to be the midplane, so that for hits in |
21 | * the middle of the counter the t values measure time-of-flight from |
22 | * the t=0 of the event. |
23 | */ |
24 | |
25 | #include <stdlib.h> |
26 | #include <stdio.h> |
27 | #include <math.h> |
28 | |
29 | #include <hddm_s.h> |
30 | #include <geant3.h> |
31 | #include <bintree.h> |
32 | #include "calibDB.h" |
33 | extern s_HDDM_t* thisInputEvent; |
34 | |
35 | static float ATTEN_LENGTH = 150.; |
36 | static float C_EFFECTIVE = 15.; |
37 | static float TWO_HIT_RESOL = 25.; |
38 | static int MAX_HITS = 100; |
39 | static float THRESH_MEV = 0.150; |
40 | static float LIGHT_GUIDE = 0.; |
41 | static float ANGLE_COR = 1.038; |
42 | static float BENT_REGION = 39.465; |
43 | |
44 | binTree_t* startCntrTree = 0; |
45 | static int paddleCount = 0; |
46 | static int pointCount = 0; |
47 | static int initialized = 0; |
48 | |
49 | |
50 | /* register hits during tracking (from gustep) */ |
51 | |
52 | void hitStartCntr (float xin[4], float xout[4], |
53 | float pin[5], float pout[5], float dEsum, |
54 | int track, int stack, int history, int ipart) |
55 | { |
56 | float x[3], t; |
57 | float dx[3], dr; |
58 | float dEdx; |
59 | float xlocal[3]; |
60 | float xvrtx[3]; |
61 | |
62 | if (!initialized) { |
63 | |
64 | mystr_t strings[50]; |
65 | float values[50]; |
66 | int nvalues = 50; |
67 | int status = GetConstants("START_COUNTER/start_parms", &nvalues, values, strings); |
68 | |
69 | if (!status) { |
70 | int ncounter = 0; |
71 | int i; |
72 | for ( i=0;i<(int)nvalues;i++){ |
73 | //printf("%d %s \n",i,strings[i].str); |
74 | if (!strcmp(strings[i].str,"START_ATTEN_LENGTH")) { |
75 | ATTEN_LENGTH = values[i]; |
76 | ncounter++; |
77 | } |
78 | if (!strcmp(strings[i].str,"START_C_EFFECTIVE")) { |
79 | C_EFFECTIVE = values[i]; |
80 | ncounter++; |
81 | } |
82 | if (!strcmp(strings[i].str,"START_TWO_HIT_RESOL")) { |
83 | TWO_HIT_RESOL = values[i]; |
84 | ncounter++; |
85 | } |
86 | if (!strcmp(strings[i].str,"START_MAX_HITS")) { |
87 | MAX_HITS = (int)values[i]; |
88 | ncounter++; |
89 | } |
90 | if (!strcmp(strings[i].str,"START_THRESH_MEV")) { |
91 | TWO_HIT_RESOL = values[i]; |
92 | ncounter++; |
93 | } |
94 | if (!strcmp(strings[i].str,"START_LIGHT_GUIDE")) { |
95 | LIGHT_GUIDE = values[i]; |
96 | ncounter++; |
97 | } |
98 | if (!strcmp(strings[i].str,"START_ANGLE_COR")) { |
99 | ANGLE_COR = values[i]; |
100 | ncounter++; |
101 | } |
102 | if (!strcmp(strings[i].str,"START_BENT_REGION")) { |
103 | BENT_REGION = values[i]; |
104 | ncounter++; |
105 | } |
106 | } |
107 | if (ncounter==8){ |
108 | printf("START: ALL parameters loaded from Data Base\n"); |
109 | } else if (ncounter<8){ |
110 | printf("START: NOT ALL necessary parameters found in Data Base %d out of 8\n",ncounter); |
111 | } else { |
112 | printf("START: SOME parameters found more than once in Data Base\n"); |
113 | } |
114 | } |
115 | initialized = 1; |
116 | } |
117 | |
118 | x[0] = (xin[0] + xout[0])/2; |
119 | x[1] = (xin[1] + xout[1])/2; |
120 | x[2] = (xin[2] + xout[2])/2; |
121 | t = (xin[3] + xout[3])/2 * 1e9; |
122 | transformCoord(x,"global",xlocal,"local")transformcoord_(x,"global",xlocal,"local",strlen("global"),strlen ("local")); |
123 | dx[0] = xin[0] - xout[0]; |
124 | dx[1] = xin[1] - xout[1]; |
125 | dx[2] = xin[2] - xout[2]; |
126 | dr = sqrt(dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]); |
127 | if (dr > 1e-3) |
128 | { |
129 | dEdx = dEsum/dr; |
130 | } |
131 | else |
132 | { |
133 | dEdx = 0; |
Value stored to 'dEdx' is never read | |
134 | } |
135 | |
136 | float dbent = 0.0; |
137 | float dpath = 0.0; |
138 | if(xlocal[2] >= BENT_REGION){ |
139 | dbent = ( xlocal[2] - BENT_REGION )*ANGLE_COR; |
140 | dpath =BENT_REGION + dbent; |
141 | } else { |
142 | dpath = xlocal[2]; |
143 | } |
144 | |
145 | float dEcorr = dEsum * exp(-dpath/ATTEN_LENGTH); |
146 | float tcorr = t + dpath/C_EFFECTIVE; |
147 | |
148 | |
149 | // printf("x_gl, z_gl, x_l, z_l %f %f %f\n", |
150 | // xin[0],xin[1],xin[2]); |
151 | |
152 | // printf("x_gl, z_gl, x_l, z_l %f %f %f %f %f %f %f\n", |
153 | // x[0],x[1],x[2], xlocal[0],xlocal[1],xlocal[2],dpath); |
154 | |
155 | |
156 | /* post the hit to the truth tree */ |
157 | |
158 | if (history == 0) |
159 | { |
160 | int mark = (1<<30) + pointCount; |
161 | void** twig = getTwig(&startCntrTree, mark); |
162 | if (*twig == 0) |
163 | { |
164 | s_StartCntr_t* stc = *twig = make_s_StartCntr(); |
165 | s_StcTruthPoints_t* points = make_s_StcTruthPoints(1); |
166 | stc->stcTruthPoints = points; |
167 | int a = thisInputEvent->physicsEvents->in[0].reactions->in[0].vertices->in[0].products->mult; |
168 | points->in[0].primary = (stack <= a); |
169 | points->in[0].track = track; |
170 | points->in[0].t = t; |
171 | points->in[0].z = x[2]; |
172 | points->in[0].r = sqrt(x[0]*x[0]+x[1]*x[1]); |
173 | points->in[0].phi = atan2(x[1],x[0]); |
174 | points->in[0].px = pin[0]*pin[4]; |
175 | points->in[0].py = pin[1]*pin[4]; |
176 | points->in[0].pz = pin[2]*pin[4]; |
177 | points->in[0].E = pin[3]; |
178 | points->in[0].dEdx = dEcorr; |
179 | points->in[0].ptype = ipart; |
180 | points->in[0].sector = getsector_(); |
181 | points->mult = 1; |
182 | pointCount++; |
183 | |
184 | |
185 | } |
186 | } |
187 | |
188 | /* post the hit to the hits tree, mark sector as hit */ |
189 | |
190 | // if( (ipart==8) && (x[2]<90.)){ |
191 | // printf("x_gl, z_gl, x_l, z_l %f %f %f %f %f %f\n", |
192 | // x[0],x[1],x[2], xlocal[0],xlocal[1],xlocal[2]); |
193 | // } |
194 | |
195 | |
196 | if (dEsum > 0) |
197 | { |
198 | int nhit; |
199 | s_StcTruthHits_t* hits; |
200 | int sector = getsector_(); |
201 | float phim = atan2(xvrtx[1],xvrtx[0]); |
202 | |
203 | |
204 | |
205 | // printf("x_gl, z_gl, x_l, z_l %f %f %f %f %f %f\n", |
206 | // x[0],x[1],x[2], xlocal[0],xlocal[1],xlocal[2]); |
207 | |
208 | |
209 | |
210 | // float dpath = xlocal[2]+(10.2-xlocal[0])*0.4; |
211 | // float tcorr = t + dpath/C_EFFECTIVE; |
212 | // float dEcorr = dEsum * exp(-dpath/ATTEN_LENGTH); |
213 | int mark = sector; |
214 | void** twig = getTwig(&startCntrTree, mark); |
215 | if (*twig == 0) |
216 | { |
217 | s_StartCntr_t* stc = *twig = make_s_StartCntr(); |
218 | s_StcPaddles_t* paddles = make_s_StcPaddles(1); |
219 | paddles->mult = 1; |
220 | paddles->in[0].sector = sector; |
221 | paddles->in[0].stcTruthHits = hits = make_s_StcTruthHits(MAX_HITS); |
222 | stc->stcPaddles = paddles; |
223 | paddleCount++; |
224 | } |
225 | else |
226 | { |
227 | s_StartCntr_t* stc = *twig; |
228 | hits = stc->stcPaddles->in[0].stcTruthHits; |
229 | } |
230 | |
231 | for (nhit = 0; nhit < hits->mult; nhit++) |
232 | { |
233 | if (fabs(hits->in[nhit].t - tcorr) < TWO_HIT_RESOL) |
234 | { |
235 | break; |
236 | } |
237 | } |
238 | if (nhit < hits->mult) /* merge with former hit */ |
239 | { |
240 | hits->in[nhit].t = |
241 | (hits->in[nhit].t * hits->in[nhit].dE + tcorr * dEcorr) |
242 | / (hits->in[nhit].dE += dEcorr); |
243 | } |
244 | else if (nhit < MAX_HITS) /* create new hit */ |
245 | { |
246 | hits->in[nhit].t = tcorr ; |
247 | hits->in[nhit].dE = dEcorr; |
248 | hits->mult++; |
249 | } |
250 | else |
251 | { |
252 | fprintf(stderrstderr,"HDGeant error in hitStart: "); |
253 | fprintf(stderrstderr,"max hit count %d exceeded, truncating!\n",MAX_HITS); |
254 | exit(2); |
255 | } |
256 | } |
257 | } |
258 | |
259 | /* entry point from fortran */ |
260 | |
261 | void hitstartcntr_(float* xin, float* xout, |
262 | float* pin, float* pout, float* dEsum, |
263 | int* track, int* stack, int* history, int* ipart) |
264 | { |
265 | hitStartCntr(xin,xout,pin,pout,*dEsum,*track,*stack,*history,*ipart); |
266 | } |
267 | |
268 | |
269 | /* pick and package the hits for shipping */ |
270 | |
271 | s_StartCntr_t* pickStartCntr () |
272 | { |
273 | s_StartCntr_t* box; |
274 | s_StartCntr_t* item; |
275 | |
276 | if ((paddleCount == 0) && (pointCount == 0)) |
277 | { |
278 | return HDDM_NULL(void*)&hddm_s_nullTarget; |
279 | } |
280 | |
281 | box = make_s_StartCntr(); |
282 | box->stcPaddles = make_s_StcPaddles(paddleCount); |
283 | box->stcTruthPoints = make_s_StcTruthPoints(pointCount); |
284 | while (item = (s_StartCntr_t*) pickTwig(&startCntrTree)) |
285 | { |
286 | s_StcPaddles_t* paddles = item->stcPaddles; |
287 | int paddle; |
288 | s_StcTruthPoints_t* points = item->stcTruthPoints; |
289 | int point; |
290 | |
291 | for (paddle=0; paddle < paddles->mult; ++paddle) |
292 | { |
293 | int m = box->stcPaddles->mult; |
294 | |
295 | s_StcTruthHits_t* hits = paddles->in[paddle].stcTruthHits; |
296 | |
297 | /* compress out the hits below threshold */ |
298 | int i,iok; |
299 | for (iok=i=0; i < hits->mult; i++) |
300 | { |
301 | if (hits->in[i].dE >= THRESH_MEV/1e3) |
302 | { |
303 | if (iok < i) |
304 | { |
305 | hits->in[iok] = hits->in[i]; |
306 | } |
307 | ++iok; |
308 | } |
309 | } |
310 | if (iok) |
311 | { |
312 | hits->mult = iok; |
313 | box->stcPaddles->in[m] = paddles->in[paddle]; |
314 | box->stcPaddles->mult++; |
315 | } |
316 | else if (hits != HDDM_NULL(void*)&hddm_s_nullTarget) |
317 | { |
318 | FREE(hits)free(hits); |
319 | } |
320 | } |
321 | if (paddles != HDDM_NULL(void*)&hddm_s_nullTarget) |
322 | { |
323 | FREE(paddles)free(paddles); |
324 | } |
325 | |
326 | for (point=0; point < points->mult; ++point) |
327 | { |
328 | int m = box->stcTruthPoints->mult++; |
329 | box->stcTruthPoints->in[m] = item->stcTruthPoints->in[point]; |
330 | } |
331 | if (points != HDDM_NULL(void*)&hddm_s_nullTarget) |
332 | { |
333 | FREE(points)free(points); |
334 | } |
335 | FREE(item)free(item); |
336 | } |
337 | |
338 | paddleCount = pointCount = 0; |
339 | |
340 | if ((box->stcPaddles != HDDM_NULL(void*)&hddm_s_nullTarget) && |
341 | (box->stcPaddles->mult == 0)) |
342 | { |
343 | FREE(box->stcPaddles)free(box->stcPaddles); |
344 | box->stcPaddles = HDDM_NULL(void*)&hddm_s_nullTarget; |
345 | } |
346 | if ((box->stcTruthPoints != HDDM_NULL(void*)&hddm_s_nullTarget) && |
347 | (box->stcTruthPoints->mult == 0)) |
348 | { |
349 | FREE(box->stcTruthPoints)free(box->stcTruthPoints); |
350 | box->stcTruthPoints = HDDM_NULL(void*)&hddm_s_nullTarget; |
351 | } |
352 | if ((box->stcPaddles->mult == 0) && |
353 | (box->stcTruthPoints->mult == 0)) |
354 | { |
355 | FREE(box)free(box); |
356 | box = HDDM_NULL(void*)&hddm_s_nullTarget; |
357 | } |
358 | return box; |
359 | } |