1 | |
2 | |
3 | |
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 | |
69 | extern bool ADD_NOISE; |
70 | |
71 | |
72 | extern bool SMEAR_HITS; |
73 | |
74 | |
75 | extern bool IGNORE_SEEDS; |
76 | |
77 | |
78 | |
79 | |
80 | extern double CDC_TDRIFT_SIGMA; |
81 | |
82 | |
83 | |
84 | |
85 | extern double CDC_TIME_WINDOW; |
86 | |
87 | |
88 | extern double CDC_PEDESTAL_SIGMA; |
89 | |
90 | |
91 | |
92 | |
93 | extern bool FDC_USE_PARAMETERIZED_SIGMA; |
94 | |
95 | |
96 | |
97 | |
98 | extern double FDC_TDRIFT_SIGMA; |
99 | |
100 | |
101 | |
102 | |
103 | |
104 | |
105 | extern double FDC_CATHODE_SIGMA; |
106 | |
107 | |
108 | |
109 | |
110 | extern double FDC_PED_NOISE; |
111 | |
112 | |
113 | |
114 | |
115 | |
116 | |
117 | extern bool FDC_ELOSS_OFF; |
118 | |
119 | |
120 | extern double FDC_TIME_WINDOW; |
121 | |
122 | |
123 | extern double FDC_HIT_DROP_FRACTION; |
124 | |
125 | |
126 | |
127 | extern double FCAL_PHOT_STAT_COEF; |
128 | |
129 | |
130 | extern double FCAL_BLOCK_THRESHOLD; |
131 | |
132 | |
133 | |
134 | double CCAL_PHOT_STAT_COEF = 0.035/2.0; |
135 | |
136 | |
137 | |
138 | double CCAL_BLOCK_THRESHOLD = 20.0*k_MeV; |
139 | |
140 | |
141 | |
142 | extern double TOF_SIGMA; |
143 | extern double TOF_PHOTONS_PERMEV; |
144 | |
145 | |
146 | extern double START_SIGMA ; |
147 | extern double START_PHOTONS_PERMEV; |
148 | |
149 | extern bool SMEAR_BCAL; |
150 | |
151 | |
152 | |
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++){ |
| 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 169 | |
|
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++){ |
| 7 | | Loop condition is true. Entering loop body | |
|
172 | for (i=1;i<=n-m;i++){ |
| 8 | | Loop condition is true. Entering loop body | |
|
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; |
| |
| 10 | | Memory is never released; potential leak of memory pointed to by 'd' |
|
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 | |
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 | |
211 | |
212 | void GetAndSetSeeds(s_HDDM_t *hddm_s) |
213 | { |
214 | |
215 | |
216 | |
217 | |
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 | |
229 | my_rand = pe->reactions->in[0].random = make_s_Random(); |
230 | |
231 | }else{ |
232 | if(!IGNORE_SEEDS){ |
233 | |
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 | |
239 | |
240 | if((seed1+seed2+seed3) != 0)gDRandom.SetSeeds(seed1, seed2, seed3); |
241 | } |
242 | } |
243 | |
244 | |
245 | gDRandom.GetSeeds(seed1, seed2, seed3); |
246 | |
247 | |
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 | |
255 | |
256 | void SmearCDC(s_HDDM_t *hddm_s) |
257 | { |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | s_PhysicsEvents_t* allEvents = hddm_s->physicsEvents; |
265 | if(!allEvents)return; |
266 | |
267 | for (unsigned int m=0; m < allEvents->mult; m++) { |
268 | |
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 | |
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 | |
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 | |
299 | |
300 | double delta_t = SampleGaussian(sigma_t)*1.0E9; |
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 | |
311 | strawhit->dE = strawtruthhit->dE+SampleGaussian(CDC_PEDESTAL_SIGMA); |
312 | |
313 | |
314 | strawhit->d = 0.; |
315 | strawhit->itrack = strawtruthhit->itrack; |
316 | strawhit->ptype = strawtruthhit->ptype; |
317 | |
318 | |
319 | } |
320 | } |
321 | } |
322 | } |
323 | |
324 | |
325 | |
326 | |
327 | void AddNoiseHitsCDC(s_HDDM_t *hddm_s) |
328 | { |
329 | if(!CDC_GEOMETRY_INITIALIZED)InitCDCGeometry(); |
330 | |
331 | |
332 | |
333 | |
334 | |
335 | |
336 | |
337 | |
338 | |
339 | |
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 | |
352 | |
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 | |
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 | |
373 | |
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 | |
386 | s_CdcStraws_t *old_cdcstraws = hits->centralDC->cdcStraws; |
387 | unsigned int Nold = old_cdcstraws->mult; |
388 | |
389 | |
390 | |
391 | s_CdcStraws_t* cdcstraws = make_s_CdcStraws((unsigned int)Nnoise_straws + Nold); |
392 | |
393 | |
394 | cdcstraws->mult = 0; |
395 | for(unsigned int j=0; j<Nold; j++){ |
396 | cdcstraws->in[cdcstraws->mult++] = old_cdcstraws->in[j]; |
397 | |
398 | |
399 | |
400 | s_CdcStraw_t *cdcstraw = &old_cdcstraws->in[j]; |
401 | cdcstraw->cdcStrawHits = (s_CdcStrawHits_t *)HDDM_NULL(void*)&hddm_s_nullTarget; |
402 | } |
403 | |
404 | |
405 | |
406 | free(old_cdcstraws); |
407 | hits->centralDC->cdcStraws = cdcstraws; |
408 | |
409 | |
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.; |
423 | } |
424 | } |
425 | } |
426 | } |
427 | |
428 | |
429 | |
430 | |
431 | void SmearFDC(s_HDDM_t *hddm_s) |
432 | { |
433 | |
434 | |
435 | |
436 | FDC_PED_NOISE=-0.0938+0.0485*FDC_CATHODE_SIGMA; |
437 | if(FDC_ELOSS_OFF)FDC_PED_NOISE*=7.0; |
438 | |
439 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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.; |
513 | |
514 | double dt=hit->t-truthhit->t_unsmeared-2.; |
515 | |
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 | |
531 | |
532 | void AddNoiseHitsFDC(s_HDDM_t *hddm_s) |
533 | { |
534 | if(!FDC_GEOMETRY_INITIALIZED)InitFDCGeometry(); |
535 | |
536 | |
537 | |
538 | |
539 | |
540 | |
541 | |
542 | |
543 | |
544 | |
545 | |
546 | |
547 | |
548 | |
549 | |
550 | |
551 | |
552 | |
553 | |
554 | |
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 | |
567 | |
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 | |
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 | |
588 | |
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 | |
600 | s_FdcChambers_t *old_fdcchambers = hits->forwardDC->fdcChambers; |
601 | unsigned int Nold = old_fdcchambers->mult; |
602 | |
603 | |
604 | |
605 | |
606 | |
607 | |
608 | |
609 | |
610 | |
611 | |
612 | s_FdcChambers_t* fdcchambers = make_s_FdcChambers(Nwire_hits.size() + Nold); |
613 | |
614 | |
615 | fdcchambers->mult = 0; |
616 | for(unsigned int j=0; j<Nold; j++){ |
617 | fdcchambers->in[fdcchambers->mult++] = old_fdcchambers->in[j]; |
618 | |
619 | |
620 | |
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 | |
627 | |
628 | free(old_fdcchambers); |
629 | hits->forwardDC->fdcChambers = fdcchambers; |
630 | |
631 | |
632 | for(unsigned int j=0; j<Nwire_hits.size(); j++){ |
633 | s_FdcChamber_t *fdcchamber = &fdcchambers->in[fdcchambers->mult++]; |
634 | |
635 | |
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 | |
643 | s_FdcAnodeWire_t *fdcAnodeWire = &fdcAnodeWires->in[fdcAnodeWires->mult++]; |
644 | |
645 | |
646 | s_FdcAnodeHits_t *fdcanodehits = make_s_FdcAnodeHits(1); |
647 | fdcAnodeWire->fdcAnodeHits = fdcanodehits; |
648 | |
649 | |
650 | fdcanodehits->mult = 1; |
651 | s_FdcAnodeHit_t *fdcanodehit = &fdcanodehits->in[0]; |
652 | |
653 | fdcanodehit->dE = 0.1; |
654 | fdcanodehit->t = SampleRange(-FDC_TIME_WINDOW/2., +FDC_TIME_WINDOW/2.)*1.e9; |
655 | fdcanodehit->d = 0.; |
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 | |
668 | |
669 | void SmearFCAL(s_HDDM_t *hddm_s) |
670 | { |
671 | |
672 | |
673 | |
674 | |
675 | |
676 | |
677 | |
678 | |
679 | |
680 | |
681 | |
682 | |
683 | |
684 | |
685 | |
686 | |
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 | |
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 | |
712 | fcalhit->E = fcaltruthhit->E; |
713 | fcalhit->t = fcaltruthhit->t; |
714 | |
715 | |
716 | |
717 | |
718 | |
719 | if(!fcalGeom->isBlockActive( block->row, block->column ))continue; |
720 | |
721 | |
722 | double sigma = FCAL_PHOT_STAT_COEF/sqrt(fcalhit->E) ; |
723 | fcalhit->E *= 1.0 + SampleGaussian(sigma); |
724 | fcalhit->t += SampleGaussian(200.0E-3); |
725 | |
726 | |
727 | |
728 | if(fcalhit->E < FCAL_BLOCK_THRESHOLD){fcalhit->E = fcalhit->t = 0.0;} |
729 | |
730 | } |
731 | } |
732 | } |
733 | |
734 | } |
735 | |
736 | |
737 | |
738 | |
739 | void SmearCCAL(s_HDDM_t *hddm_s) |
740 | { |
741 | |
742 | |
743 | |
744 | |
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 | |
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 | |
770 | ccalhit->E = ccaltruthhit->E; |
771 | ccalhit->t = ccaltruthhit->t; |
772 | |
773 | |
774 | |
775 | |
776 | |
777 | if(!ccalGeom->isBlockActive( block->row, block->column ))continue; |
778 | |
779 | |
780 | double sigma = CCAL_PHOT_STAT_COEF/sqrt(ccalhit->E) ; |
781 | ccalhit->E *= 1.0 + SampleGaussian(sigma); |
782 | ccalhit->t += SampleGaussian(200.0E-3); |
783 | |
784 | |
785 | |
786 | if(ccalhit->E < CCAL_BLOCK_THRESHOLD){ccalhit->E = ccalhit->t = 0.0;} |
787 | |
788 | } |
789 | } |
790 | } |
791 | |
792 | } |
793 | |
794 | |
795 | |
796 | |
797 | void SmearTOF(s_HDDM_t *hddm_s) |
798 | { |
799 | |
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 | |
813 | |
814 | for(unsigned int j=0;j<ftofCounters->mult; j++){ |
815 | s_FtofCounter_t *ftofCounter = &(ftofCounters->in[j]); |
816 | |
817 | |
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 | |
827 | ftofHit->t = ftofNorthTruthHit->t + SampleGaussian(TOF_SIGMA); |
828 | |
829 | |
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 | |
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 | |
846 | ftofHit->t = ftofSouthTruthHit->t + SampleGaussian(TOF_SIGMA); |
847 | |
848 | |
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 | } |
856 | } |
857 | } |
858 | |
859 | |
860 | |
861 | |
862 | void SmearSTC(s_HDDM_t *hddm_s) |
863 | { |
864 | |
865 | |
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 | |
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 | |
888 | stcHit->t=stcTruthHit->t + SampleGaussian(START_SIGMA); |
889 | |
890 | |
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 | |
902 | |
903 | void SmearCherenkov(s_HDDM_t *hddm_s) |
904 | { |
905 | |
906 | |
907 | |
908 | } |
909 | |
910 | |
911 | |
912 | |
913 | void InitCDCGeometry(void) |
914 | { |
915 | CDC_GEOMETRY_INITIALIZED = true; |
916 | |
917 | CDC_MAX_RINGS = 28; |
918 | |
919 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
991 | |
992 | void InitFDCGeometry(void) |
993 | { |
994 | FDC_GEOMETRY_INITIALIZED = true; |
995 | |
996 | int FDC_NUM_LAYERS = 24; |
997 | |
998 | |
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 | |
1041 | |
1042 | FDC_RATE_COEFFICIENT = exp(-log(4.0)/23.0)/2.0/log(24.0)*FDC_TIME_WINDOW/1000.0E-9; |
1043 | |
1044 | |
1045 | |
1046 | FDC_RATE_COEFFICIENT *= 0.353; |
1047 | } |
1048 | |