Bug Summary

File:programs/Simulation/mcsmear/smear.cc
Location:line 406, column 3
Description:Argument to free() is the address of the global variable 'hddm_s_nullTarget', which is not memory allocated by malloc()

Annotated Source Code

1// $Id: smear.cc 8693 2012-01-09 16:04:25Z davidl $
2//
3// Created June 22, 2005 David Lawrence
4
5#include <iostream>
6#include <iomanip>
7#include <vector>
8using namespace std;
9
10#include <FCAL/DFCALGeometry.h>
11#include <CCAL/DCCALGeometry.h>
12#include <BCAL/DBCALGeometry.h>
13
14#include <math.h>
15#include "units.h"
16#include "HDDM/hddm_s.h"
17#include <TF1.h>
18#include <TH2.h>
19
20#include "DRandom2.h"
21
22#ifndef _DBG_std::cerr<<"smear.cc"<<":"<<22<<" "
23#define _DBG_std::cerr<<"smear.cc"<<":"<<23<<" " cout<<__FILE__"smear.cc"<<":"<<__LINE__23<<" "
24#define _DBG__std::cerr<<"smear.cc"<<":"<<24<<std::
endl cout<<__FILE__"smear.cc"<<":"<<__LINE__24<<endl
25#endif
26
27
28extern vector<vector<float> >fdc_smear_parms;
29extern TH2F *fdc_drift_time_smear_hist;
30extern TH2F *fdc_drift_dist_smear_hist;
31extern TH2F *fdc_drift_time;
32extern TH2F *cdc_drift_time;
33extern TH1F *fdc_anode_mult;
34extern TH2F *cdc_drift_smear;
35
36void GetAndSetSeeds(s_HDDM_t *hddm_s);
37void SmearCDC(s_HDDM_t *hddm_s);
38void AddNoiseHitsCDC(s_HDDM_t *hddm_s);
39void SmearFDC(s_HDDM_t *hddm_s);
40void AddNoiseHitsFDC(s_HDDM_t *hddm_s);
41void SmearFCAL(s_HDDM_t *hddm_s);
42void SmearCCAL(s_HDDM_t *hddm_s);
43void SmearBCAL(s_HDDM_t *hddm_s);
44void SmearTOF(s_HDDM_t *hddm_s);
45void SmearSTC(s_HDDM_t *hddm_s);
46void SmearCherenkov(s_HDDM_t *hddm_s);
47void InitCDCGeometry(void);
48void InitFDCGeometry(void);
49
50pthread_mutex_t mutex_fdc_smear_function = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } };
51
52bool CDC_GEOMETRY_INITIALIZED = false;
53int CDC_MAX_RINGS=0;
54vector<unsigned int> NCDC_STRAWS;
55vector<double> CDC_RING_RADIUS;
56
57DFCALGeometry *fcalGeom = NULL__null;
58DCCALGeometry *ccalGeom = NULL__null;
59bool FDC_GEOMETRY_INITIALIZED = false;
60unsigned int NFDC_WIRES_PER_PLANE;
61vector<double> FDC_LAYER_Z;
62double FDC_RATE_COEFFICIENT;
63
64double SampleGaussian(double sigma);
65double SamplePoisson(double lambda);
66double SampleRange(double x1, double x2);
67
68// Do we or do we not add noise hits
69extern bool ADD_NOISE;
70
71// Do we or do we not smear real hits
72extern bool SMEAR_HITS;
73
74// Use or ignore random number seeds found in HDDM file
75extern bool IGNORE_SEEDS;
76
77// The error on the drift time in the CDC. The drift times
78// for the actual CDC hits coming from the input file
79// are smeared by a gaussian with this sigma.
80extern double CDC_TDRIFT_SIGMA;
81
82// The time window for which CDC hits are accumulated.
83// This is used to determine the number of background
84// hits in the CDC for a given event.
85extern double CDC_TIME_WINDOW;
86
87// The error in the energy deposition measurement in the CDC due to pedestal noise
88extern double CDC_PEDESTAL_SIGMA;
89
90// If the following flag is true, then include the drift-distance
91// dependency on the error in the FDC position. Otherwise, use a
92// flat distribution given by the FDC_TDRIFT_SIGMA below.
93extern bool FDC_USE_PARAMETERIZED_SIGMA;
94
95// The error on the drift time in the FDC. The drift times
96// for the actual FDC hits coming from the input file
97// are smeared by a gaussian with this sigma.
98extern double FDC_TDRIFT_SIGMA;
99
100// The error in the distance along the wire as measured by
101// the cathodes. This should NOT include the Lorentz
102// effect which is already included in hdgeant. It
103// should include any fluctuations due to ion trail density
104// etc.
105extern double FDC_CATHODE_SIGMA;
106
107// The FDC pedestal noise is used to smear the cathode ADC
108// values such that the position along the wire has the resolution
109// specified by FDC_CATHODE_SIGMA.
110extern double FDC_PED_NOISE; //pC (calculated from FDC_CATHODE_SIGMA in SmearFDC)
111
112// If energy loss was turned off in the FDC then the pedestal
113// noise will not be scaled properly to give the nominal 200 micron
114// resolution along the wires. This flag is used to indicated
115// the magic scale factor should be applied to FDC_PED_NOISE
116// when it is calculated below to get the correct resolution.
117extern bool FDC_ELOSS_OFF;
118
119// Time window for acceptance of FDC hits
120extern double FDC_TIME_WINDOW;
121
122// Fraction of FDC hits to randomly drop (0=drop nothing 1=drop everything)
123extern double FDC_HIT_DROP_FRACTION;
124
125// Photon-statistics factor for smearing hit energy (from Criss's MC)
126// (copied from DFCALMCResponse_factory.cc 7/2/2009 DL)
127extern double FCAL_PHOT_STAT_COEF;
128
129// Single block energy threshold (applied after smearing)
130extern double FCAL_BLOCK_THRESHOLD;
131
132// Photon-statistics factor for smearing hit energy for CompCal
133// (This is just a rough estimate 11/30/2010 DL)
134double CCAL_PHOT_STAT_COEF = 0.035/2.0;
135
136// Single block energy threshold (applied after smearing)
137// (This is just a rough estimate 11/30/2010 DL)
138double CCAL_BLOCK_THRESHOLD = 20.0*k_MeV;
139
140
141// Forward TOF resolution
142extern double TOF_SIGMA;
143extern double TOF_PHOTONS_PERMEV;
144
145// Start counter resolution
146extern double START_SIGMA ;
147extern double START_PHOTONS_PERMEV;
148
149extern bool SMEAR_BCAL;
150
151// Polynomial interpolation on a grid.
152// Adapted from Numerical Recipes in C (2nd Edition), pp. 121-122.
153void polint(float *xa, float *ya,int n,float x, float *y,float *dy){
154 int i,m,ns=0;
155 float den,dif,dift,ho,hp,w;
156
157 float *c=(float *)calloc(n,sizeof(float));
158 float *d=(float *)calloc(n,sizeof(float));
159
160 dif=fabs(x-xa[0]);
161 for (i=0;i<n;i++){
162 if ((dift=fabs(x-xa[i]))<dif){
163 ns=i;
164 dif=dift;
165 }
166 c[i]=ya[i];
167 d[i]=ya[i];
168 }
169 *y=ya[ns--];
170
171 for (m=1;m<n;m++){
172 for (i=1;i<=n-m;i++){
173 ho=xa[i-1]-x;
174 hp=xa[i+m-1]-x;
175 w=c[i+1-1]-d[i-1];
176 if ((den=ho-hp)==0.0) return;
177
178 den=w/den;
179 d[i-1]=hp*den;
180 c[i-1]=ho*den;
181
182 }
183
184 *y+=(*dy=(2*ns<(n-m) ?c[ns+1]:d[ns--]));
185 }
186 free(c);
187 free(d);
188}
189
190//-----------
191// Smear
192//-----------
193void Smear(s_HDDM_t *hddm_s)
194{
195 GetAndSetSeeds(hddm_s);
196
197 if(SMEAR_HITS)SmearCDC(hddm_s);
198 if(ADD_NOISE)AddNoiseHitsCDC(hddm_s);
199 if(SMEAR_HITS)SmearFDC(hddm_s);
200 if(ADD_NOISE)AddNoiseHitsFDC(hddm_s);
201 if(SMEAR_HITS)SmearFCAL(hddm_s);
202 if(SMEAR_HITS)SmearCCAL(hddm_s);
203 if(SMEAR_BCAL)SmearBCAL(hddm_s);
204 if(SMEAR_HITS)SmearTOF(hddm_s);
205 if(SMEAR_HITS)SmearSTC(hddm_s);
206 if(SMEAR_HITS)SmearCherenkov(hddm_s);
207}
208
209//-----------
210// GetAndSetSeeds
211//-----------
212void GetAndSetSeeds(s_HDDM_t *hddm_s)
213{
214 // Check if a non-zero seed for mcsmear exists in the input
215 // HDDM file. If so, use it to set the seeds for the random number
216 // generator. Then, make sure the seeds that are used are stored
217 // in the output event.
218
219 if(hddm_s == NULL__null)return;
220 if(hddm_s->physicsEvents == NULL__null || hddm_s->physicsEvents==HDDM_NULL(void*)&hddm_s_nullTarget)return;
221 if(hddm_s->physicsEvents->mult<1)return;
222 s_PhysicsEvent_t *pe = &hddm_s->physicsEvents->in[0];
223 if(pe->reactions == NULL__null || pe->reactions==HDDM_NULL(void*)&hddm_s_nullTarget)return;
224 if(pe->reactions->mult<1)return;
225 s_Random_t *my_rand = pe->reactions->in[0].random;
226 UInt_t seed1, seed2, seed3;
227 if(my_rand == NULL__null || my_rand==HDDM_NULL(void*)&hddm_s_nullTarget){
228 // No seeds stored in event. Add them
229 my_rand = pe->reactions->in[0].random = make_s_Random();
230
231 }else{
232 if(!IGNORE_SEEDS){
233 // Copy seeds from file to local variables
234 seed1 = *((UInt_t*)&my_rand->seed_mcsmear1);
235 seed2 = *((UInt_t*)&my_rand->seed_mcsmear2);
236 seed3 = *((UInt_t*)&my_rand->seed_mcsmear3);
237
238 // Set the seeds in the random generator with those found
239 // in the file, but ONLY if they are not all zeros
240 if((seed1+seed2+seed3) != 0)gDRandom.SetSeeds(seed1, seed2, seed3);
241 }
242 }
243
244 // Copy seeds from generator to local variables
245 gDRandom.GetSeeds(seed1, seed2, seed3);
246
247 // Copy seeds from local variables to HDDM file
248 my_rand->seed_mcsmear1 = *((Int_t*)&seed1);
249 my_rand->seed_mcsmear2 = *((Int_t*)&seed2);
250 my_rand->seed_mcsmear3 = *((Int_t*)&seed3);
251}
252
253//-----------
254// SmearCDC
255//-----------
256void SmearCDC(s_HDDM_t *hddm_s)
257{
258 /// Smear the drift times of all CDC hits.
259 /// This will add cdcStrawHit objects generated by smearing values in the
260 /// cdcStrawTruthHit objects that hdgeant outputs. Any existing cdcStrawHit
261 /// objects will be replaced.
262
263 // Acquire the pointer to the physics events
264 s_PhysicsEvents_t* allEvents = hddm_s->physicsEvents;
265 if(!allEvents)return;
266
267 for (unsigned int m=0; m < allEvents->mult; m++) {
268 // Acquire the pointer to the overall hits section of the data
269 s_HitView_t *hits = allEvents->in[m].hitView;
270
271 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget)return;
272 if (hits->centralDC == HDDM_NULL(void*)&hddm_s_nullTarget)return;
273 if (hits->centralDC->cdcStraws == HDDM_NULL(void*)&hddm_s_nullTarget)return;
274 for(unsigned int k=0; k<hits->centralDC->cdcStraws->mult; k++){
275 s_CdcStraw_t *cdcstraw = &hits->centralDC->cdcStraws->in[k];
276
277 // If the event already contains a s_CdcStrawHits_t object then free it.
278 if(cdcstraw->cdcStrawHits!=HDDM_NULL(void*)&hddm_s_nullTarget){
279 static bool warned = false;
280 FREE(cdcstraw->cdcStrawHits)free(cdcstraw->cdcStrawHits);
281 if(!warned){
282 warned = true;
283 cerr<<endl;
284 cerr<<"WARNING: CDC hits already exist in input file! Overwriting!"<<endl;
285 cerr<<endl;
286 }
287 }
288
289 // Create new s_CdcStrawHits_t
290 cdcstraw->cdcStrawHits = make_s_CdcStrawHits(cdcstraw->cdcStrawTruthHits->mult);
291 cdcstraw->cdcStrawHits->mult = cdcstraw->cdcStrawTruthHits->mult;
292
293 for(unsigned int j=0; j<cdcstraw->cdcStrawTruthHits->mult; j++){
294 s_CdcStrawHit_t *strawhit = &cdcstraw->cdcStrawHits->in[j];
295 s_CdcStrawTruthHit_t *strawtruthhit = &cdcstraw->cdcStrawTruthHits->in[j];
296
297 double sigma_t = CDC_TDRIFT_SIGMA;
298 // Smear out the CDC drift time using the specified sigma.
299 // This is for timing resolution from the electronics; diffusion is handled in hdgeant.
300 double delta_t = SampleGaussian(sigma_t)*1.0E9; // delta_t is in ns
301 strawhit->t = strawtruthhit->t + delta_t;
302
303 if (j==0){
304 cdc_drift_time->Fill(strawhit->t,strawtruthhit->d);
305
306 cdc_drift_smear->Fill(strawtruthhit->d,
307 strawtruthhit->d-(0.02887*sqrt(strawhit->t)-1.315e-5*strawtruthhit->t));
308 }
309
310 // Smear energy deposition
311 strawhit->dE = strawtruthhit->dE+SampleGaussian(CDC_PEDESTAL_SIGMA);
312
313 // to be consistent initialize the value d=DOCA to zero
314 strawhit->d = 0.;
315 strawhit->itrack = strawtruthhit->itrack;
316 strawhit->ptype = strawtruthhit->ptype;
317 // If the time is negative, reject this smear and try again
318 //if(strawhit->t<0)j--;
319 }
320 }
321 }
322}
323
324//-----------
325// AddNoiseHitsCDC
326//-----------
327void AddNoiseHitsCDC(s_HDDM_t *hddm_s)
328{
329 if(!CDC_GEOMETRY_INITIALIZED)InitCDCGeometry();
1
Taking false branch
→
330
331 // Calculate the number of noise hits for each straw and store
332 // them in a sparse map. We must do it this way since we have to know
333 // the total number of CdcStraw_t structures to allocate in our
334 // call to make_s_CdcStraws.
335 //
336 // The straw rates are obtained using a parameterization done
337 // to calculate the event size for the August 29, 2007 online
338 // meeting. This parameterization is almost already obsolete.
339 // 10/12/2007 D. L.
340 vector<int> Nstraw_hits;
341 vector<int> straw_number;
342 vector<int> ring_number;
343 int Nnoise_straws = 0;
344 int Nnoise_hits = 0;
345 for(unsigned int ring=1; ring<=NCDC_STRAWS.size(); ring++){
2
←
Loop condition is false. Execution continues on line 365
→
346 double p[2] = {10.4705, -0.103046};
347 double r_prime = (double)(ring+3);
348 double N = exp(p[0] + r_prime*p[1]);
349 N *= CDC_TIME_WINDOW;
350 for(unsigned int straw=1; straw<=NCDC_STRAWS[ring-1]; straw++){
351 // Indivdual straw rates should be way less than 1/event so
352 // we just use the rate as a probablity.
353 double Nhits = SampleRange(0.0, 1.0)<N ? 1.0:0.0;
354 if(Nhits<1.0)continue;
355 int iNhits = (int)floor(Nhits);
356 Nstraw_hits.push_back(iNhits);
357 straw_number.push_back(straw);
358 ring_number.push_back(ring);
359 Nnoise_straws++;
360 Nnoise_hits+=iNhits;
361 }
362 }
363
364 // Loop over Physics Events
365 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
366 if(!PE) return;
3
←
Taking false branch
→
367
368 for(unsigned int i=0; i<PE->mult; i++){
4
←
Loop condition is true. Entering loop body
→
10
←
Loop condition is true. Entering loop body
→
16
←
Loop condition is true. Entering loop body
→
369 s_HitView_t *hits = PE->in[i].hitView;
370 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
5
←
Taking false branch
→
11
←
Taking false branch
→
17
←
Taking false branch
→
371
372 // If no CDC hits were produced by HDGeant, then we need to add
373 // the branches needed to the HDDM tree
374 if(hits->centralDC == HDDM_NULL(void*)&hddm_s_nullTarget){
6
←
Taking false branch
→
12
←
Taking false branch
→
18
←
Taking true branch
→
375 hits->centralDC = make_s_CentralDC();
376 hits->centralDC->cdcStraws = (s_CdcStraws_t*)HDDM_NULL(void*)&hddm_s_nullTarget;
377 hits->centralDC->cdcTruthPoints = (s_CdcTruthPoints_t*)HDDM_NULL(void*)&hddm_s_nullTarget;
378 }
379
380 if(hits->centralDC->cdcStraws == HDDM_NULL(void*)&hddm_s_nullTarget){
7
←
Taking false branch
→
13
←
Taking false branch
→
19
←
Taking false branch
→
381 hits->centralDC->cdcStraws = make_s_CdcStraws(0);
382 hits->centralDC->cdcStraws->mult=0;
383 }
384
385 // Get existing hits
386 s_CdcStraws_t *old_cdcstraws = hits->centralDC->cdcStraws;
387 unsigned int Nold = old_cdcstraws->mult;
388
389 // Create CdcStraws structure that has enough slots for
390 // both the real and noise hits.
391 s_CdcStraws_t* cdcstraws = make_s_CdcStraws((unsigned int)Nnoise_straws + Nold);
392
393 // Add real hits back in first
394 cdcstraws->mult = 0;
395 for(unsigned int j=0; j<Nold; j++){
8
←
Loop condition is false. Execution continues on line 406
→
14
←
Loop condition is false. Execution continues on line 406
→
20
←
Loop condition is false. Execution continues on line 406
→
396 cdcstraws->in[cdcstraws->mult++] = old_cdcstraws->in[j];
397
398 // We need to transfer ownership of the hits to the new cdcstraws
399 // branch so they don't get deleted when old_cdcstraws is freed.
400 s_CdcStraw_t *cdcstraw = &old_cdcstraws->in[j];
401 cdcstraw->cdcStrawHits = (s_CdcStrawHits_t *)HDDM_NULL(void*)&hddm_s_nullTarget;
402 }
403
404 // Delete memory used for old hits structure and
405 // replace pointer in HDDM tree with ours
406 free(old_cdcstraws);
21
←
Argument to free() is the address of the global variable 'hddm_s_nullTarget', which is not memory allocated by malloc()
407 hits->centralDC->cdcStraws = cdcstraws;
408
409 // Loop over straws with noise hits
410 for(unsigned int j=0; j<Nstraw_hits.size(); j++){
9
←
Loop condition is false. Execution continues on line 368
→
15
←
Loop condition is false. Execution continues on line 368
→
411 s_CdcStraw_t *cdcstraw = &cdcstraws->in[cdcstraws->mult++];
412 s_CdcStrawHits_t *strawhits = make_s_CdcStrawHits(Nstraw_hits[j]);
413 cdcstraw->cdcStrawHits = strawhits;
414 cdcstraw->ring = ring_number[j];
415 cdcstraw->straw = straw_number[j];
416
417 strawhits->mult = 0;
418 for(int k=0; k<Nstraw_hits[j]; k++){
419 s_CdcStrawHit_t *strawhit = &strawhits->in[strawhits->mult++];
420 strawhit->dE = 1.0;
421 strawhit->t = SampleRange(-CDC_TIME_WINDOW/2.0, +CDC_TIME_WINDOW/2.0)*1.e9;
422 strawhit->d = 0.; // be consistent to initialize d=DOCA to zero
423 }
424 }
425 }
426}
427
428//-----------
429// SmearFDC
430//-----------
431void SmearFDC(s_HDDM_t *hddm_s)
432{
433
434 // Calculate ped noise level based on position resolution
435 // FDC_PED_NOISE=-0.004594+0.008711*FDC_CATHODE_SIGMA+0.000010*FDC_CATHODE_SIGMA*FDC_CATHODE_SIGMA; //pC
436 FDC_PED_NOISE=-0.0938+0.0485*FDC_CATHODE_SIGMA;
437 if(FDC_ELOSS_OFF)FDC_PED_NOISE*=7.0; // empirical 4/29/2009 DL
438
439 // Loop over Physics Events
440 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
441 if(!PE) return;
442
443 for(unsigned int i=0; i<PE->mult; i++){
444 s_HitView_t *hits = PE->in[i].hitView;
445 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget ||
446 hits->forwardDC == HDDM_NULL(void*)&hddm_s_nullTarget ||
447 hits->forwardDC->fdcChambers == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
448
449 s_FdcChambers_t* fdcChambers = hits->forwardDC->fdcChambers;
450 s_FdcChamber_t *fdcChamber = fdcChambers->in;
451 for(unsigned int j=0; j<fdcChambers->mult; j++, fdcChamber++){
452
453 // Add pedestal noise to strip charge data
454 s_FdcCathodeStrips_t *strips= fdcChamber->fdcCathodeStrips;
455 if (strips!=HDDM_NULL(void*)&hddm_s_nullTarget){
456 s_FdcCathodeStrip_t *strip=strips->in;
457 for (unsigned int k=0;k<strips->mult;k++,strip++){
458
459 // If a s_FdcCathodeHits_t object already exists delete it
460 if(strip->fdcCathodeHits!=HDDM_NULL(void*)&hddm_s_nullTarget){
461 FREE(strip->fdcCathodeHits)free(strip->fdcCathodeHits);
462 strip->fdcCathodeHits = (s_FdcCathodeHits_t*)HDDM_NULL(void*)&hddm_s_nullTarget;
463 }
464
465 s_FdcCathodeTruthHits_t *truthhits=strip->fdcCathodeTruthHits;
466 if (truthhits==HDDM_NULL(void*)&hddm_s_nullTarget)continue;
467
468 // Allocate s_FdcCathodeHits_t objects corresponding to this s_FdcCathodeTruthHits_t
469 s_FdcCathodeHits_t *hits = strip->fdcCathodeHits = make_s_FdcCathodeHits(truthhits->mult);
470 hits->mult = truthhits->mult;
471
472 s_FdcCathodeHit_t *hit=hits->in;
473 s_FdcCathodeTruthHit_t *truthhit=truthhits->in;
474 for (unsigned int s=0;s<hits->mult;s++,hit++,truthhit++){
475 if(SampleRange(0.0, 1.0)<=FDC_HIT_DROP_FRACTION){
476 hit->q = 0.0;
477 hit->t = 1.0E6;
478 }else{
479 hit->q = truthhit->q + SampleGaussian(FDC_PED_NOISE);
480 hit->t = truthhit->t + SampleGaussian(FDC_TDRIFT_SIGMA)*1.0E9;;
481 }
482 }
483 }
484 }
485
486 // Add drift time varation to the anode data
487 s_FdcAnodeWires_t *wires=fdcChamber->fdcAnodeWires;
488
489 if (wires!=HDDM_NULL(void*)&hddm_s_nullTarget){
490 s_FdcAnodeWire_t *wire=wires->in;
491 for (unsigned int k=0;k<wires->mult;k++,wire++){
492
493 // If a s_FdcAnodeHits_t object exists already delete it
494 if(wire->fdcAnodeHits!=HDDM_NULL(void*)&hddm_s_nullTarget){
495 FREE(wire->fdcAnodeHits)free(wire->fdcAnodeHits);
496 wire->fdcAnodeHits = (s_FdcAnodeHits_t*)HDDM_NULL(void*)&hddm_s_nullTarget;
497 }
498
499 s_FdcAnodeTruthHits_t *truthhits=wire->fdcAnodeTruthHits;
500 if (truthhits==HDDM_NULL(void*)&hddm_s_nullTarget)continue;
501
502 // Allocate s_FdcAnodeHits_t object corresponding to this s_FdcAnodeTruthHits_t
503 s_FdcAnodeHits_t *hits = wire->fdcAnodeHits = make_s_FdcAnodeHits(truthhits->mult);
504 hits->mult = truthhits->mult;
505
506 s_FdcAnodeHit_t *hit=hits->in;
507 s_FdcAnodeTruthHit_t *truthhit=truthhits->in;
508 fdc_anode_mult->Fill(hits->mult);
509 for (unsigned int s=0;s<hits->mult;s++, hit++, truthhit++){
510 hit->t = truthhit->t + SampleGaussian(FDC_TDRIFT_SIGMA)*1.0E9;
511 hit->dE = truthhit->dE;
512 hit->d = 0.; // initialize d=DOCA to zero for consistency.
513
514 double dt=hit->t-truthhit->t_unsmeared-2.;
515 // Fill diagnostic histograms for first hit
516 if (s==0){
517 fdc_drift_time_smear_hist->Fill(truthhit->d,dt);
518 fdc_drift_dist_smear_hist->Fill(truthhit->d,dt*(0.5*0.02421/sqrt(truthhit->t_unsmeared)+5.09e-4));
519
520 fdc_drift_time->Fill(hit->t,truthhit->d);
521 }
522 }
523 }
524 }
525 }
526 }
527}
528
529//-----------
530// AddNoiseHitsFDC
531//-----------
532void AddNoiseHitsFDC(s_HDDM_t *hddm_s)
533{
534 if(!FDC_GEOMETRY_INITIALIZED)InitFDCGeometry();
535
536 // Calculate the number of noise hits for each FDC wire and store
537 // them in a sparse map. We must do it this way since we have to know
538 // the total number of s_FdcAnodeWire_t structures to allocate in our
539 // call to make_s_FdcAnodeWires.
540 //
541 // We do this using the individual wire rates to calculate the probability
542 // of the wire firing for a single event. For the FDC, we calculate the
543 // wire rates as a function of both wire number (distance from beam line)
544 // and layer (position in z). We want a roughly 1/r distribution in the
545 // radial direction and a roughly exponential rise in rate in the
546 // +z direction.
547 //
548 // The wire rates are obtained using a parameterization done
549 // to calculate the event size for the August 29, 2007 online
550 // meeting. This parameterization is almost already obsolete.
551 // In rough terms, the layer rate (integrated over all wires)
552 // is about 1 MHz. For a 24 layer chamber with a 1us time window,
553 // we should have approximately 24 background hits per event.
554 // 11/9/2007 D. L.
555 vector<int> Nwire_hits;
556 vector<int> wire_number;
557 vector<int> layer_number;
558 int Nnoise_wires = 0;
559 int Nnoise_hits = 0;
560 for(unsigned int layer=1; layer<=FDC_LAYER_Z.size(); layer++){
561 double No = FDC_RATE_COEFFICIENT*exp((double)layer*log(4.0)/24.0);
562 for(unsigned int wire=1; wire<=96; wire++){
563 double rwire = fabs(96.0/2.0 - (double)wire);
564 double N = No*log((rwire+0.5)/(rwire-0.5));
565
566 // Indivdual wire rates should be way less than 1/event so
567 // we just use the rate as a probablity.
568 double Nhits = SampleRange(0.0, 1.0)<N ? 1.0:0.0;
569 if(Nhits<1.0)continue;
570 int iNhits = (int)floor(Nhits);
571 Nwire_hits.push_back(iNhits);
572 wire_number.push_back(wire);
573 layer_number.push_back(layer);
574 Nnoise_wires++;
575 Nnoise_hits+=iNhits;
576 }
577 }
578
579 // Loop over Physics Events
580 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
581 if(!PE) return;
582
583 for(unsigned int i=0; i<PE->mult; i++){
584 s_HitView_t *hits = PE->in[i].hitView;
585 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
586
587 // If no FDC hits were produced by HDGeant, then we need to add
588 // the branches needed to the HDDM tree
589 if(hits->forwardDC == HDDM_NULL(void*)&hddm_s_nullTarget){
590 hits->forwardDC = make_s_ForwardDC();
591 hits->forwardDC->fdcChambers = (s_FdcChambers_t*)HDDM_NULL(void*)&hddm_s_nullTarget;
592 }
593
594 if(hits->forwardDC->fdcChambers == HDDM_NULL(void*)&hddm_s_nullTarget){
595 hits->forwardDC->fdcChambers = make_s_FdcChambers(0);
596 hits->forwardDC->fdcChambers->mult=0;
597 }
598
599 // Get existing hits
600 s_FdcChambers_t *old_fdcchambers = hits->forwardDC->fdcChambers;
601 unsigned int Nold = old_fdcchambers->mult;
602
603 // If we were doing this "right" we'd conglomerate all of the noise
604 // hits from the same chamber into the same s_FdcChamber_t structure.
605 // That's a pain in the butt and really may only save a tiny bit of disk
606 // space so we just add each noise hit back in as another chamber
607 // structure.
608
609
610 // Create FdcChambers structure that has enough slots for
611 // both the real and noise hits.
612 s_FdcChambers_t* fdcchambers = make_s_FdcChambers(Nwire_hits.size() + Nold);
613
614 // Add real hits back in first
615 fdcchambers->mult = 0;
616 for(unsigned int j=0; j<Nold; j++){
617 fdcchambers->in[fdcchambers->mult++] = old_fdcchambers->in[j];
618
619 // We need to transfer ownership of the hits to the new fdcchambers
620 // branch so they don't get deleted when old_fdcchambers is freed.
621 s_FdcChamber_t *fdcchamber = &old_fdcchambers->in[j];
622 fdcchamber->fdcAnodeWires = (s_FdcAnodeWires_t *)HDDM_NULL(void*)&hddm_s_nullTarget;
623 fdcchamber->fdcCathodeStrips = (s_FdcCathodeStrips_t *)HDDM_NULL(void*)&hddm_s_nullTarget;
624 }
625
626 // Delete memory used for old hits structure and
627 // replace pointer in HDDM tree with ours
628 free(old_fdcchambers);
629 hits->forwardDC->fdcChambers = fdcchambers;
630
631 // Loop over wires with noise hits
632 for(unsigned int j=0; j<Nwire_hits.size(); j++){
633 s_FdcChamber_t *fdcchamber = &fdcchambers->in[fdcchambers->mult++];
634
635 // Create structure for anode wires
636 s_FdcAnodeWires_t *fdcAnodeWires = make_s_FdcAnodeWires(Nwire_hits[j]);
637 fdcchamber->fdcAnodeWires = fdcAnodeWires;
638
639 fdcAnodeWires->mult = 0;
640
641 for(int k=0; k<Nwire_hits[j]; k++){
642 // Get pointer to anode wire structure
643 s_FdcAnodeWire_t *fdcAnodeWire = &fdcAnodeWires->in[fdcAnodeWires->mult++];
644
645 // Create anode hits structure
646 s_FdcAnodeHits_t *fdcanodehits = make_s_FdcAnodeHits(1);
647 fdcAnodeWire->fdcAnodeHits = fdcanodehits;
648
649 // Get pointer to anode hit structure
650 fdcanodehits->mult = 1;
651 s_FdcAnodeHit_t *fdcanodehit = &fdcanodehits->in[0];
652
653 fdcanodehit->dE = 0.1; // what should this be?
654 fdcanodehit->t = SampleRange(-FDC_TIME_WINDOW/2., +FDC_TIME_WINDOW/2.)*1.e9;
655 fdcanodehit->d = 0.; // d=DOCA initialize to avoid any NAN
656
657 fdcAnodeWire->wire = wire_number[j];
658
659 fdcchamber->layer = (layer_number[j]-1)%3 + 1;
660 fdcchamber->module = (layer_number[j]-1)/3 + 1;
661 }
662 }
663 }
664}
665
666//-----------
667// SmearFCAL
668//-----------
669void SmearFCAL(s_HDDM_t *hddm_s)
670{
671 /// Smear the FCAL hits using the nominal resolution of the individual blocks.
672 /// The way this works is a little funny and warrants a little explanation.
673 /// The information coming from hdgeant is truth information indexed by
674 /// row and column, but containing energy deposited and time. The mcsmear
675 /// program will copy the truth information from the FcalTruthBlock to the
676 /// FcalBlock branch, smearing the values with the appropriate detector
677 /// resolution.
678 ///
679 /// To access the "truth" values in DANA, get the DFCALHit objects using the
680 /// "TRUTH" tag.
681
682 // The code below is perhaps slightly odd in that it simultaneously creates
683 // and fills the mirror (truth) branch while smearing the values in the
684 // nominal hit branch.
685
686 // Loop over Physics Events
687 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
688 if(!PE) return;
689
690 if(!fcalGeom)fcalGeom = new DFCALGeometry();
691
692 for(unsigned int i=0; i<PE->mult; i++){
693 s_HitView_t *hits = PE->in[i].hitView;
694 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget ||
695 hits->forwardEMcal == HDDM_NULL(void*)&hddm_s_nullTarget ||
696 hits->forwardEMcal->fcalBlocks == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
697
698 s_FcalBlocks_t *blocks = hits->forwardEMcal->fcalBlocks;
699 for(unsigned int j=0; j<blocks->mult; j++){
700 s_FcalBlock_t *block = &blocks->in[j];
701
702 // Create FCAL hits structures to put smeared data into
703 if(block->fcalHits!=HDDM_NULL(void*)&hddm_s_nullTarget)free(block->fcalHits);
704 block->fcalHits = make_s_FcalHits(block->fcalTruthHits->mult);
705 block->fcalHits->mult = block->fcalTruthHits->mult;
706
707 for(unsigned int k=0; k<block->fcalTruthHits->mult; k++){
708 s_FcalTruthHit_t *fcaltruthhit = &block->fcalTruthHits->in[k];
709 s_FcalHit_t *fcalhit = &block->fcalHits->in[k];
710
711 // Copy info from truth stream before doing anything else
712 fcalhit->E = fcaltruthhit->E;
713 fcalhit->t = fcaltruthhit->t;
714
715 // Simulation simulates a grid of blocks for simplicity.
716 // Do not bother smearing inactive blocks. They will be
717 // discarded in DEventSourceHDDM.cc while being read in
718 // anyway.
719 if(!fcalGeom->isBlockActive( block->row, block->column ))continue;
720
721 // Smear the energy and timing of the hit
722 double sigma = FCAL_PHOT_STAT_COEF/sqrt(fcalhit->E) ;
723 fcalhit->E *= 1.0 + SampleGaussian(sigma);
724 fcalhit->t += SampleGaussian(200.0E-3); // smear by 200 ps fixed for now 7/2/2009 DL
725
726 // Apply a single block threshold. If the (smeared) energy is below this,
727 // then set the energy and time to zero.
728 if(fcalhit->E < FCAL_BLOCK_THRESHOLD){fcalhit->E = fcalhit->t = 0.0;}
729
730 } // k (fcalhits)
731 } // j (blocks)
732 } // i (physicsEvents)
733
734}
735
736//-----------
737// SmearCCAL
738//-----------
739void SmearCCAL(s_HDDM_t *hddm_s)
740{
741 /// Smear the CCAL hits using the same procedure as the FCAL above.
742 /// See those comments for details.
743
744 // Loop over Physics Events
745 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
746 if(!PE) return;
747
748 if(!ccalGeom)ccalGeom = new DCCALGeometry();
749
750 for(unsigned int i=0; i<PE->mult; i++){
751 s_HitView_t *hits = PE->in[i].hitView;
752 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget ||
753 hits->ComptonEMcal == HDDM_NULL(void*)&hddm_s_nullTarget ||
754 hits->ComptonEMcal->ccalBlocks == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
755
756 s_CcalBlocks_t *blocks = hits->ComptonEMcal->ccalBlocks;
757 for(unsigned int j=0; j<blocks->mult; j++){
758 s_CcalBlock_t *block = &blocks->in[j];
759
760 // Create FCAL hits structures to put smeared data into
761 if(block->ccalHits!=HDDM_NULL(void*)&hddm_s_nullTarget)free(block->ccalHits);
762 block->ccalHits = make_s_CcalHits(block->ccalTruthHits->mult);
763 block->ccalHits->mult = block->ccalTruthHits->mult;
764
765 for(unsigned int k=0; k<block->ccalTruthHits->mult; k++){
766 s_CcalTruthHit_t *ccaltruthhit = &block->ccalTruthHits->in[k];
767 s_CcalHit_t *ccalhit = &block->ccalHits->in[k];
768
769 // Copy info from truth stream before doing anything else
770 ccalhit->E = ccaltruthhit->E;
771 ccalhit->t = ccaltruthhit->t;
772
773 // Simulation simulates a grid of blocks for simplicity.
774 // Do not bother smearing inactive blocks. They will be
775 // discarded in DEventSourceHDDM.cc while being read in
776 // anyway.
777 if(!ccalGeom->isBlockActive( block->row, block->column ))continue;
778
779 // Smear the energy and timing of the hit
780 double sigma = CCAL_PHOT_STAT_COEF/sqrt(ccalhit->E) ;
781 ccalhit->E *= 1.0 + SampleGaussian(sigma);
782 ccalhit->t += SampleGaussian(200.0E-3); // smear by 200 ps fixed for now 7/2/2009 DL
783
784 // Apply a single block threshold. If the (smeared) energy is below this,
785 // then set the energy and time to zero.
786 if(ccalhit->E < CCAL_BLOCK_THRESHOLD){ccalhit->E = ccalhit->t = 0.0;}
787
788 } // k (fcalhits)
789 } // j (blocks)
790 } // i (physicsEvents)
791
792}
793
794//-----------
795// SmearTOF
796//-----------
797void SmearTOF(s_HDDM_t *hddm_s)
798{
799 // Loop over Physics Events
800 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
801 if(!PE) return;
802
803 for(unsigned int i=0; i<PE->mult; i++){
804 s_HitView_t *hits = PE->in[i].hitView;
805 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget ||
806 hits->forwardTOF == HDDM_NULL(void*)&hddm_s_nullTarget ||
807 hits->forwardTOF->ftofCounters == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
808
809
810 s_FtofCounters_t* ftofCounters = hits->forwardTOF->ftofCounters;
811
812 // Loop over counters
813
814 for(unsigned int j=0;j<ftofCounters->mult; j++){
815 s_FtofCounter_t *ftofCounter = &(ftofCounters->in[j]);
816
817 // take care of North Hits
818 s_FtofNorthTruthHits_t *ftofNorthTruthHits = ftofCounter->ftofNorthTruthHits;
819 ftofCounter->ftofNorthHits = make_s_FtofNorthHits(ftofNorthTruthHits->mult);
820 ftofCounter->ftofNorthHits->mult = ftofNorthTruthHits->mult;
821
822 for (unsigned int m=0;m<ftofNorthTruthHits->mult;m++){
823 s_FtofNorthTruthHit_t *ftofNorthTruthHit = &(ftofNorthTruthHits->in[m]);
824 s_FtofNorthHit_t *ftofHit = &(ftofCounter->ftofNorthHits->in[m]);
825
826 // Smear the time
827 ftofHit->t = ftofNorthTruthHit->t + SampleGaussian(TOF_SIGMA);
828
829 // Smear the energy
830 double npe = (double)ftofNorthTruthHit->dE * 1000. * TOF_PHOTONS_PERMEV;
831 npe = npe + SampleGaussian(sqrt(npe));
832 float NewE = npe/TOF_PHOTONS_PERMEV/1000.;
833 ftofHit->dE = NewE;
834 }
835
836 // take care of South Hits
837 s_FtofSouthTruthHits_t *ftofSouthTruthHits = ftofCounter->ftofSouthTruthHits;
838 ftofCounter->ftofSouthHits = make_s_FtofSouthHits(ftofSouthTruthHits->mult);
839 ftofCounter->ftofSouthHits->mult = ftofSouthTruthHits->mult;
840
841 for (unsigned int m=0;m<ftofSouthTruthHits->mult;m++){
842 s_FtofSouthTruthHit_t *ftofSouthTruthHit = &(ftofSouthTruthHits->in[m]);
843 s_FtofSouthHit_t *ftofHit = &(ftofCounter->ftofSouthHits->in[m]);
844
845 // Smear the time
846 ftofHit->t = ftofSouthTruthHit->t + SampleGaussian(TOF_SIGMA);
847
848 // Smear the energy
849 double npe = (double)ftofSouthTruthHit->dE * 1000. * TOF_PHOTONS_PERMEV;
850 npe = npe + SampleGaussian(sqrt(npe));
851 float NewE = npe/TOF_PHOTONS_PERMEV/1000.;
852 ftofHit->dE = NewE;
853
854 }
855 } // end loop over all counters
856 }
857}
858
859//-----------
860// SmearSTC - smear hits in the start counter
861//-----------
862void SmearSTC(s_HDDM_t *hddm_s)
863{
864
865 // Loop over Physics Events
866 s_PhysicsEvents_t* PE = hddm_s->physicsEvents;
867 if(!PE) return;
868
869 for(unsigned int i=0; i<PE->mult; i++){
870 s_HitView_t *hits = PE->in[i].hitView;
871 if (hits == HDDM_NULL(void*)&hddm_s_nullTarget ||
872 hits->startCntr == HDDM_NULL(void*)&hddm_s_nullTarget ||
873 hits->startCntr->stcPaddles == HDDM_NULL(void*)&hddm_s_nullTarget)continue;
874 s_StcPaddles_t *stcPaddles= hits->startCntr->stcPaddles;
875
876 // Loop over counters
877 for(unsigned int j=0;j<stcPaddles->mult; j++){
878 s_StcPaddle_t *stcPaddle = &(stcPaddles->in[j]);
879
880 s_StcTruthHits_t *stcTruthHits=stcPaddle->stcTruthHits;
881 stcPaddle->stcHits=make_s_StcHits(stcTruthHits->mult);
882 stcPaddle->stcHits->mult=stcTruthHits->mult;
883 for (unsigned int m=0;m<stcTruthHits->mult;m++){
884 s_StcTruthHit_t *stcTruthHit=&(stcTruthHits->in[m]);
885 s_StcHit_t *stcHit=&(stcPaddle->stcHits->in[m]);
886
887 // smear the time
888 stcHit->t=stcTruthHit->t + SampleGaussian(START_SIGMA);
889
890 // Smear the energy
891 double npe = (double)stcTruthHit->dE * 1000. * START_PHOTONS_PERMEV;
892 npe = npe + SampleGaussian(sqrt(npe));
893 float NewE = npe/START_PHOTONS_PERMEV/1000.;
894 stcHit->dE = NewE;
895 }
896 }
897 }
898}
899
900//-----------
901// SmearCherenkov
902//-----------
903void SmearCherenkov(s_HDDM_t *hddm_s)
904{
905
906
907
908}
909
910//-----------
911// InitCDCGeometry
912//-----------
913void InitCDCGeometry(void)
914{
915 CDC_GEOMETRY_INITIALIZED = true;
916
917 CDC_MAX_RINGS = 28;
918
919 //-- This was cut and pasted from DCDCTrackHit_factory.cc on 10/11/2007 --
920
921 float degrees0 = 0.0;
922 float degrees6 = 6.0*M_PI3.14159265358979323846/180.0;
923
924 for(int ring=1; ring<=CDC_MAX_RINGS; ring++){
925 int myNstraws=0;
926 float radius = 0.0;
927 float stereo=0.0;
928 float phi_shift=0.0;
929 float deltaX=0.0, deltaY=0.0;
930 float rotX=0.0, rotY=0.0;
931 switch(ring){
932 // axial
933 case 1: myNstraws= 42; radius= 10.7219; stereo= degrees0; phi_shift= 0.00000; break;
934 case 2: myNstraws= 42; radius= 12.097; stereo= degrees0; phi_shift= 4.285714; break;
935 case 3: myNstraws= 54; radius= 13.7803; stereo= degrees0; phi_shift= 2.00000; break;
936 case 4: myNstraws= 54; radius= 15.1621; stereo= degrees0; phi_shift= 5.3333333; break;
937
938 // -stereo
939 case 5: myNstraws= 66; radius= 16.9321; stereo= -degrees6; phi_shift= 0.33333; break;
940 case 6: myNstraws= 66; phi_shift= 0.33333; deltaX= 18.2948 ; deltaY= 0.871486; rotX=-6.47674; rotY=-0.302853; break;
941 case 7: myNstraws= 80; radius= 20.5213; stereo= -degrees6; phi_shift= -0.5000; break;
942 case 8: myNstraws= 80; phi_shift= -0.5000; deltaX= 21.8912; deltaY= 0.860106; rotX=-6.39548; rotY=-0.245615; break;
943
944 // +stereo
945 case 9: myNstraws= 93; radius= 23.8544; stereo= +degrees6; phi_shift= 1.1000; break;
946 case 10: myNstraws= 93; phi_shift= 1.1000; deltaX= 25.229; deltaY= 0.852573; rotX=+6.34142; rotY=+0.208647; break;
947 case 11: myNstraws= 106; radius= 27.1877; stereo= +degrees6; phi_shift= -1.40; break;
948 case 12: myNstraws= 106; phi_shift= -1.400; deltaX= 28.5658; deltaY= 0.846871; rotX=+6.30035; rotY=+0.181146; break;
949
950 // axial
951 case 13: myNstraws= 123; radius= 31.3799; stereo= degrees0; phi_shift= 0.5000000; break;
952 case 14: myNstraws= 123; radius= 32.7747; stereo= degrees0; phi_shift= 1.9634146; break;
953 case 15: myNstraws= 135; radius= 34.4343; stereo= degrees0; phi_shift= 1.0000000; break;
954 case 16: myNstraws= 135; radius= 35.8301; stereo= degrees0; phi_shift= 2.3333333; break;
955
956 // -stereo
957 case 17: myNstraws= 146; radius= 37.4446; stereo= -degrees6; phi_shift= 0.2; break;
958 case 18: myNstraws= 146; phi_shift= 0.2; deltaX= 38.8295; deltaY= 0.835653; rotX=-6.21919 ; rotY=-0.128247; break;
959 case 19: myNstraws= 158; radius= 40.5364; stereo= -degrees6; phi_shift= 0.7; break;
960 case 20: myNstraws= 158; phi_shift= 0.7; deltaX=41.9225 ; deltaY= 0.833676; rotX=-6.20274; rotY=-0.118271; break;
961
962 // +stereo
963 case 21: myNstraws= 170; radius= 43.6152; stereo= +degrees6; phi_shift= 1.1000; break;
964 case 22: myNstraws= 170; phi_shift= 1.1000; deltaX=45.0025 ; deltaY= 0.831738; rotX=+6.18859; rotY=+0.109325; break;
965 case 23: myNstraws= 182; radius= 46.6849; stereo= +degrees6; phi_shift= 1.40; break;
966 case 24: myNstraws= 182; phi_shift= 1.400; deltaX= 48.0733; deltaY= 0.829899; rotX=+6.1763; rotY=+0.101315; break;
967
968 // axial
969 case 25: myNstraws= 197; radius= 50.37; stereo= degrees0; phi_shift= 0.200000000; break;
970 case 26: myNstraws= 197; radius= 51.77; stereo= degrees0; phi_shift= 1.113705000; break;
971 case 27: myNstraws= 209; radius= 53.363; stereo= degrees0; phi_shift= 0.800000000; break;
972 case 28: myNstraws= 209; radius= 54.76; stereo= degrees0; phi_shift= 1.661244; break;
973 default:
974 cerr<<__FILE__"smear.cc"<<":"<<__LINE__974<<" Invalid value for CDC ring ("<<ring<<") should be 1-28 inclusive!"<<endl;
975 }
976 NCDC_STRAWS.push_back(myNstraws);
977 CDC_RING_RADIUS.push_back(radius);
978 }
979
980 double Nstraws = 0;
981 double alpha = 0.0;
982 for(unsigned int i=0; i<NCDC_STRAWS.size(); i++){
983 Nstraws += (double)NCDC_STRAWS[i];
984 alpha += (double)NCDC_STRAWS[i]/CDC_RING_RADIUS[i];
985 }
986}
987
988
989//-----------
990// InitFDCGeometry
991//-----------
992void InitFDCGeometry(void)
993{
994 FDC_GEOMETRY_INITIALIZED = true;
995
996 int FDC_NUM_LAYERS = 24;
997 //int WIRES_PER_PLANE = 96;
998 //int WIRE_SPACING = 1.116;
999
1000 for(int layer=1; layer<=FDC_NUM_LAYERS; layer++){
1001
1002 float degrees00 = 0.0;
1003 float degrees60 = M_PI3.14159265358979323846*60.0/180.0;
1004
1005 float angle=0.0;
1006 float z_anode=212.0+95.5;
1007 switch(layer){
1008 case 1: z_anode+= -92.5-2.0; angle= degrees00; break;
1009 case 2: z_anode+= -92.5+0.0; angle= +degrees60; break;
1010 case 3: z_anode+= -92.5+2.0; angle= -degrees60; break;
1011 case 4: z_anode+= -86.5-2.0; angle= degrees00; break;
1012 case 5: z_anode+= -86.5+0.0; angle= +degrees60; break;
1013 case 6: z_anode+= -86.5+2.0; angle= -degrees60; break;
1014
1015 case 7: z_anode+= -32.5-2.0; angle= degrees00; break;
1016 case 8: z_anode+= -32.5+0.0; angle= +degrees60; break;
1017 case 9: z_anode+= -32.5+2.0; angle= -degrees60; break;
1018 case 10: z_anode+= -26.5-2.0; angle= degrees00; break;
1019 case 11: z_anode+= -26.5+0.0; angle= +degrees60; break;
1020 case 12: z_anode+= -26.5+2.0; angle= -degrees60; break;
1021
1022 case 13: z_anode+= +26.5-2.0; angle= degrees00; break;
1023 case 14: z_anode+= +26.5+0.0; angle= +degrees60; break;
1024 case 15: z_anode+= +26.5+2.0; angle= -degrees60; break;
1025 case 16: z_anode+= +32.5-2.0; angle= degrees00; break;
1026 case 17: z_anode+= +32.5+0.0; angle= +degrees60; break;
1027 case 18: z_anode+= +32.5+2.0; angle= -degrees60; break;
1028
1029 case 19: z_anode+= +86.5-2.0; angle= degrees00; break;
1030 case 20: z_anode+= +86.5+0.0; angle= +degrees60; break;
1031 case 21: z_anode+= +86.5+2.0; angle= -degrees60; break;
1032 case 22: z_anode+= +92.5-2.0; angle= degrees00; break;
1033 case 23: z_anode+= +92.5+0.0; angle= +degrees60; break;
1034 case 24: z_anode+= +92.5+2.0; angle= -degrees60; break;
1035 }
1036
1037 FDC_LAYER_Z.push_back(z_anode);
1038 }
1039
1040 // Coefficient used to calculate FDCsingle wire rate. We calculate
1041 // it once here just to save calculating it for every wire in every event
1042 FDC_RATE_COEFFICIENT = exp(-log(4.0)/23.0)/2.0/log(24.0)*FDC_TIME_WINDOW/1000.0E-9;
1043
1044 // Something is a little off in my calculation above so I scale it down via
1045 // an emprical factor:
1046 FDC_RATE_COEFFICIENT *= 0.353;
1047}
1048