File: | programs/Simulation/mcsmear/smear.cc |
Location: | line 946, column 65 |
Description: | Value stored to 'deltaY' is never read |
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> |
8 | using 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 | |
28 | extern vector<vector<float> >fdc_smear_parms; |
29 | extern TH2F *fdc_drift_time_smear_hist; |
30 | extern TH2F *fdc_drift_dist_smear_hist; |
31 | extern TH2F *fdc_drift_time; |
32 | extern TH2F *cdc_drift_time; |
33 | extern TH1F *fdc_anode_mult; |
34 | extern TH2F *cdc_drift_smear; |
35 | |
36 | void GetAndSetSeeds(s_HDDM_t *hddm_s); |
37 | void SmearCDC(s_HDDM_t *hddm_s); |
38 | void AddNoiseHitsCDC(s_HDDM_t *hddm_s); |
39 | void SmearFDC(s_HDDM_t *hddm_s); |
40 | void AddNoiseHitsFDC(s_HDDM_t *hddm_s); |
41 | void SmearFCAL(s_HDDM_t *hddm_s); |
42 | void SmearCCAL(s_HDDM_t *hddm_s); |
43 | void SmearBCAL(s_HDDM_t *hddm_s); |
44 | void SmearTOF(s_HDDM_t *hddm_s); |
45 | void SmearSTC(s_HDDM_t *hddm_s); |
46 | void SmearCherenkov(s_HDDM_t *hddm_s); |
47 | void InitCDCGeometry(void); |
48 | void InitFDCGeometry(void); |
49 | |
50 | pthread_mutex_t mutex_fdc_smear_function = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
51 | |
52 | bool CDC_GEOMETRY_INITIALIZED = false; |
53 | int CDC_MAX_RINGS=0; |
54 | vector<unsigned int> NCDC_STRAWS; |
55 | vector<double> CDC_RING_RADIUS; |
56 | |
57 | DFCALGeometry *fcalGeom = NULL__null; |
58 | DCCALGeometry *ccalGeom = NULL__null; |
59 | bool FDC_GEOMETRY_INITIALIZED = false; |
60 | unsigned int NFDC_WIRES_PER_PLANE; |
61 | vector<double> FDC_LAYER_Z; |
62 | double FDC_RATE_COEFFICIENT; |
63 | |
64 | double SampleGaussian(double sigma); |
65 | double SamplePoisson(double lambda); |
66 | double SampleRange(double x1, double x2); |
67 | |
68 | // Do we or do we not add noise hits |
69 | extern bool ADD_NOISE; |
70 | |
71 | // Do we or do we not smear real hits |
72 | extern bool SMEAR_HITS; |
73 | |
74 | // Use or ignore random number seeds found in HDDM file |
75 | extern 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. |
80 | extern 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. |
85 | extern double CDC_TIME_WINDOW; |
86 | |
87 | // The error in the energy deposition measurement in the CDC due to pedestal noise |
88 | extern 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. |
93 | extern 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. |
98 | extern 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. |
105 | extern 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. |
110 | extern 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. |
117 | extern bool FDC_ELOSS_OFF; |
118 | |
119 | // Time window for acceptance of FDC hits |
120 | extern double FDC_TIME_WINDOW; |
121 | |
122 | // Fraction of FDC hits to randomly drop (0=drop nothing 1=drop everything) |
123 | extern 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) |
127 | extern double FCAL_PHOT_STAT_COEF; |
128 | |
129 | // Single block energy threshold (applied after smearing) |
130 | extern 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) |
134 | double 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) |
138 | double CCAL_BLOCK_THRESHOLD = 20.0*k_MeV; |
139 | |
140 | |
141 | // Forward TOF resolution |
142 | extern double TOF_SIGMA; |
143 | extern double TOF_PHOTONS_PERMEV; |
144 | |
145 | // Start counter resolution |
146 | extern double START_SIGMA ; |
147 | extern double START_PHOTONS_PERMEV; |
148 | |
149 | extern bool SMEAR_BCAL; |
150 | |
151 | // Polynomial interpolation on a grid. |
152 | // Adapted from Numerical Recipes in C (2nd Edition), pp. 121-122. |
153 | void 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 | //----------- |
193 | void 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 | //----------- |
212 | void 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 | //----------- |
256 | void 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 | //----------- |
327 | void AddNoiseHitsCDC(s_HDDM_t *hddm_s) |
328 | { |
329 | if(!CDC_GEOMETRY_INITIALIZED)InitCDCGeometry(); |
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++){ |
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; |
367 | |
368 | for(unsigned int i=0; i<PE->mult; i++){ |
369 | s_HitView_t *hits = PE->in[i].hitView; |
370 | if (hits == HDDM_NULL(void*)&hddm_s_nullTarget)continue; |
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){ |
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){ |
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++){ |
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); |
407 | hits->centralDC->cdcStraws = cdcstraws; |
408 | |
409 | // Loop over straws with noise hits |
410 | for(unsigned int j=0; j<Nstraw_hits.size(); j++){ |
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 | //----------- |
431 | void 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 | //----------- |
532 | void 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 | //----------- |
669 | void 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 | //----------- |
739 | void 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 | //----------- |
797 | void 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 | //----------- |
862 | void 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 | //----------- |
903 | void SmearCherenkov(s_HDDM_t *hddm_s) |
904 | { |
905 | |
906 | |
907 | |
908 | } |
909 | |
910 | //----------- |
911 | // InitCDCGeometry |
912 | //----------- |
913 | void 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; |
Value stored to 'deltaY' is never read | |
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 | //----------- |
992 | void 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 |