Bug Summary

File:programs/Simulation/HDGeant/hitFTOF.c
Location:line 275, column 7
Description:Value stored to 'soMCHits' is never read

Annotated Source Code

1/*
2 * hitFTOF - registers hits for forward Time-Of-Flight
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: -B. Zihlmann June 19. 2007
10 * add hit position to north and south hit structure
11 * set THRESH_MEV to zero to NOT concatenate hits.
12 *
13 * changes: Wed Jun 20 13:19:56 EDT 2007 B. Zihlmann
14 * add ipart to the function hitforwardTOF
15 *
16 * Programmer's Notes:
17 * -------------------
18 * 1) In applying the attenuation to light propagating down to both ends
19 * of the counters, there has to be some point where the attenuation
20 * factor is 1. I chose it to be the midplane, so that in the middle
21 * of the counter both ends see the unattenuated dE values. Closer to
22 * either end, that end has a larger dE value and the opposite end a
23 * lower dE value than the actual deposition.
24 * 2) In applying the propagation delay to light propagating down to the
25 * ends of the counters, there has to be some point where the timing
26 * offset is 0. I chose it to be the midplane, so that for hits in
27 * the middle of the counter the t values measure time-of-flight from
28 * the t=0 of the event. For hits closer to one end, that end sees
29 * a t value smaller than its true time-of-flight, and the other end
30 * sees a value correspondingly larger. The average is the true tof.
31 */
32
33#include <stdlib.h>
34#include <stdio.h>
35#include <math.h>
36
37#include <hddm_s.h>
38#include <geant3.h>
39#include <bintree.h>
40
41#include "calibDB.h"
42
43// plastic scintillator specific constants
44static float ATTEN_LENGTH = 150 ;
45static float C_EFFECTIVE = 15.0 ;
46static float BAR_LENGTH = 252.0; // length of the bar
47
48// kinematic constants
49static float TWO_HIT_RESOL = 25. ;// separation time between two different hits
50
51static float THRESH_MEV = 0. ; // do not through away any hits, one can do that later
52
53// maximum particle tracks per counter
54static int MAX_HITS = 25; // was 100 changed to 25
55
56// maximum MC hits per paddle
57static int MAX_PAD_HITS = 25 ;
58
59// top level pointer of FTOF hit tree
60binTree_t* forwardTOFTree = 0;
61
62static int counterCount = 0;
63static int pointCount = 0;
64static int initialized = 0;
65
66
67/* register hits during tracking (from gustep) */
68// track is ITRA from GEANT
69// stack is ISTAK from GEANT
70// history is ISTORY from GEANT User flag for current track history (reset to 0 in GLTRAC)
71
72void hitForwardTOF (float xin[4], float xout[4],
73 float pin[5], float pout[5], float dEsum,
74 int track, int stack, int history, int ipart) {
75 float x[3], t;
76 float dx[3], dr;
77 float dEdx;
78 float xlocal[3];
79 float xftof[3];
80 float zeroHat[] = {0,0,0};
81
82 if (!initialized) {
83
84
85 mystr_t strings[50];
86 float values[50];
87 int nvalues = 50;
88 int status = GetConstants("TOF/tof_parms", &nvalues, values, strings);
89
90 if (!status) {
91
92 int ncounter = 0;
93 int i;
94 for ( i=0;i<(int)nvalues;i++){
95 //printf("%d %s \n",i,strings[i].str);
96 if (!strcmp(strings[i].str,"TOF_ATTEN_LENGTH")) {
97 ATTEN_LENGTH = values[i];
98 ncounter++;
99 }
100 if (!strcmp(strings[i].str,"TOF_C_EFFECTIVE")) {
101 C_EFFECTIVE = values[i];
102 ncounter++;
103 }
104 if (!strcmp(strings[i].str,"TOF_PADDLE_LENGTH")) {
105 BAR_LENGTH = values[i];
106 ncounter++;
107 }
108 if (!strcmp(strings[i].str,"TOF_TWO_HIT_RESOL")) {
109 TWO_HIT_RESOL = values[i];
110 ncounter++;
111 }
112 if (!strcmp(strings[i].str,"TOF_THRESH_MEV")) {
113 THRESH_MEV = values[i];
114 ncounter++;
115 }
116 if (!strcmp(strings[i].str,"TOF_MAX_HITS")){
117 MAX_HITS = (int)values[i];
118 ncounter++;
119 }
120 if (!strcmp(strings[i].str,"TOF_MAX_PAD_HITS")) {
121 MAX_PAD_HITS = (int)values[i];
122 ncounter++;
123 }
124 }
125 if (ncounter==7){
126 printf("TOF: ALL parameters loaded from Data Base\n");
127 } else if (ncounter<7){
128 printf("TOF: NOT ALL necessary parameters found in Data Base %d out of 7\n",ncounter);
129 } else {
130 printf("TOF: SOME parameters found more than once in Data Base\n");
131 }
132
133 }
134 initialized = 1;
135
136 }
137
138
139
140 // getplane is coded in
141 // src/programs/Simulation/HDGeant/hddsGeant3.F
142 // this file is automatically generated from the geometry file
143 // written in xml format
144 // NOTE: there are three files hddsGeant3.F with the same name in
145 // the source code tree namely
146 // 1) src/programs/Utilities/geantbfield2root/hddsGeant3.F
147 // 2) src/programs/Simulation/HDGeant/hddsGeant3.F
148 // 3) src/programs/Simulation/hdds/hddsGeant3.F
149 //
150 // while 2) and 3) are identical 1) is a part of 2) and 3)
151 int plane = getplane_();
152
153 // calculate mean location of track and mean time in [ns] units
154 // the units of xin xout and x are in [cm]
155 x[0] = (xin[0] + xout[0])/2;
156 x[1] = (xin[1] + xout[1])/2;
157 x[2] = (xin[2] + xout[2])/2;
158 t = (xin[3] + xout[3])/2 * 1e9;
159
160 // tranform the the global x coordinate into the local coordinate of the top_volume FTOF
161 // defined in the geometry file src/programs/Simulation/hdds/ForwardTOF_HDDS.xml
162 // the function transform Coord is defined in src/programs/Simulation/HDGeant/hitutil/hitutil.F
163 transformCoord(x,"global",xlocal,"FTOF")transformcoord_(x,"global",xlocal,"FTOF",strlen("global"),strlen
("FTOF"));
164 transformCoord(zeroHat,"local",xftof,"FTOF")transformcoord_(zeroHat,"local",xftof,"FTOF",strlen("local"),
strlen("FTOF"));
165
166 // track vector of this step
167 dx[0] = xin[0] - xout[0];
168 dx[1] = xin[1] - xout[1];
169 dx[2] = xin[2] - xout[2];
170 // length of the track of this step
171 dr = sqrt(dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]);
172 // calculate dEdx only if track length is >0.001 cm
173 if (dr > 1e-3) {
174 dEdx = dEsum/dr;
175 }
176 else {
177 dEdx = 0;
178 }
179
180 /* post the hit to the truth tree */
181 // in other words: store the GENERATED track information
182
183 if ((history == 0) && (plane == 0)) {
184
185 // save all tracks from particles that hit the first plane of FTOF
186 // save the generated "true" values
187
188 int mark = (1<<30) + pointCount;
189
190 // getTwig is defined in src/programs/Simulation/HDGeant/bintree.c
191 // the two arguments are a pointer to a pointer and an integer
192
193 void** twig = getTwig(&forwardTOFTree, mark);
194 if (*twig == 0) {
195 // make_s_ForwardTOF is defined in src/programs/Analysis/hddm/hddm_s.h
196 // and coded in src/programs/Analysis/hddm/hddm_s.c
197 // the same holds for make_s_FtofTruthPoints
198
199 // make_s_ForwardTOF returns pointer to structure s_ForwardTOF generated memory
200 // tof->ftofCoutners and tof-> ftofTruthPoints are initialized already
201
202 s_ForwardTOF_t* tof = *twig = make_s_ForwardTOF();
203 s_FtofTruthPoints_t* points = make_s_FtofTruthPoints(1);
204 tof->ftofTruthPoints = points;
205 points->in[0].primary = (stack == 0);
206 points->in[0].track = track;
207 points->in[0].x = x[0];
208 points->in[0].y = x[1];
209 points->in[0].z = x[2];
210 points->in[0].t = t;
211 points->in[0].px = pin[0]*pin[4];
212 points->in[0].py = pin[1]*pin[4];
213 points->in[0].pz = pin[2]*pin[4];
214 points->in[0].E = pin[3];
215 points->in[0].ptype = ipart;
216 points->mult = 1;
217 pointCount++;
218 }
219 }
220
221 /* post the hit to the hits tree, mark slab as hit */
222 // in other words now store the simulated detector response
223 if (dEsum > 0) {
224 int nhit;
225 s_FtofNorthTruthHits_t* northHits;
226 s_FtofSouthTruthHits_t* southHits;
227
228 s_FtofMCHits_t* noMCHits;
229 s_FtofMCHits_t* soMCHits;
230
231 // getrow and getcolumn are both coded in hddsGeant3.F
232 // see above for function getplane()
233
234 int row = getrow_();
235 int column = getcolumn_();
236
237 // distance of hit from PMT north w.r.t. center and similar for PMT south
238 // this means positive x points north. to get a right handed system y must
239 // point vertically up as z is the beam axis.
240 // plane==0 horizontal plane, plane==1 vertical plane
241 // float dist = xlocal[0];
242
243 float dist = x[1]; // do not use local coordinate for x and y
244 if (plane==1) dist = x[0];
245 float dxnorth = BAR_LENGTH/2.-dist;
246 float dxsouth = BAR_LENGTH/2.+dist;
247
248 // calculate time at the PMT "normalized" to the center, so a hit in the
249 // center will have time "t" at both PMTs
250 // the speed of signal travel is C_EFFECTIVE
251 // propagte time to the end of the bar
252 // column = 0 is a full paddle column ==1,2 is a half paddle
253
254 float tnorth = (column == 2) ? 0 : t + dxnorth/C_EFFECTIVE;
255 float tsouth = (column == 1) ? 0 : t + dxsouth/C_EFFECTIVE;
256
257 // calculate energy seen by PM for this track step using attenuation factor
258 float dEnorth = (column == 2) ? 0 : dEsum * exp(-dxnorth/ATTEN_LENGTH);
259 float dEsouth = (column == 1) ? 0 : dEsum * exp(-dxsouth/ATTEN_LENGTH);
260
261 int mark = (plane<<20) + (row<<10) + column;
262 void** twig = getTwig(&forwardTOFTree, mark);
263
264 if (*twig == 0) { // this paddle has not been hit yet by any particle track
265 // get space and store it
266
267 s_ForwardTOF_t* tof = *twig = make_s_ForwardTOF();
268 s_FtofCounters_t* counters = make_s_FtofCounters(1);
269 counters->mult = 1;
270 counters->in[0].plane = plane;
271 counters->in[0].bar = row;
272 northHits = HDDM_NULL(void*)&hddm_s_nullTarget;
273 southHits = HDDM_NULL(void*)&hddm_s_nullTarget;
274 noMCHits = HDDM_NULL(void*)&hddm_s_nullTarget;
275 soMCHits = HDDM_NULL(void*)&hddm_s_nullTarget;
Value stored to 'soMCHits' is never read
276
277 // get space for the left/top or right/down PMT data for a total
278 // of MAX_HITS possible hits in a single paddle
279 // and space for up to MAX_PAD_HITS hits in a paddle to store MC track information
280 // Note: column=0 means paddle read out on both ends coumn=1or2 means single ended readout
281
282 if (column == 0 || column == 1) {
283 counters->in[0].ftofNorthTruthHits = northHits
284 = make_s_FtofNorthTruthHits(MAX_HITS);
285 }
286 if (column == 0 || column == 2) {
287 counters->in[0].ftofSouthTruthHits = southHits
288 = make_s_FtofSouthTruthHits(MAX_HITS);
289 }
290
291 tof->ftofCounters = counters;
292 counterCount++;
293
294 } else {
295
296 // this paddle is already registered (was hit before)
297 // get the hit list back
298 s_ForwardTOF_t* tof = *twig;
299 northHits = tof->ftofCounters->in[0].ftofNorthTruthHits;
300 southHits = tof->ftofCounters->in[0].ftofSouthTruthHits;
301
302 }
303
304 if (northHits != HDDM_NULL(void*)&hddm_s_nullTarget) {
305
306 // loop over hits in this PM to find correct time slot
307
308 for (nhit = 0; nhit < northHits->mult; nhit++) {
309
310 if (fabs(northHits->in[nhit].t - t) < TWO_HIT_RESOL) {
311 break;
312 }
313 }
314
315 // this hit is within the time frame of a previous hit
316 // combine the times of this weighted by the energy of the hit
317
318 if (nhit < northHits->mult) { /* merge with former hit */
319 northHits->in[nhit].t =
320 (northHits->in[nhit].t * northHits->in[nhit].dE
321 + tnorth * dEnorth)
322 / (northHits->in[nhit].dE += dEnorth);
323
324 // now add MC tracking information
325 // first get MC pointer of this paddle
326
327 noMCHits = northHits->in[nhit].ftofMCHits;
328 unsigned int nMChit = noMCHits->mult;
329 if (nMChit<MAX_PAD_HITS) {
330 noMCHits->in[nMChit].x = x[0];
331 noMCHits->in[nMChit].y = x[1];
332 noMCHits->in[nMChit].z = x[2];
333 noMCHits->in[nMChit].E = pin[3];
334 noMCHits->in[nMChit].px = pin[0]*pin[4];
335 noMCHits->in[nMChit].py = pin[1]*pin[4];
336 noMCHits->in[nMChit].pz = pin[2]*pin[4];
337 noMCHits->in[nMChit].ptype = ipart;
338 noMCHits->in[nMChit].itrack = track;
339 noMCHits->in[nMChit].dist = dist;
340 noMCHits->mult++;
341 }
342
343 } else if (nhit < MAX_HITS){ // hit in new time window
344 northHits->in[nhit].t = tnorth;
345 northHits->in[nhit].dE = dEnorth;
346 northHits->mult++;
347
348 // create memory for MC track hit information
349 northHits->in[nhit].ftofMCHits = noMCHits = make_s_FtofMCHits(MAX_PAD_HITS);
350
351 noMCHits->in[0].x = x[0];
352 noMCHits->in[0].y = x[1];
353 noMCHits->in[0].z = x[2];
354 noMCHits->in[0].E = pin[3];
355 noMCHits->in[0].px = pin[0]*pin[4];
356 noMCHits->in[0].py = pin[1]*pin[4];
357 noMCHits->in[0].pz = pin[2]*pin[4];
358 noMCHits->in[0].ptype = ipart;
359 noMCHits->in[0].itrack = track;
360 noMCHits->in[0].dist = dist;
361 noMCHits->mult = 1;
362
363 } else {
364 fprintf(stderrstderr,"HDGeant error in hitForwardTOF (file hitFTOF.c): ");
365 fprintf(stderrstderr,"max hit count %d exceeded, truncating!\n",MAX_HITS);
366 }
367 }
368
369 if (southHits != HDDM_NULL(void*)&hddm_s_nullTarget) {
370
371 // loop over hits in this PM to find correct time slot
372
373 for (nhit = 0; nhit < southHits->mult; nhit++) {
374
375 if (fabs(southHits->in[nhit].t - t) < TWO_HIT_RESOL) {
376 break;
377 }
378 }
379
380 // this hit is within the time frame of a previous hit
381 // combine the times of this weighted by the energy of the hit
382
383 if (nhit < southHits->mult) { /* merge with former hit */
384 southHits->in[nhit].t =
385 (southHits->in[nhit].t * southHits->in[nhit].dE
386 + tsouth * dEsouth)
387 / (southHits->in[nhit].dE += dEsouth);
388
389 soMCHits = southHits->in[nhit].ftofMCHits;
390
391 // now add MC tracking information
392 unsigned int nMChit = soMCHits->mult;
393 if (nMChit<MAX_PAD_HITS) {
394
395 soMCHits->in[nMChit].x = x[0];
396 soMCHits->in[nMChit].y = x[1];
397 soMCHits->in[nMChit].z = x[2];
398 soMCHits->in[nMChit].E = pin[3];
399 soMCHits->in[nMChit].px = pin[0]*pin[4];
400 soMCHits->in[nMChit].py = pin[1]*pin[4];
401 soMCHits->in[nMChit].pz = pin[2]*pin[4];
402 soMCHits->in[nMChit].ptype = ipart;
403 soMCHits->in[nMChit].itrack = track;
404 soMCHits->in[nMChit].dist = dist;
405 soMCHits->mult++;
406 }
407
408 } else if (nhit < MAX_HITS){ // hit in new time window
409 southHits->in[nhit].t = tsouth;
410 southHits->in[nhit].dE = dEsouth;
411 southHits->mult++;
412
413 // create memory space for MC track hit information
414 southHits->in[nhit].ftofMCHits = soMCHits = make_s_FtofMCHits(MAX_PAD_HITS);
415
416 soMCHits->in[0].x = x[0];
417 soMCHits->in[0].y = x[1];
418 soMCHits->in[0].z = x[2];
419 soMCHits->in[0].E = pin[3];
420 soMCHits->in[0].px = pin[0]*pin[4];
421 soMCHits->in[0].py = pin[1]*pin[4];
422 soMCHits->in[0].pz = pin[2]*pin[4];
423 soMCHits->in[0].ptype = ipart;
424 soMCHits->in[0].itrack = track;
425 soMCHits->in[0].dist = dist;
426 soMCHits->mult = 1;
427
428 } else {
429 fprintf(stderrstderr,"HDGeant error in hitForwardTOF (file hitFTOF.c): ");
430 fprintf(stderrstderr,"max hit count %d exceeded, truncating!\n",MAX_HITS);
431 }
432 }
433 }
434}
435
436/* entry point from fortran */
437
438void hitforwardtof_ (float* xin, float* xout,
439 float* pin, float* pout, float* dEsum,
440 int* track, int* stack, int* history, int* ipart)
441{
442 hitForwardTOF(xin,xout,pin,pout,*dEsum,*track,*stack,*history, *ipart);
443}
444
445
446/* pick and package the hits for shipping */
447// this function is called by loadoutput() (coded in hddmOutput.c)
448// which in turn is called by GUOUT at the end of each event
449
450s_ForwardTOF_t* pickForwardTOF ()
451{
452 s_ForwardTOF_t* box;
453 s_ForwardTOF_t* item;
454
455 if ((counterCount == 0) && (pointCount == 0))
456 {
457 return HDDM_NULL(void*)&hddm_s_nullTarget;
458 }
459
460 box = make_s_ForwardTOF();
461 box->ftofCounters = make_s_FtofCounters(counterCount);
462 box->ftofTruthPoints = make_s_FtofTruthPoints(pointCount);
463
464 while (item = (s_ForwardTOF_t*) pickTwig(&forwardTOFTree)){
465 s_FtofCounters_t* counters = item->ftofCounters;
466 int counter;
467 s_FtofTruthPoints_t* points = item->ftofTruthPoints;
468 int point;
469
470 for (counter=0; counter < counters->mult; ++counter) {
471 s_FtofNorthTruthHits_t* northHits = counters->in[counter].ftofNorthTruthHits;
472 s_FtofSouthTruthHits_t* southHits = counters->in[counter].ftofSouthTruthHits;
473
474 /* compress out the hits below threshold */
475 // cut off parameter is THRESH_MEV
476 int iok,i;
477 int mok=0;
478 // loop over all hits in a counter for the left/up PMT
479 for (iok=i=0; i < northHits->mult; i++) {
480
481 // check threshold
482 if (northHits->in[i].dE >= THRESH_MEV/1e3) {
483
484 if (iok < i) {
485 northHits->in[iok] = northHits->in[i];
486 }
487 ++mok;
488 ++iok;
489 }
490 }
491
492 if (iok) {
493 northHits->mult = iok;
494 } else if (northHits != HDDM_NULL(void*)&hddm_s_nullTarget){ // no hits left over for this PMT
495 counters->in[counter].ftofNorthHits = HDDM_NULL(void*)&hddm_s_nullTarget;
496 FREE(northHits)free(northHits);
497 }
498
499 // now same loop for the right/bottom PMT of a paddle
500 for (iok=i=0; i < southHits->mult; i++){
501
502 if (southHits->in[i].dE >= THRESH_MEV/1e3){
503 if (iok < i){
504 southHits->in[iok] = southHits->in[i];
505 }
506 ++mok;
507 ++iok;
508 }
509 }
510
511 if (iok){
512 southHits->mult = iok;
513 }
514 else if (southHits != HDDM_NULL(void*)&hddm_s_nullTarget) {
515 counters->in[counter].ftofSouthHits = HDDM_NULL(void*)&hddm_s_nullTarget;
516 FREE(southHits)free(southHits);
517 }
518 if (mok){ // total number of time independent FTOF hits in this counter
519 int m = box->ftofCounters->mult++;
520 // add the hit list of this counter to the list
521 box->ftofCounters->in[m] = counters->in[counter];
522 }
523 } // end of loop over all counters
524
525 if (counters != HDDM_NULL(void*)&hddm_s_nullTarget) {
526 FREE(counters)free(counters);
527 }
528
529 // keep also the MC generated primary track particles
530 for (point=0; point < points->mult; ++point) {
531 int m = box->ftofTruthPoints->mult++;
532 box->ftofTruthPoints->in[m] = points->in[point];
533 }
534 if (points != HDDM_NULL(void*)&hddm_s_nullTarget) {
535 FREE(points)free(points);
536 }
537 FREE(item)free(item);
538 }
539
540 // reset the counters
541 counterCount = pointCount = 0;
542
543 // free the hit list memory used by this event
544 if ((box->ftofCounters != HDDM_NULL(void*)&hddm_s_nullTarget) &&
545 (box->ftofCounters->mult == 0)) {
546 FREE(box->ftofCounters)free(box->ftofCounters);
547 box->ftofCounters = HDDM_NULL(void*)&hddm_s_nullTarget;
548 }
549 if ((box->ftofTruthPoints != HDDM_NULL(void*)&hddm_s_nullTarget) &&
550 (box->ftofTruthPoints->mult == 0)) {
551 FREE(box->ftofTruthPoints)free(box->ftofTruthPoints);
552 box->ftofTruthPoints = HDDM_NULL(void*)&hddm_s_nullTarget;
553 }
554 if ((box->ftofCounters->mult == 0) &&
555 (box->ftofTruthPoints->mult == 0)) {
556 FREE(box)free(box);
557 box = HDDM_NULL(void*)&hddm_s_nullTarget;
558 }
559 return box;
560}
561