Bug Summary

File:programs/Simulation/HDGeant/hitFDC.c
Location:line 143, column 29
Description:Memory is never released; potential leak of memory pointed to by 'd'

Annotated Source Code

1/*
2 * hitFDC - registers hits for forward drift chambers
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 hitForwardDC
11 */
12
13#include <stdlib.h>
14#include <stdio.h>
15#include <math.h>
16
17#include <hddm_s.h>
18#include <geant3.h>
19#include <bintree.h>
20
21#include "calibDB.h"
22extern s_HDDM_t* thisInputEvent;
23extern double asic_response(double t);
24extern double Ei(double x);
25
26typedef struct{
27 int writeenohits;
28 int showersincol;
29 int driftclusters;
30}controlparams_t;
31
32extern controlparams_t controlparams_;
33
34
35const float wire_dead_zone_radius[4]={3.0,3.0,3.9,3.9};
36const float strip_dead_zone_radius[4]={1.3,1.3,1.3,1.3};
37
38#define CATHODE_ROT_ANGLE1.309 1.309 // 75 degrees
39
40// Drift speed 2.2cm/us is appropriate for a 90/10 Argon/Methane mixture
41static float DRIFT_SPEED =.0055;
42static float ACTIVE_AREA_OUTER_RADIUS =48.5;
43static float ANODE_CATHODE_SPACING =0.5;
44static float TWO_HIT_RESOL =1.;
45static int WIRES_PER_PLANE =96;
46static float WIRE_SPACING =1.0;
47static float U_OF_WIRE_ZERO =0;//(-((WIRES_PER_PLANE-1.)*WIRE_SPACING)/2)
48static float STRIPS_PER_PLANE =192;
49static float STRIP_SPACING =0.5;
50static float U_OF_STRIP_ZERO =0;// (-((STRIPS_PER_PLANE-1.)*STRIP_SPACING)/2)
51static float STRIP_GAP =0.1;
52static int MAX_HITS =100;
53static float K2 =1.15;
54static float STRIP_NODES = 3;
55static float THRESH_KEV =1. ;
56static float THRESH_ANODE = 1.;
57static float THRESH_STRIPS =5. ; /* pC */
58static float ELECTRON_CHARGE =1.6022e-4; /* fC */
59static float DIFFUSION_COEFF = 1.1e-6; // cm^2/s --> 200 microns at 1 cm
60static float FDC_TIME_WINDOW = 1000.0; //time window for accepting FDC hits, ns
61static float GAS_GAIN = 8e4;
62
63
64#if 0
65static float wire_dx_offset[2304];
66static float wire_dz_offset[2304];
67#endif
68
69binTree_t* forwardDCTree = 0;
70static int stripCount = 0;
71static int wireCount = 0;
72static int pointCount = 0;
73static int initializedx=0;
74
75void gpoiss_(float*,int*,const int*); // avoid solaris compiler warnings
76void rnorml_(float*,int*);
77
78typedef int (*compfn)(const void*, const void*);
79
80// Sort functions for sorting clusters
81int fdc_anode_cluster_sort(const void *a,const void *b){
82 const s_FdcAnodeTruthHit_t *ca=a;
83 const s_FdcAnodeTruthHit_t *cb=b;
84 if (ca->t<cb->t) return -1;
85 else if (ca->t>cb->t) return 1;
86 else return 0;
87}
88int fdc_cathode_cluster_sort(const void *a,const void *b){
89 const s_FdcCathodeTruthHit_t *ca=a;
90 const s_FdcCathodeTruthHit_t *cb=b;
91 if (ca->t<cb->t) return -1;
92 else if (ca->t>cb->t) return 1;
93 else return 0;
94}
95
96
97
98// Locate a position in array xx given x
99void locate(float *xx,int n,float x,int *j){
100 int ju,jm,jl;
101 int ascnd;
102
103 jl=-1;
104 ju=n;
105 ascnd=(xx[n-1]>=xx[0]);
106 while(ju-jl>1){
107 jm=(ju+jl)>>1;
108 if (x>=xx[jm]==ascnd)
109 jl=jm;
110 else
111 ju=jm;
112 }
113 if (x==xx[0]) *j=0;
114 else if (x==xx[n-1]) *j=n-2;
115 else *j=jl;
116}
117
118// Polynomial interpolation on a grid.
119// Adapted from Numerical Recipes in C (2nd Edition), pp. 121-122.
120void polint(float *xa, float *ya,int n,float x, float *y,float *dy){
121 int i,m,ns=0;
122 float den,dif,dift,ho,hp,w;
123
124 float *c=(float *)calloc(n,sizeof(float));
125 float *d=(float *)calloc(n,sizeof(float));
1
Memory is allocated
→
126
127 dif=fabs(x-xa[0]);
128 for (i=0;i<n;i++){
2
←
Loop condition is true. Entering loop body
→
4
←
Loop condition is true. Entering loop body
→
6
←
Loop condition is false. Execution continues on line 136
→
129 if ((dift=fabs(x-xa[i]))<dif){
3
←
Taking false branch
→
5
←
Taking false branch
→
130 ns=i;
131 dif=dift;
132 }
133 c[i]=ya[i];
134 d[i]=ya[i];
135 }
136 *y=ya[ns--];
137
138 for (m=1;m<n;m++){
7
←
Loop condition is true. Entering loop body
→
139 for (i=1;i<=n-m;i++){
8
←
Loop condition is true. Entering loop body
→
140 ho=xa[i-1]-x;
141 hp=xa[i+m-1]-x;
142 w=c[i+1-1]-d[i-1];
143 if ((den=ho-hp)==0.0) return;
9
←
Taking true branch
→
10
←
Memory is never released; potential leak of memory pointed to by 'd'
144
145 den=w/den;
146 d[i-1]=hp*den;
147 c[i-1]=ho*den;
148
149 }
150
151 *y+=(*dy=(2*ns<(n-m) ?c[ns+1]:d[ns--]));
152 }
153 free(c);
154 free(d);
155}
156
157// Simulation of signal on a wire
158double wire_signal(double t,s_FdcAnodeTruthHits_t* ahits){
159 double t0=1.0; // ns; rough order of magnitude
160 int m;
161 double asic_gain=0.76; // mV/fC
162 double func=0;
163 for (m=0;m<ahits->mult;m++){
164 if (t>ahits->in[m].t){
165 double my_time=t-ahits->in[m].t;
166 func+=asic_gain*ahits->in[m].dE*asic_response(my_time);
167 }
168 }
169 return func;
170}
171
172// Simulation of signal on a cathode strip (ASIC output)
173double cathode_signal(double t,s_FdcCathodeTruthHits_t* chits){
174 double t0=1.0; // ns; rough order of magnitude
175 int m;
176 double asic_gain=2.3;
177 double func=0;
178 for (m=0;m<chits->mult;m++){
179 if (t>chits->in[m].t){
180 double my_time=t-chits->in[m].t;
181 func+=asic_gain*chits->in[m].q*asic_response(my_time);
182 }
183 }
184 return func;
185}
186
187// Generate hits in two cathode planes flanking the wire plane
188void AddFDCCathodeHits(int PackNo,float xwire,float avalanche_y,float tdrift,
189 int n_p,int track,int ipart,int chamber,int module,
190 int layer, int global_wire_number){
191
192 s_FdcCathodeTruthHits_t* chits;
193
194 // Anode charge
195 float q_anode;
196 int n_t;
197 // Average number of secondary ion pairs for 40/60 Ar/CO2 mixture
198 float n_s_per_p=1.89;
199 if (controlparams_.driftclusters==0){
200 /* Total number of ion pairs. On average for each primary ion
201 pair produced there are n_s secondary ion pairs produced. The
202 probability distribution is a compound poisson distribution
203 that requires generating two Poisson variables.
204 */
205 int n_s,one=1;
206 float n_s_mean = ((float)n_p)*n_s_per_p;
207 gpoiss_(&n_s_mean,&n_s,&one);
208 n_t = n_s+n_p;
209 q_anode=((float)n_t)*GAS_GAIN*ELECTRON_CHARGE;
210 }
211 else{
212 // Distribute the number of secondary ionizations for this primary
213 // ionization according to a Poisson distribution with mean n_s_over_p.
214 // For simplicity we assume these secondary electrons and the primary
215 // electron stay together as a cluster.
216 int n_s;
217 int one=1;
218 gpoiss_(&n_s_per_p,&n_s,&one);
219 // Anode charge in units of fC
220 n_t=1+n_s;
221 q_anode=GAS_GAIN*ELECTRON_CHARGE*((float)n_t);
222 }
223
224 /* Mock-up of cathode strip charge distribution */
225 int plane, node;
226 for (plane=1; plane<4; plane+=2){
227 float theta = (plane == 1)? -CATHODE_ROT_ANGLE1.309 : +CATHODE_ROT_ANGLE1.309;
228 float cathode_u = xwire*cos(theta)+avalanche_y*sin(theta);
229 int strip1 = ceil((cathode_u-U_OF_STRIP_ZERO)/STRIP_SPACING +0.5);
230 float cathode_u1 = (strip1-1)*STRIP_SPACING + U_OF_STRIP_ZERO;
231 float delta = cathode_u-cathode_u1;
232 float half_gap=ANODE_CATHODE_SPACING;
233
234#if 0
235 half_gap+=(plane==1)?+wire_dz_offset[global_wire_number]:
236 -wire_dz_offset[global_wire_number];
237#endif
238
239 for (node=-STRIP_NODES; node<=STRIP_NODES; node++){
240 /* Induce charge on the strips according to the Mathieson
241 function tuned to results from FDC prototype
242 */
243 float lambda1=(((float)node-0.5)*STRIP_SPACING+STRIP_GAP/2.
244 -delta)/half_gap;
245 float lambda2=(((float)node+0.5)*STRIP_SPACING-STRIP_GAP/2.
246 -delta)/half_gap;
247 float factor=0.25*M_PI3.14159265358979323846*K2;
248 float q = 0.25*q_anode*(tanh(factor*lambda2)-tanh(factor*lambda1));
249
250 int strip = strip1+node;
251 /* Throw away hits on strips falling within a certain dead-zone
252 radius */
253 float strip_outer_u=cathode_u1
254 +(STRIP_SPACING+STRIP_GAP/2.)*(int)node;
255 float cathode_v=-xwire*sin(theta)+avalanche_y*cos(theta);
256 float check_radius=sqrt(strip_outer_u*strip_outer_u
257 +cathode_v*cathode_v);
258
259 if ((strip > 0)
260 && (check_radius>strip_dead_zone_radius[PackNo])
261 && (strip <= STRIPS_PER_PLANE)){
262 int mark = (chamber<<20) + (plane<<10) + strip;
263 void** cathodeTwig = getTwig(&forwardDCTree, mark);
264 if (*cathodeTwig == 0){
265 s_ForwardDC_t* fdc = *cathodeTwig = make_s_ForwardDC();
266 s_FdcChambers_t* chambers = make_s_FdcChambers(1);
267 s_FdcCathodeStrips_t* strips = make_s_FdcCathodeStrips(1);
268 strips->mult = 1;
269 strips->in[0].plane = plane;
270 strips->in[0].strip = strip;
271 strips->in[0].fdcCathodeTruthHits = chits
272 = make_s_FdcCathodeTruthHits(MAX_HITS);
273 chambers->mult = 1;
274 chambers->in[0].module = module;
275 chambers->in[0].layer = layer;
276 chambers->in[0].fdcCathodeStrips = strips;
277 fdc->fdcChambers = chambers;
278 stripCount++;
279 }
280 else{
281 s_ForwardDC_t* fdc = *cathodeTwig;
282 chits = fdc->fdcChambers->in[0].fdcCathodeStrips
283 ->in[0].fdcCathodeTruthHits;
284 }
285
286 int nhit;
287 for (nhit = 0; nhit < chits->mult; nhit++){
288 // To cut down on the number of output clusters, combine
289 // those that would be indistiguishable in time given the
290 // expected timing resolution
291 if (fabs(chits->in[nhit].t - tdrift) <TWO_HIT_RESOL)
292 {
293 break;
294 }
295 }
296 if (nhit < chits->mult) /* merge with former hit */
297 {
298 /* Use the time from the earlier hit but add the charge */
299 chits->in[nhit].q += q;
300 if(chits->in[nhit].t>tdrift){
301 chits->in[nhit].t = tdrift;
302 chits->in[nhit].itrack = track;
303 chits->in[nhit].ptype = ipart;
304 }
305 }
306 else if (nhit < MAX_HITS){ /* create new hit */
307 chits->in[nhit].t = tdrift;
308 chits->in[nhit].q = q;
309 chits->in[nhit].itrack = track;
310 chits->in[nhit].ptype = ipart;
311 chits->mult++;
312 }
313 else{
314 // supress warning
315 /*
316 fprintf(stderr,"HDGeant error in hitForwardDC: ");
317 fprintf(stderr,"max hit count %d exceeded, truncating!\n",
318 MAX_HITS);
319 */
320 }
321
322 }
323 } // loop over cathode strips
324 } // loop over cathode views
325}
326
327
328
329
330
331
332
333
334
335// Add wire information
336int AddFDCAnodeHit(s_FdcAnodeTruthHits_t* ahits,int layer,int ipart,int track,
337 float xwire,float xyz[3],float dE,float t,float *tdrift){
338
339 // Generate 2 random numbers from a Gaussian distribution
340 //
341 float rndno[2];
342 int two=2;
343
344 // Only and always use the built-in Geant random generator,
345 // otherwise debugging is a problem because sequences are not
346 // reproducible from a given pair of random seeds. [rtj]
347
348 /* rnorml_(rndno,&two); */ {
349 float rho,phi1;
350 grndm_(rndno,&two);
351 rho = sqrt(-2*log(rndno[0]));
352 phi1 = rndno[1]*2*M_PI3.14159265358979323846;
353 rndno[0] = rho*cos(phi1);
354 rndno[1] = rho*sin(phi1);
355 }
356
357 // Get the magnetic field at this cluster position
358 float x[3],B[3];
359 transformCoord(xyz,"local",x,"global")transformcoord_(xyz,"local",x,"global",strlen("local"),strlen
("global"));
360 gufld2_(x,B);
361
362 // Find the angle between the wire direction and the direction of the
363 // magnetic field in the x-y plane
364 float wire_dir[2];
365 float wire_theta=1.0472*(float)((layer%3)-1);
366 float phi=0.;;
367 float Br=sqrt(B[0]*B[0]+B[1]*B[1]);
368
369 wire_dir[0]=sin(wire_theta);
370 wire_dir[1]=cos(wire_theta);
371 if (Br>0.) phi= acos((B[0]*wire_dir[0]+B[1]*wire_dir[1])/Br);
372
373 // useful combinations of dx and dz
374 float dx=xyz[0]-xwire;
375 float dx2=dx*dx;
376 float dx4=dx2*dx2;
377 float dz2=xyz[2]*xyz[2];
378 float dz4=dz2*dz2;
379
380 // Next compute the avalanche position along wire.
381 // Correct avalanche position with deflection along wire due to
382 // Lorentz force.
383 xyz[1]+=( 0.1458*B[2]*(1.-0.048*Br) )*dx
384 +( 0.1717+0.01227*B[2] )*(Br*cos(phi))*xyz[2]
385 +( -0.000176 )*dx*dx2/(dz2+0.001);
386 // Add transverse diffusion
387 xyz[1]+=(( 0.01 )*pow(dx2+dz2,0.125)+( 0.0061 )*dx2)*rndno[1];
388
389 // Do not use this cluster if the Lorentz force would deflect
390 // the electrons outside the active region of the detector
391 if (sqrt(xyz[1]*xyz[1]+xwire*xwire)>ACTIVE_AREA_OUTER_RADIUS)
392 return 0;
393
394 // Model the drift time and longitudinal diffusion as a function of
395 // position of the cluster within the cell
396 float tdrift_unsmeared=( 1086.0-106.7*B[2] )*dx2+( 1068.0 )*dz2
397 +dx4*(( -2.675 )/(dz2+0.001)+( 2.4e4 )*dz2);
398 float dt=(( 39.44 )*dx4/(0.5-dz2)+( 56.0 )*dz4/(0.5-dx2)
399 +( 0.01566 )*dx4/(dz4+0.002)/(0.251-dx2))*rndno[1];
400 // Only allow deviations to longer times if near the cell border
401 double dradius=sqrt(dx2+dz2);
402 if (dradius-0.5>=0. && dt<0){
403 dt=0.;
404 }
405
406 // Minimum drift time for docas near wire (very crude approximation)
407 double v_max=0.08; // guess for now based on Garfield, near wire
408 double tmin=dradius/v_max;
409 double tdrift_smeared=tdrift_unsmeared+dt;
410 if (tdrift_smeared<tmin){
411 tdrift_smeared=tmin;
412 }
413
414 // Avalanche time
415 *tdrift=t+tdrift_smeared;
416
417 // Skip cluster if the time would go beyond readout window
418 if (t>FDC_TIME_WINDOW) return 0;
419
420 int nhit;
421
422 // Record the anode hit
423 for (nhit = 0; nhit < ahits->mult; nhit++)
424 {
425 if (fabs(ahits->in[nhit].t - *tdrift) < TWO_HIT_RESOL)
426 {
427 break;
428 }
429 }
430 if (nhit < ahits->mult) /* merge with former hit */
431 {
432 /* use the time from the earlier hit but add the energy */
433 ahits->in[nhit].dE += dE;
434 if(ahits->in[nhit].t>*tdrift){
435 ahits->in[nhit].t = *tdrift;
436 ahits->in[nhit].t_unsmeared=tdrift_unsmeared;
437 ahits->in[nhit].d = sqrt(dx2+dz2);
438
439 ahits->in[nhit].itrack = track;
440 ahits->in[nhit].ptype = ipart;
441 }
442
443 /*ahits->in[nhit].t =
444 (ahits->in[nhit].t * ahits->in[nhit].dE + tdrift * dE)
445 / (ahits->in[nhit].dE += dE);
446 */
447 }
448 else if (nhit < MAX_HITS) /* create new hit */
449 {
450 ahits->in[nhit].t = *tdrift;
451 ahits->in[nhit].t_unsmeared=tdrift_unsmeared;
452 ahits->in[nhit].dE = dE;
453 ahits->in[nhit].d = sqrt(dx2+dz2);
454 ahits->in[nhit].itrack = track;
455 ahits->in[nhit].ptype = ipart;
456 ahits->mult++;
457 }
458 else
459 {
460 fprintf(stderrstderr,"HDGeant error in hitForwardDC: ");
461 fprintf(stderrstderr,"max hit count %d exceeded, truncating!\n",MAX_HITS);
462 }
463
464 return 1;
465}
466
467/* register hits during tracking (from gustep) */
468
469void hitForwardDC (float xin[4], float xout[4],
470 float pin[5], float pout[5], float dEsum,
471 int track, int stack, int history, int ipart)
472{
473 float x[3], t;
474 float dx[3], dr;
475 float dEdx;
476 float xlocal[3];
477 float xinlocal[3];
478 float xoutlocal[3];
479 float dradius;
480 float alpha,sinalpha,cosalpha;
481 int i,j;
482
483 if (!initializedx){
484 mystr_t strings[250];
485 float values[250];
486 int nvalues = 250;
487
488 // Get parameters related to the geometry and the signals
489 int status = GetConstants("FDC/fdc_parms", &nvalues,values,strings);
490 if (!status) {
491 int ncounter = 0;
492 int i;
493 for ( i=0;i<(int)nvalues;i++){
494 //printf("%d %s %f\n",i,strings[i].str,values[i]);
495 if (!strcmp(strings[i].str,"FDC_DRIFT_SPEED")) {
496 DRIFT_SPEED = values[i];
497 ncounter++;
498 }
499 if (!strcmp(strings[i].str,"FDC_ACTIVE_AREA_OUTER_RADIUS")) {
500 ACTIVE_AREA_OUTER_RADIUS = values[i];
501 ncounter++;
502 }
503 if (!strcmp(strings[i].str,"FDC_ANODE_CATHODE_SPACING")) {
504 ANODE_CATHODE_SPACING = values[i];
505 ncounter++;
506 }
507 if (!strcmp(strings[i].str,"FDC_TWO_HIT_RESOL")) {
508 TWO_HIT_RESOL = values[i];
509 ncounter++;
510 }
511 if (!strcmp(strings[i].str,"FDC_WIRES_PER_PLANE")) {
512 WIRES_PER_PLANE = values[i];
513 ncounter++;
514 }
515 if (!strcmp(strings[i].str,"FDC_WIRE_SPACING")) {
516 WIRE_SPACING = values[i];
517 ncounter++;
518 }
519 if (!strcmp(strings[i].str,"FDC_STRIPS_PER_PLANE")) {
520 STRIPS_PER_PLANE = values[i];
521 ncounter++;
522 }
523 if (!strcmp(strings[i].str,"FDC_STRIP_SPACING")) {
524 STRIP_SPACING = values[i];
525 ncounter++;
526 }
527 if (!strcmp(strings[i].str,"FDC_STRIP_GAP")) {
528 STRIP_GAP = values[i];
529 ncounter++;
530 }
531 if (!strcmp(strings[i].str,"FDC_MAX_HITS")) {
532 MAX_HITS = values[i];
533 ncounter++;
534 }
535 if (!strcmp(strings[i].str,"FDC_K2")) {
536 K2 = values[i];
537 ncounter++;
538 }
539 if (!strcmp(strings[i].str,"FDC_STRIP_NODES")) {
540 STRIP_NODES = values[i];
541 ncounter++;
542 }
543 if (!strcmp(strings[i].str,"FDC_THRESH_KEV")) {
544 THRESH_KEV = values[i];
545 ncounter++;
546 }
547 if (!strcmp(strings[i].str,"FDC_THRESH_STRIPS")) {
548 THRESH_STRIPS = values[i];
549 ncounter++;
550 }
551 if (!strcmp(strings[i].str,"FDC_ELECTRON_CHARGE")) {
552 ELECTRON_CHARGE = values[i];
553 ncounter++;
554 }
555 if (!strcmp(strings[i].str,"FDC_DIFFUSION_COEFF")) {
556 DIFFUSION_COEFF = values[i];
557 ncounter++;
558 }
559 }
560 U_OF_WIRE_ZERO = (-((WIRES_PER_PLANE-1.)*WIRE_SPACING)/2);
561 U_OF_STRIP_ZERO = (-((STRIPS_PER_PLANE-1.)*STRIP_SPACING)/2);
562
563 if (ncounter==16){
564 printf("FDC: ALL parameters loaded from Data Base\n");
565 } else if (ncounter<16){
566 printf("FDC: NOT ALL necessary parameters found in Data Base %d out of 16\n",ncounter);
567 } else {
568 printf("FDC: SOME parameters found more than once in Data Base\n");
569 }
570#if 0
571 {
572 int num_values=2304*2;
573 float my_values[2304*2];
574 mystr_t my_strings[2304*2];
575 status=GetArrayConstants("FDC/fdc_wire_offsets",&num_values,
576 my_values,my_strings);
577 if (!status){
578 int i;
579 for (i=0;i<num_values;i+=2){
580 int j=i/2;
581 wire_dx_offset[j]=my_values[i];
582 wire_dz_offset[j]=my_values[i+1];
583 }
584 }
585 }
586#endif
587 }
588 initializedx = 1 ;
589 }
590
591 /* Get chamber information */
592 int layer = getlayer_();
593 if (layer==0){
594 printf("hitFDC: FDC layer number evaluates to zero! THIS SHOULD NEVER HAPPEN! drop this particle.\n");
595 return;
596 }
597 //int module = getmodule_();
598 //int chamber = (module*10)+layer;
599 //int PackNo = (chamber-11)/20;
600 int PackNo = getpackage_()-1;
601 if (PackNo==-1){
602 printf("hitFDC: FDC package number evaluates to zero! THIS SHOULD NEVER HAPPEN! drop this particle.\n");
603 return;
604 }
605 int glayer=6*PackNo+layer-1;
606 int module = 2*(PackNo)+(layer-1)/3+1;
607 int chamber = (module*10)+(layer-1)%3+1;
608 int wire1,wire2;
609 int wire,dwire;
610
611 // transform layer number into Richard's scheme
612 layer=(layer-1)%3+1;
613
614 transformCoord(xin,"global",xinlocal,"local")transformcoord_(xin,"global",xinlocal,"local",strlen("global"
),strlen("local"));
615
616 wire1 = ceil((xinlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
617 transformCoord(xout,"global",xoutlocal,"local")transformcoord_(xout,"global",xoutlocal,"local",strlen("global"
),strlen("local"));
618 wire2 = ceil((xoutlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
619 // Check that wire numbers are not out of range
620 if ((wire1>WIRES_PER_PLANE && wire2==WIRES_PER_PLANE) ||
621 (wire2>WIRES_PER_PLANE && wire1==WIRES_PER_PLANE))
622 wire1=wire2=WIRES_PER_PLANE;
623 if ((wire1==0 && wire2 == 1) || (wire1==1 && wire2== 0)){
624 wire1=wire2=1;
625 }
626 // Make sure at least one wire number is valid
627 if (wire1>WIRES_PER_PLANE&&wire2>WIRES_PER_PLANE) return;
628 if (wire1==0 && wire2==0) return;
629
630 if (wire1>WIRES_PER_PLANE) wire1=wire2;
631 else if (wire2>WIRES_PER_PLANE) wire2=wire1;
632 if (wire1==0) wire1=wire2;
633 else if (wire2==0) wire2=wire1;
634
635 dwire = (wire1 < wire2)? 1 : -1;
636 alpha = atan2(xoutlocal[0]-xinlocal[0],xoutlocal[2]-xinlocal[2]);
637 sinalpha=sin(alpha);
638 cosalpha=cos(alpha);
639 xlocal[0] = (xinlocal[0] + xoutlocal[0])/2;
640 xlocal[1] = (xinlocal[1] + xoutlocal[1])/2;
641 xlocal[2] = (xinlocal[2] + xoutlocal[2])/2;
642
643 wire = ceil((xlocal[0] - U_OF_WIRE_ZERO)/WIRE_SPACING +0.5);
644 x[0] = (xin[0] + xout[0])/2;
645 x[1] = (xin[1] + xout[1])/2;
646 x[2] = (xin[2] + xout[2])/2;
647 t = (xin[3] + xout[3])/2 * 1e9;
648 dx[0] = xin[0] - xout[0];
649 dx[1] = xin[1] - xout[1];
650 dx[2] = xin[2] - xout[2];
651 dr = sqrt(dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2]);
652 if (dr > 1e-3)
653 {
654 dEdx = dEsum/dr;
655 }
656 else
657 {
658 dEdx = 0;
659 }
660
661 /* Make a fuzzy boundary around the forward dead region
662 * by killing any track segment whose midpoint is within the boundary */
663
664 if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
665 < wire_dead_zone_radius[PackNo])
666 {
667 return;
668 }
669
670 /* post the hit to the truth tree */
671
672 if (history == 0)
673 {
674 int mark = (1<<16) + (chamber<<20) + pointCount;
675 void** twig = getTwig(&forwardDCTree, mark);
676 if (*twig == 0)
677 {
678 s_ForwardDC_t* fdc = *twig = make_s_ForwardDC();
679 s_FdcChambers_t* chambers = make_s_FdcChambers(1);
680 s_FdcTruthPoints_t* points = make_s_FdcTruthPoints(1);
681 float xwire = U_OF_WIRE_ZERO + (wire-1)*WIRE_SPACING;
682 float u[2];
683 u[0] = xinlocal[2];
684 u[1] = xinlocal[0]-xwire;
685 dradius = fabs(u[1]*cosalpha-u[0]*sinalpha);
686 points->mult = 1;
687 int a = thisInputEvent->physicsEvents->in[0].reactions->in[0].vertices->in[0].products->mult;
688 points->in[0].primary = (stack <= a);
689 points->in[0].track = track;
690 points->in[0].x = x[0];
691 points->in[0].y = x[1];
692 points->in[0].z = x[2];
693 points->in[0].t = t;
694 points->in[0].px = pin[0]*pin[4];
695 points->in[0].py = pin[1]*pin[4];
696 points->in[0].pz = pin[2]*pin[4];
697 points->in[0].E = pin[3];
698 points->in[0].dradius = dradius;
699 points->in[0].dEdx = dEdx;
700 points->in[0].ptype = ipart;
701 chambers->mult = 1;
702 chambers->in[0].module = module;
703 chambers->in[0].layer = layer;
704 chambers->in[0].fdcTruthPoints = points;
705 fdc->fdcChambers = chambers;
706 pointCount++;
707 }
708 }
709
710 /* post the hit to the hits tree, mark cell as hit */
711
712 if (dEsum > 0)
713 {
714 float sign=1.; // for dealing with the y-position for tracks crossing two cells
715
716 for (wire=wire1; wire-dwire != wire2; wire+=dwire)
717 {
718 int valid_hit=1;
719 float dE,dt;
720 float u[2];
721 float x0[3],x1[3];
722 float avalanche_y;
723 float xwire = U_OF_WIRE_ZERO + (wire-1)*WIRE_SPACING;
724 int global_wire_number=96*glayer+wire-1;
725
726#if 0
727 xwire+=wire_dx_offset[global_wire_number];
728#endif
729
730 if (wire1==wire2){
731 dE=dEsum;
732 x0[0] = xinlocal[0];
733 x0[1] = xinlocal[1];
734 x0[2] = xinlocal[2];
735 x1[0] = xoutlocal[0];
736 x1[1] = xoutlocal[1];
737 x1[2] = xoutlocal[2];
738 }
739 else{
740 x0[0] = xwire-0.5*dwire*WIRE_SPACING;
741 x0[1] = xinlocal[1] + (x0[0]-xinlocal[0]+1e-20)*
742 (xoutlocal[1]-xinlocal[1])/(xoutlocal[0]-xinlocal[0]+1e-20);
743 x0[2] = xinlocal[2] + (x0[0]-xinlocal[0]+1e-20)*
744 (xoutlocal[2]-xinlocal[2])/(xoutlocal[0]-xinlocal[0]+1e-20);
745 if (fabs(x0[2]-xoutlocal[2]) > fabs(xinlocal[2]-xoutlocal[2]))
746 {
747 x0[0] = xinlocal[0];
748 x0[1] = xinlocal[1];
749 x0[2] = xinlocal[2];
750 }
751 x1[0] = xwire+0.5*dwire*WIRE_SPACING;
752 x1[1] = xinlocal[1] + (x1[0]-xinlocal[0]+1e-20)*
753 (xoutlocal[1]-xinlocal[1])/(xoutlocal[0]-xinlocal[0]+1e-20);
754 x1[2] = xinlocal[2] + (x1[0]-xinlocal[0]+1e-20)*
755 (xoutlocal[2]-xinlocal[2])/(xoutlocal[0]-xinlocal[0]+1e-20);
756 if (fabs(x1[2]-xinlocal[2]) > fabs(xoutlocal[2]-xinlocal[2]))
757 {
758 x1[0] = xoutlocal[0];
759 x1[1] = xoutlocal[1];
760 x1[2] = xoutlocal[2];
761 }
762 dE = dEsum*(x1[2]-x0[2])/(xoutlocal[2]-xinlocal[2]);
763 }
764
765 if (dE > 0){
766 s_FdcAnodeTruthHits_t* ahits;
767
768 // Create (or grab) an entry in the tree for the anode wire
769 int mark = (chamber<<20) + (2<<10) + wire;
770 void** twig = getTwig(&forwardDCTree, mark);
771
772 if (*twig == 0)
773 {
774 s_ForwardDC_t* fdc = *twig = make_s_ForwardDC();
775 s_FdcChambers_t* chambers = make_s_FdcChambers(1);
776 s_FdcAnodeWires_t* wires = make_s_FdcAnodeWires(1);
777 wires->mult = 1;
778 wires->in[0].wire = wire;
779 wires->in[0].fdcAnodeTruthHits = ahits = make_s_FdcAnodeTruthHits(MAX_HITS);
780 chambers->mult = 1;
781 chambers->in[0].module = module;
782 chambers->in[0].layer = layer;
783 chambers->in[0].fdcAnodeWires = wires;
784 fdc->fdcChambers = chambers;
785 wireCount++;
786 }
787 else
788 {
789 s_ForwardDC_t* fdc = *twig;
790 ahits = fdc->fdcChambers->in[0].fdcAnodeWires->in[0].fdcAnodeTruthHits;
791 }
792
793
794 float rndno[2];
795 int two=2;
796
797 // Find the number of primary ion pairs:
798 /* The total number of ion pairs depends on the energy deposition
799 and the effective average energy to produce a pair, w_eff.
800 On average for each primary ion pair produced there are n_s_per_p
801 secondary ion pairs produced.
802 */
803 int one=1;
804 // Average number of secondary ion pairs for 40/60 Ar/CO2 mixture
805 float n_s_per_p=1.89;
806 //Average energy needed to produce an ion pair for 50/50 mixture
807 float w_eff=30.2e-9; // GeV
808 // Average number of primary ion pairs
809 float n_p_mean = dE/w_eff/(1.+n_s_per_p);
810 int n_p; // number of primary ion pairs
811 gpoiss_(&n_p_mean,&n_p,&one);
812
813 // Drift time
814 float tdrift=0;
815
816 if (controlparams_.driftclusters==0){
817 float zrange=x1[2]-x0[2];
818 float tany=(x1[1]-x0[1])/zrange;
819 float tanx=(x1[0]-x0[0])/zrange;
820 float dz=ANODE_CATHODE_SPACING-dradius*sign*sinalpha;
821#if 0
822 dz+=wire_dz_offset[global_wire_number];
823#endif
824 xlocal[0]=x0[0]+tanx*dz;
825 if (fabs(xlocal[0]-xwire)>0.5){
826 xlocal[0]=x1[0];
827 xlocal[1]=x1[1];
828 xlocal[2]=x1[2];
829 }
830 else{
831 xlocal[1]=x0[1]+tany*dz;
832 xlocal[2]=x0[2]+dz;
833 }
834
835 /* If the cluster position is within the wire-deadened region of the
836 detector, skip this cluster
837 */
838 if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
839 >=wire_dead_zone_radius[PackNo]){
840 if (AddFDCAnodeHit(ahits,layer,ipart,track,xwire,xlocal,dE,t,
841 &tdrift)){
842 AddFDCCathodeHits(PackNo,xwire,xlocal[1],tdrift,n_p,track,ipart,
843 chamber,module,layer,global_wire_number);
844 }
845
846 }
847 }
848 else{
849 // Loop over the number of primary ion pairs
850 int n;
851 for (n=0;n<n_p;n++){
852 // Generate a cluster at a random position along the path with cell
853 float rndno[2];
854 grndm_(rndno,&two);
855 float dzrand=(x1[2]-x0[2])*rndno[0];
856 // Position of the cluster
857 xlocal[0]=x0[0]+(x1[0]-x0[0])*rndno[0];
858 xlocal[1]=x0[1]+(x1[1]-x0[1])*rndno[0];
859 xlocal[2]=x0[2]+dzrand;
860 /* If the cluster position is within the wire-deadened region of the
861 detector, skip this cluster
862 */
863 if (sqrt(xlocal[0]*xlocal[0]+xlocal[1]*xlocal[1])
864 >=wire_dead_zone_radius[PackNo]){
865 if (AddFDCAnodeHit(ahits,layer,ipart,track,xwire,xlocal,dE,t,
866 &tdrift)){
867 AddFDCCathodeHits(PackNo,xwire,xlocal[1],tdrift,n_p,track,ipart,
868 chamber,module,layer,global_wire_number);
869 }
870 }
871
872 } // loop over primary ion pairs
873 }
874 } // Check for non-zero energy
875
876 sign*=-1; // for dealing with the y-position for tracks crossing two cells
877 } // loop over wires
878 } // Check that total energy deposition is not zero
879}
880
881/* entry points from fortran */
882
883void hitforwarddc_(float* xin, float* xout,
884 float* pin, float* pout, float* dEsum,
885 int* track, int* stack, int* history, int* ipart)
886{
887 hitForwardDC(xin,xout,pin,pout,*dEsum,*track,*stack,*history,*ipart);
888}
889
890
891/* pick and package the hits for shipping */
892
893s_ForwardDC_t* pickForwardDC ()
894{
895 s_ForwardDC_t* box;
896 s_ForwardDC_t* item;
897
898 if ((stripCount == 0) && (wireCount == 0) && (pointCount == 0))
899 {
900 return HDDM_NULL(void*)&hddm_s_nullTarget;
901 }
902
903 box = make_s_ForwardDC();
904 box->fdcChambers = make_s_FdcChambers(32);
905 box->fdcChambers->mult = 0;
906 while (item = (s_ForwardDC_t*) pickTwig(&forwardDCTree))
907 {
908 s_FdcChambers_t* chambers = item->fdcChambers;
909 int module = chambers->in[0].module;
910 int layer = chambers->in[0].layer;
911 int m = box->fdcChambers->mult;
912
913 /* compress out the hits below threshold */
914 s_FdcAnodeWires_t* wires = chambers->in[0].fdcAnodeWires;
915 int wire;
916 s_FdcCathodeStrips_t* strips = chambers->in[0].fdcCathodeStrips;
917 int strip;
918 s_FdcTruthPoints_t* points = chambers->in[0].fdcTruthPoints;
919 int point;
920 int mok=0;
921 for (wire=0; wire < wires->mult; wire++)
922 {
923 s_FdcAnodeTruthHits_t* ahits = wires->in[wire].fdcAnodeTruthHits;
924
925 // Sort the clusters by time
926 qsort(ahits->in,ahits->mult,sizeof(s_FdcAnodeTruthHit_t),
927 (compfn)fdc_anode_cluster_sort);
928
929 int i,iok=0;
930
931 if (controlparams_.driftclusters==0){
932 for (iok=i=0; i < ahits->mult; i++)
933 {
934 if (ahits->in[i].dE >= THRESH_KEV/1e6)
935 {
936 if (iok < i)
937 {
938 ahits->in[iok] = ahits->in[i];
939 }
940 ++iok;
941 ++mok;
942 }
943 }
944 }
945 else{ // Simulate clusters within the cell
946
947 // printf("-------------\n");
948
949
950 // Temporary histogram in 1 ns bins to store waveform data
951 int num_samples=(int)FDC_TIME_WINDOW;
952 float *samples=(float *)malloc(num_samples*sizeof(float));
953 for (i=0;i<num_samples;i++){
954 samples[i]=wire_signal((float)i,ahits);
955 //printf("%f %f\n",(float)i,samples[i]);
956 }
957
958 int returned_to_baseline=0;
959 float q=0;
960 for (i=0;i<num_samples;i++){
961 if (samples[i]>=THRESH_ANODE){
962 if (returned_to_baseline==0){
963 ahits->in[iok].itrack = ahits->in[0].itrack;
964 ahits->in[iok].ptype = ahits->in[0].ptype;
965
966 // Do an interpolation to find the time at which the threshold
967 // was crossed.
968 float t_array[4];
969 int k;
970 float my_t,my_terr;
971 for (k=0;k<4;k++) t_array[k]=i-1+k;
972 polint(&samples[i-1],t_array,4,THRESH_ANODE,&my_t,&my_terr);
973 ahits->in[iok].t=my_t;
974
975 returned_to_baseline=1;
976 iok++;
977 mok++;
978 }
979 q+=samples[i];
980 }
981 if (returned_to_baseline
982 && (samples[i]<THRESH_ANODE)){
983 returned_to_baseline=0;
984 }
985 }
986 free(samples);
987 } // Simulation of clusters within cell
988
989 if (iok)
990 {
991 ahits->mult = iok;
992 }
993 else if (ahits != HDDM_NULL(void*)&hddm_s_nullTarget)
994 {
995 FREE(ahits)free(ahits);
996 }
997 }
998 if ((wires != HDDM_NULL(void*)&hddm_s_nullTarget) && (mok == 0))
999 {
1000 FREE(wires)free(wires);
1001 wires = HDDM_NULL(void*)&hddm_s_nullTarget;
1002 }
1003
1004
1005 mok = 0;
1006 for (strip=0; strip < strips->mult; strip++)
1007 {
1008 s_FdcCathodeTruthHits_t* chits = strips->in[strip].fdcCathodeTruthHits;
1009
1010 // Sort the clusters by time
1011 qsort(chits->in,chits->mult,sizeof(s_FdcCathodeTruthHit_t),
1012 (compfn)fdc_cathode_cluster_sort);
1013
1014 int i,iok=0;
1015
1016 if (controlparams_.driftclusters==0){
1017 for (iok=i=0; i < chits->mult; i++)
1018 {
1019 if (chits->in[i].q >= THRESH_STRIPS)
1020 {
1021 if (iok < i)
1022 {
1023 chits->in[iok] = chits->in[i];
1024 }
1025 ++iok;
1026 ++mok;
1027 }
1028 }
1029
1030 }
1031 else{
1032
1033 // Temporary histogram in 1 ns bins to store waveform data
1034 int num_samples=(int)(FDC_TIME_WINDOW);
1035 float *samples=(float *)malloc(num_samples*sizeof(float));
1036 for (i=0;i<num_samples;i++){
1037 samples[i]=cathode_signal((float)i,chits);
1038 //printf("t %f V %f\n",(float)i,samples[i]);
1039 }
1040
1041 int threshold_toggle=0;
1042 int istart=0;
1043 float q=0;
1044 int FADC_BIN_SIZE=1;
1045 for (i=0;i<num_samples;i+=FADC_BIN_SIZE){
1046 if (samples[i]>=THRESH_STRIPS){
1047 if (threshold_toggle==0){
1048 chits->in[iok].itrack = chits->in[0].itrack;
1049 chits->in[iok].ptype = chits->in[0].ptype;
1050 chits->in[iok].t=(float) i;
1051 //chits->in[iok].q=samples[i];
1052 istart=i-1;
1053 threshold_toggle=1;
1054 //iok++;
1055 //mok++;
1056 }
1057 }
1058 if (threshold_toggle &&
1059 (samples[i]<THRESH_STRIPS)){
1060 int j;
1061 // Find the first peak
1062 for (j=istart+1;j<i-1;j++){
1063 if (samples[j]>samples[j-1] && samples[j]>samples[j+1]){
1064 chits->in[iok].q=samples[j];
1065 break;
1066 }
1067 }
1068 threshold_toggle=0;
1069 iok++;
1070 mok++;
1071 //break;
1072 }
1073 }
1074 i=num_samples-1;
1075 if (samples[i]>=THRESH_STRIPS&&threshold_toggle){
1076 int j;
1077 for (j=istart+1;j<i-1;j++){
1078 if (samples[j]>samples[j-1] && samples[j]>samples[j+1]){
1079 chits->in[iok].q=samples[j];
1080 break;
1081 }
1082 }
1083 }
1084
1085 free(samples);
1086 }// Simulate clusters within cell
1087
1088 if (iok)
1089 {
1090 chits->mult = iok;
1091 //chits->mult=1;
1092 }
1093 else if (chits != HDDM_NULL(void*)&hddm_s_nullTarget)
1094 {
1095 FREE(chits)free(chits);
1096 }
1097
1098 }
1099 if ((strips != HDDM_NULL(void*)&hddm_s_nullTarget) && (mok == 0))
1100 {
1101 FREE(strips)free(strips);
1102 strips = HDDM_NULL(void*)&hddm_s_nullTarget;
1103 }
1104
1105 if ((wires != HDDM_NULL(void*)&hddm_s_nullTarget) ||
1106 (strips != HDDM_NULL(void*)&hddm_s_nullTarget) ||
1107 (points != HDDM_NULL(void*)&hddm_s_nullTarget))
1108 {
1109 if ((m == 0) || (module > box->fdcChambers->in[m-1].module)
1110 || (layer > box->fdcChambers->in[m-1].layer))
1111 {
1112 box->fdcChambers->in[m] = chambers->in[0];
1113
1114 box->fdcChambers->in[m].fdcCathodeStrips =
1115 make_s_FdcCathodeStrips(stripCount);
1116 box->fdcChambers->in[m].fdcAnodeWires =
1117 make_s_FdcAnodeWires(wireCount);
1118 box->fdcChambers->in[m].fdcTruthPoints =
1119 make_s_FdcTruthPoints(pointCount);
1120 box->fdcChambers->mult++;
1121 }
1122 else
1123 {
1124 m--;
1125 }
1126 for (strip=0; strip < strips->mult; ++strip)
1127 {
1128 int mm = box->fdcChambers->in[m].fdcCathodeStrips->mult++;
1129 box->fdcChambers->in[m].fdcCathodeStrips->in[mm] = strips->in[strip];
1130 }
1131 if (strips != HDDM_NULL(void*)&hddm_s_nullTarget)
1132 {
1133 FREE(strips)free(strips);
1134 }
1135 for (wire=0; wire < wires->mult; ++wire)
1136 {
1137 int mm = box->fdcChambers->in[m].fdcAnodeWires->mult++;
1138 box->fdcChambers->in[m].fdcAnodeWires->in[mm] = wires->in[wire];
1139 }
1140 if (wires != HDDM_NULL(void*)&hddm_s_nullTarget)
1141 {
1142 FREE(wires)free(wires);
1143 }
1144 for (point=0; point < points->mult; ++point)
1145 {
1146 int mm = box->fdcChambers->in[m].fdcTruthPoints->mult++;
1147 box->fdcChambers->in[m].fdcTruthPoints->in[mm] = points->in[point];
1148 }
1149 if (points != HDDM_NULL(void*)&hddm_s_nullTarget)
1150 {
1151 FREE(points)free(points);
1152 }
1153 }
1154 FREE(chambers)free(chambers);
1155 FREE(item)free(item);
1156 }
1157
1158 stripCount = wireCount = pointCount = 0;
1159
1160 if ((box->fdcChambers != HDDM_NULL(void*)&hddm_s_nullTarget) &&
1161 (box->fdcChambers->mult == 0))
1162 {
1163 FREE(box->fdcChambers)free(box->fdcChambers);
1164 box->fdcChambers = HDDM_NULL(void*)&hddm_s_nullTarget;
1165 }
1166 if (box->fdcChambers->mult == 0)
1167 {
1168 FREE(box)free(box);
1169 box = HDDM_NULL(void*)&hddm_s_nullTarget;
1170 }
1171 return box;
1172}