/[PAMELA software]/tof/flight/ToFPatch/src/ToFPatch.cpp
ViewVC logotype

Contents of /tof/flight/ToFPatch/src/ToFPatch.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Thu Aug 30 17:11:33 2018 UTC (6 years, 3 months ago) by mayorov
Branch: MAIN
CVS Tags: v10REDr01, HEAD
ToF patch from Wolfgang. Duplicates the simple method of the 8th reduction to calculate the ToF dEdx

1 #include <ToFPatch.h>
2 #include <TGraph.h>
3 #include <TSpline.h>
4 #include <TMVA/TSpline1.h>
5
6 //------------------------------------------------------------------------
7 void ToFPatch::Clear(Option_t *option)
8 {
9 //
10 L2 = 0;
11 OBT = 0;
12 PKT = 0;
13 atime = 0;
14 tr = 0;
15 //sntr = 0;
16 // Set arrays and initialize structure
17
18 eDEDXpmt.Set(48); eDEDXpmt.Reset(-1); // Set array size and reset structure
19 eDEDXpmtraw.Set(48); eDEDXpmtraw.Reset(-1); // Set array size and reset structure
20 eDEDXpad.Set(24); eDEDXpad.Reset(-1);
21 eDEDXlayer.Set(6); eDEDXlayer.Reset(-1);
22 INFOpmt.Set(48); INFOpmt.Reset(0);
23 INFOlayer.Set(6); INFOlayer.Reset(0);
24
25 eDEDXpmtstd.Set(48); eDEDXpmtstd.Reset(-1); // Set array size and reset structure
26 eDEDXpmtrawstd.Set(48); eDEDXpmtrawstd.Reset(-1); // Set array size and reset structure
27 eDEDXpadstd.Set(24); eDEDXpadstd.Reset(-1);
28 eDEDXlayerstd.Set(6); eDEDXlayerstd.Reset(-1);
29 INFOpmtstd.Set(48); INFOpmtstd.Reset(0);
30 INFOlayerstd.Set(6); INFOlayerstd.Reset(0);
31
32
33 //
34 };
35
36 //------------------------------------------------------------------------
37 void ToFPatch::Print(Option_t *option)
38 {
39 //
40 printf("========================================================================\n");
41 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
42 };
43
44
45 //------------------------------------------------------------------------
46 void ToFPatch::InitPar(const char *pardir, const char *param)
47 {
48 // expensive function - call it once/run
49 Define_PMTsat();
50 ReadParAtt(pardir, param);
51 };
52
53
54 //------------------------------------------------------------------------
55 void ToFPatch::Process( PamLevel2 *l2p, const char* alg)
56 {
57 // the parameters should be already initialised by InitPar()
58
59 Float_t thelp1,thelp2;
60 Float_t dedx_corr[48];
61 Float_t xv[6],yv[6],zin[6];
62
63 Clear();
64 L2 = l2p;
65 if(!L2) return;
66 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
67 TrkLevel2 *trkL2 = L2->GetTrkLevel2();
68 if(!trkL2) return;
69 ToFLevel2 *tofL2 = L2->GetToFLevel2();
70 if(!tofL2) return;
71 eNtr = trkL2->GetNTracks();
72 if(eNtr<1) return;
73
74
75 trkAlg = alg;
76
77 OBT = L2->GetOrbitalInfo()->OBT;
78 PKT = L2->GetOrbitalInfo()->pkt_num;
79 atime = L2->GetOrbitalInfo()->absTime;
80
81 //----------------------------------------------------------------------
82 //cout<<atime<<" "<<tmin_atten<<" "<<tmax_atten<<endl;
83
84 // tmin_atten & tmax_atten time interval limits for attenuation parameters
85 if ((atime < tmin_atten) || (atime > tmax_atten)) {
86 ical_atten = 0;
87 for (Int_t ii=0; ii<100; ii++) {
88 //cout<<ii<<" tmin "<<T_int_min[ii]<<" tmax "<<T_int_max[ii]<<endl;
89 if ((atime>T_int_min[ii]) && (atime<=T_int_max[ii])) { // time interval found
90 ical_atten = ii;
91 tmin_atten = T_int_min[ii];
92 tmax_atten = T_int_max[ii];
93 cout<<"Time Interval found: "<<ii<<" tmin "<<T_int_min[ii]<<" tmax "<<T_int_max[ii]<<endl;
94 break;
95 }
96 } // ii
97
98 TArrayF &Apar0 = A0_array[ical_atten];
99 TArrayF &Apar1 = A1_array[ical_atten];
100 TArrayF &Apar2 = A2_array[ical_atten];
101 TArrayF &Apar3 = A3_array[ical_atten];
102
103 for (Int_t ii=0; ii<48;ii++) {
104 A0[ii]=Apar0[ii];
105 A1[ii]=Apar1[ii];
106 A2[ii]=Apar2[ii];
107 A3[ii]=Apar3[ii];
108 // cout<<A0[ii]<<" "<<A1[ii]<<" "<<A2[ii]<<" "<<A3[ii]<<endl;
109 }
110
111
112 cout<<"attenuation corr: abs time "<<atime<<" interval "<<ical_atten<<" tmin "<<T_int_min[ical_atten]<<" tmax "<<T_int_max[ical_atten]<<endl;
113 }
114
115 //----------------------------------------------------------------------
116 // tmin_2nd & tmax_2nd time interval limits for 2nd order correction
117
118 //cout<<atime<<" "<<tmin_2nd<<" "<<tmax_2nd<<endl;
119
120 if ((atime < tmin_2nd) || (atime > tmax_2nd)) {
121
122 ical_2nd=0;
123 for (Int_t ii=0; ii<1500; ii++) {
124 //cout<<ii<<" "<<atime<<" tmin "<<mtime[ii]<<" tmax "<<mtime[ii+1]<<endl;
125 if ((atime>mtime[ii]) && (atime<=mtime[ii+1])) { // time interval found
126 ical_2nd = ii;
127 tmin_2nd = mtime[ii];
128 tmax_2nd = mtime[ii+1];
129 //cout<<"Time Interval 2nd order corr "<<ii<<" tmin "<<mtime[ii]<<" tmax "<<mtime[ii+1]<<endl;
130 break;
131 }
132 };
133 cout<<"2nd order corr: abs time "<<atime<<" interval "<<ical_2nd<<" tmin "<<mtime[ical_2nd]<<" tmax "<<mtime[ical_2nd+1]<<endl;
134 }
135 //-----------------------------------------------------------------------
136
137 //== interpolate betwen time limits
138
139 thelp1 = mtime[ical_2nd];
140 thelp2 = mtime[ical_2nd+1];
141
142 TArrayF &corr1 = dedx_corr_m[ical_2nd];
143 TArrayF &corr2 = dedx_corr_m[ical_2nd+1];
144
145 for (Int_t ii=0; ii<48;ii++) {
146 Float_t yhelp1 = corr1[ii];
147 Float_t yhelp2 = corr2[ii];
148 Float_t slope = (yhelp2-yhelp1)/(thelp2-thelp1);
149 Float_t inter = yhelp1 - slope*thelp1;
150 dedx_corr[ii] = slope*atime + inter;
151 // if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
152 // cout<<ii<<" "<<yhelp1<<" "<<yhelp2<<" => "<<dedx_corr[ii]<<endl;
153 }
154
155
156 //================================================================
157
158 if( L2->GetNTracks(trkAlg)>=1 ){
159
160
161 // PamTrack *trackTRK = L2->GetTrack(itr); // 9th
162 PamTrack *trackTRK = L2->GetTrack(0,trkAlg); // 10th
163
164
165 Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]);
166
167 if(betamean<0.05 || betamean>2){
168 for(int i=0;i<48;i++)INFOpmt[i]=1;
169 for(int i=0;i<6;i++)INFOlayer[i]=1;
170 }
171
172 // define angle:
173
174 for (Int_t jj=0; jj<6; jj++){
175 xv[jj]=trackTRK->GetToFTrack()->xtr_tof[jj];
176 yv[jj]=trackTRK->GetToFTrack()->ytr_tof[jj];
177 zin[jj] = L2->GetToFLevel2()->GetZTOF(L2->GetToFLevel2()->GetToFPlaneID(jj));
178 }
179
180 Float_t theta[6];
181
182 for (Int_t jj=0; jj<3; jj++){
183 Float_t dx=0.;
184 Float_t dy=0.;
185 Float_t dr=0.;
186 Int_t kk=2*jj;
187
188 theta[kk] =0;
189 theta[kk+1]=0;
190 Float_t dist = zin[kk] - zin[kk+1];
191
192 if ((xv[kk])<100.) {
193 dx = xv[kk] - xv[kk+1];
194 dy = yv[kk] - yv[kk+1];
195 dr = sqrt(dx*dx+dy*dy);
196 theta[kk] = atan(dr/dist);
197 theta[kk+1] = atan(dr/dist);
198 }
199
200 } // jj
201
202 // Clear ADC;
203 Float_t adc[48];
204 for(int i=0;i<48;i++){
205 // adc[i]=0;
206 adc[i]=4095;
207 }
208
209
210 for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
211 ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
212 Int_t pmt_id = tofpmt->pmt_id;
213 adc[pmt_id] = tofpmt->adc ; // raw ADC
214 } //loop on 48 pmts
215
216
217
218 Int_t paddleidoftrack[6]= {0} ;
219
220 for (Int_t ilay=0; ilay<6; ilay++) {
221 paddleidoftrack[ilay] = L2->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.2) ;
222 }
223
224 /*
225 cout<<"paddleidoftrack trk ";
226 for (Int_t ilay=0; ilay<6; ilay++) {
227 cout<<" "<<paddleidoftrack[ilay];
228 }
229 cout<<endl;
230
231
232 for (Int_t ilay=0; ilay<6; ilay++) {
233 cout<<" "<<theta[ilay];
234 }
235 cout<<endl;
236 */
237
238 //==================================================================
239 //==================================================================
240 //==================================================================
241
242 // Clear the dEdx arrays
243 for(int i=0;i<48;i++){
244 eDEDXpmtraw[i] = 1000;
245 eDEDXpmt[i] = 1000;
246 }
247
248 for(int i=0;i<24;i++) eDEDXpad[i] = 1000.;
249 for(int i=0;i<6;i++) eDEDXlayer[i] = 1000.;
250
251
252 // ratio helium to proton ca. 4
253 Float_t hepratio = 4.;
254
255 Int_t ipad,ihelp;
256 Float_t adclraw,adcrraw;
257 Float_t adcl,adcr;
258 Float_t dedxl,dedxr;
259 Float_t xkorr;
260
261 Int_t ipad_a[6] = {0,8,14,16,18,21};
262 Int_t ipmt_a[6] = {0,16,28,32,36,42};
263
264
265 for (Int_t ilay=0; ilay<6; ilay++) {
266
267 Int_t i1 = paddleidoftrack[ilay];
268 if (i1 != -1) {
269
270 ipad = ipad_a[ilay]+i1;
271 ihelp = ipmt_a[ilay]+(i1*2);
272
273 /*
274 if(adc[ihelp]==4095)INFOpmt[ihelp]=2;
275 else if(adc[iihelp]>=PMTsat[iihelp]-5)INFOpmt[iihelp]=3;
276
277 if( adc[ihelp] >= PMTsat[iihelp]-5 ) continue;
278 if( adc[ihelp] <= 0. ) continue;
279
280 double adcpC = f_adcPC( adc[iihelp] ); // - adc conversion in pC
281 double adccorr = adcpC*fabs(cos(theta[ilay]));
282
283 if(adccorr<=0)INFOpmt[ihelp]=4;
284 if(adccorr<=0.) continue;
285 */
286
287 adcl= adc[ihelp];
288 adcr= adc[ihelp+1];
289 adclraw = adcl;
290 adcrraw = adcr;
291 adcl = f_adcPC( adcl ); // - adc conversion in pC
292 adcr = f_adcPC( adcr ); // - adc conversion in pC
293
294 Float_t x;
295 if (ilay==0) x=yv[0];
296 if (ilay==1) x=xv[1];
297 if (ilay==2) x=xv[2];
298 if (ilay==3) x=yv[3];
299 if (ilay==4) x=yv[4];
300 if (ilay==5) x=xv[5];
301
302
303 if (adclraw<4095) {
304 adcl = adcl*cos(theta[ilay]);
305 xkorr = atten(A0[ihelp],A1[ihelp],A2[ihelp],A3[ihelp],x);
306 xkorr = xkorr/hepratio;
307 dedxl = adcl/xkorr;
308 eDEDXpmtraw[ihelp] = dedxl;
309 Float_t dedxl1 = dedxl*4./dedx_corr[ihelp]; // 2nd order corr
310 // cout<<A0[ihelp]<<" "<<A1[ihelp]<<" "<<A2[ihelp]<<" "<<A3[ihelp]<<" => "<<xkorr<<endl;
311 // cout<<"Layer "<<ilay<<" PMT "<<ihelp+1<<" ADClraw "<<adclraw<<" adcl "<<adcl<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxl<<" dedx_2nd "<<dedxl1<<endl;
312 eDEDXpmt[ihelp] = dedxl1;
313 // cout<<eDEDXpmt[ihelp]<<endl;
314 }
315
316 if (adcrraw<4095) {
317 adcr = adcr*cos(theta[ilay]);
318 xkorr = atten(A0[ihelp+1],A1[ihelp+1],A2[ihelp+1],A3[ihelp+1],x);
319 // cout<<A0[ihelp+1]<<" "<<A1[ihelp+1]<<" "<<A2[ihelp+1]<<" "<<A3[ihelp+1]<<" => "<<xkorr<<endl;
320 xkorr = xkorr/hepratio;
321 //cout<<adcrraw<<" "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr"<<xkorr<<endl;
322 dedxr = adcr/xkorr;
323 eDEDXpmtraw[ihelp+1] = dedxr;
324 Float_t dedxr1 = dedxr*4./dedx_corr[ihelp+1]; // 2nd order corr
325 // cout<<"Layer "<<ilay<<" PMT "<<ihelp+1+1<<" ADCrraw "<<adcrraw<<" adcr "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxr<<" dedx_2nd "<<dedxr1<<endl;
326 eDEDXpmt[ihelp+1] = dedxr1;
327 // cout<<eDEDXpmt[ihelp+1]<<endl;
328 }
329
330
331 } // i1>-1
332 } // ilay
333
334 //-----------------------------------------------------------------
335 // The GetdEdx method "handmade"
336 //-----------------------------------------------------------------
337
338
339 for (Int_t ilay=0; ilay<6; ilay++) {
340 Float_t xhelp=1000.;
341 Int_t i1 = paddleidoftrack[ilay];
342 if (i1 != -1) {
343 ipad = ipad_a[ilay]+i1;
344 ihelp = ipmt_a[ilay]+(i1*2);
345 if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]==1000.)) xhelp = eDEDXpmt[ihelp];
346 if ((eDEDXpmt[ihelp+1]<1000.) && (eDEDXpmt[ihelp]==1000.)) xhelp = eDEDXpmt[ihelp+1];
347 if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmt[ihelp]+eDEDXpmt[ihelp+1]);
348 eDEDXpad[ipad] = xhelp;
349 eDEDXlayer[ilay] = xhelp;
350 }
351 } // ilay
352
353 } // if( L2->GetNTracks(trkAlg)>=1 )...
354
355
356 //================================================================
357 //================================================================
358 //============== standalone dEdx and beta =============
359 //================================================================
360 //================================================================
361
362 if (L2->GetToFStoredTrack(-1) ) {
363
364 // ToF track
365
366 ToFTrkVar *tof = L2->GetToFStoredTrack(-1);
367
368 Float_t betamean = tof->CalcBeta(10., 10., 20.);
369
370 if(betamean<0.05 || betamean>2){
371 for(int i=0;i<48;i++)INFOpmt[i]=1;
372 for(int i=0;i<6;i++)INFOlayer[i]=1;
373 }
374
375
376
377 // define angle (simple way... all theta[i] are the same...)
378
379 Float_t theta[6] = {0., 0., 0., 0., 0., 0. };
380
381 if (fabs((tof->xtofpos[0])<100) && (fabs(tof->xtofpos[2])<100) && (fabs(tof->ytofpos[0])<100) && (fabs(tof->ytofpos[2])<100)) {
382 Float_t dx = tof->xtofpos[0] - tof->xtofpos[2];
383 Float_t dy = tof->ytofpos[0] - tof->ytofpos[2];
384 Float_t dr = sqrt(dx*dx+dy*dy);
385 for (Int_t jj=0; jj<6; jj++) theta[jj]=atan(dr/77.31);
386 // cout<<tof->xtofpos[0]<<" "<<tof->xtofpos[2]<<" => "<<dx<<endl;
387 // cout<<tof->ytofpos[0]<<" "<<tof->ytofpos[2]<<" => "<<dy<<endl;
388 } // if...
389
390 /*
391 for (Int_t ilay=0; ilay<6; ilay++) {
392 cout<<" "<<theta[ilay];
393 }
394 cout<<endl;
395 */
396
397 // Clear ADC;
398 Float_t adc[48];
399 for(int i=0;i<48;i++){
400 // adc[i]=0;
401 adc[i]=4095;
402 }
403
404
405 for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
406 ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
407 Int_t pmt_id = tofpmt->pmt_id;
408 adc[pmt_id] = tofpmt->adc ; // raw ADC
409 } //loop on 48 pmts
410
411
412 //====================================================================
413 //=========== Check ToF standalone using HitPaddle ==================
414 //====================================================================
415
416 Int_t paddleidoftrack[6] = {-1, -1, -1, -1, -1, -1};;
417
418 for(Int_t jj=0; jj<8; jj++){
419 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(0,jj);
420 if (HitPad==1) paddleidoftrack[0] = jj;
421 }
422 for(Int_t jj=0; jj<6; jj++){
423 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(1,jj);
424 if (HitPad==1) paddleidoftrack[1] = jj;
425 }
426
427 for(Int_t jj=0; jj<2; jj++){
428 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(2,jj);
429 if (HitPad==1) paddleidoftrack[2] = jj;
430 }
431 for(Int_t jj=0; jj<2; jj++){
432 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(3,jj);
433 if (HitPad==1) paddleidoftrack[3] = jj;
434 }
435
436 for(Int_t jj=0; jj<3; jj++){
437 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(4,jj);
438 if (HitPad==1) paddleidoftrack[4] = jj;
439 }
440 for(Int_t jj=0; jj<3; jj++){
441 Int_t HitPad = L2->GetToFLevel2()->HitPaddle(5,jj);
442 if (HitPad==1) paddleidoftrack[5] = jj;
443 }
444
445 /*
446 cout<<"paddleidoftrack std ";
447 for (Int_t ilay=0; ilay<6; ilay++) {
448 cout<<" "<<paddleidoftrack[ilay];
449 }
450 cout<<endl;
451 */
452 //==================================================================
453 //==================================================================
454 //==================================================================
455
456 // Clear the dEdx arrays
457 for(int i=0;i<48;i++){
458 eDEDXpmtrawstd[i] = 1000;
459 eDEDXpmtstd[i] = 1000;
460 }
461
462 for(int i=0;i<24;i++) eDEDXpadstd[i] = 1000.;
463 for(int i=0;i<6;i++) eDEDXlayerstd[i] = 1000.;
464
465
466 // ratio helium to proton ca. 4
467 Float_t hepratio = 4.;
468
469 Int_t ipad,ihelp;
470 Float_t adclraw,adcrraw;
471 Float_t adcl,adcr;
472 Float_t dedxl,dedxr;
473 Float_t xkorr;
474
475 Int_t ipad_a[6] = {0,8,14,16,18,21};
476 Int_t ipmt_a[6] = {0,16,28,32,36,42};
477
478
479 for (Int_t ilay=0; ilay<6; ilay++) {
480
481 Int_t i1 = paddleidoftrack[ilay];
482 if (i1 != -1) {
483
484 ipad = ipad_a[ilay]+i1;
485 ihelp = ipmt_a[ilay]+(i1*2);
486
487 /*
488 if(adc[ihelp]==4095)INFOpmt[ihelp]=2;
489 else if(adc[iihelp]>=PMTsat[iihelp]-5)INFOpmt[iihelp]=3;
490
491 if( adc[ihelp] >= PMTsat[iihelp]-5 ) continue;
492 if( adc[ihelp] <= 0. ) continue;
493
494 double adcpC = f_adcPC( adc[iihelp] ); // - adc conversion in pC
495 double adccorr = adcpC*fabs(cos(theta[ilay]));
496
497 if(adccorr<=0)INFOpmt[ihelp]=4;
498 if(adccorr<=0.) continue;
499 */
500
501 adcl= adc[ihelp];
502 adcr= adc[ihelp+1];
503 adclraw = adcl;
504 adcrraw = adcr;
505 adcl = f_adcPC( adcl ); // - adc conversion in pC
506 adcr = f_adcPC( adcr ); // - adc conversion in pC
507
508 Float_t x;
509 if (ilay==0) x = tof->ytofpos[0];
510 if (ilay==1) x = tof->xtofpos[0];
511 if (ilay==2) x = tof->xtofpos[1];
512 if (ilay==3) x = tof->ytofpos[1];
513 if (ilay==4) x = tof->ytofpos[2];
514 if (ilay==5) x = tof->xtofpos[2];
515
516 /*
517 Float_t xx;
518 if (ilay==0) xx=yv[0];
519 if (ilay==1) xx=xv[1];
520 if (ilay==2) xx=xv[2];
521 if (ilay==3) xx=yv[3];
522 if (ilay==4) xx=yv[4];
523 if (ilay==5) xx=xv[5];
524
525 cout<<ilay<<" xv "<<xx<<" tof std "<<x<<endl;
526 */
527
528 if (fabs(x) < 100.) {
529
530 if (adclraw<4095) {
531 adcl = adcl*cos(theta[ilay]);
532 xkorr = atten(A0[ihelp],A1[ihelp],A2[ihelp],A3[ihelp],x);
533 xkorr = xkorr/hepratio;
534 dedxl = adcl/xkorr;
535 eDEDXpmtrawstd[ihelp] = dedxl;
536 Float_t dedxl1 = dedxl*4./dedx_corr[ihelp]; // 2nd order corr
537 // cout<<A0[ihelp]<<" "<<A1[ihelp]<<" "<<A2[ihelp]<<" "<<A3[ihelp]<<" => "<<xkorr<<endl;
538 // cout<<"Layer "<<ilay<<" PMT "<<ihelp+1<<" ADClraw "<<adclraw<<" adcl "<<adcl<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxl<<" dedx_2nd "<<dedxl1<<endl;
539 eDEDXpmtstd[ihelp] = dedxl1;
540 // cout<<"trk "<<eDEDXpmt[ihelp]<<" std "<<eDEDXpmtstd[ihelp]<<endl;
541 }
542
543 if (adcrraw<4095) {
544 adcr = adcr*cos(theta[ilay]);
545 xkorr = atten(A0[ihelp+1],A1[ihelp+1],A2[ihelp+1],A3[ihelp+1],x);
546 // cout<<A0[ihelp+1]<<" "<<A1[ihelp+1]<<" "<<A2[ihelp+1]<<" "<<A3[ihelp+1]<<" => "<<xkorr<<endl;
547 xkorr = xkorr/hepratio;
548 //cout<<adcrraw<<" "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr"<<xkorr<<endl;
549 dedxr = adcr/xkorr;
550 eDEDXpmtrawstd[ihelp+1] = dedxr;
551 Float_t dedxr1 = dedxr*4./dedx_corr[ihelp+1]; // 2nd order corr
552 // cout<<"Layer "<<ilay<<" PMT "<<ihelp+1+1<<" ADCrraw "<<adcrraw<<" adcr "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxr<<" dedx_2nd "<<dedxr1<<endl;
553 eDEDXpmtstd[ihelp+1] = dedxr1;
554 // cout<<eDEDXpmt[ihelp+1]<<endl;
555 // cout<<"trk "<<eDEDXpmt[ihelp+1]<<" std "<<eDEDXpmtstd[ihelp+1]<<endl;
556 }
557 } // tofpos < 100
558
559 } // i1>-1
560 } // ilay
561
562 //-----------------------------------------------------------------
563 // The GetdEdx method "handmade"
564 //-----------------------------------------------------------------
565
566
567 for (Int_t ilay=0; ilay<6; ilay++) {
568 Float_t xhelp=1000.;
569 Int_t i1 = paddleidoftrack[ilay];
570 if (i1 != -1) {
571 ipad = ipad_a[ilay]+i1;
572 ihelp = ipmt_a[ilay]+(i1*2);
573 if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]==1000.)) xhelp = eDEDXpmtstd[ihelp];
574 if ((eDEDXpmtstd[ihelp+1]<1000.) && (eDEDXpmtstd[ihelp]==1000.)) xhelp = eDEDXpmtstd[ihelp+1];
575 if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmtstd[ihelp]+eDEDXpmtstd[ihelp+1]);
576 eDEDXpadstd[ipad] = xhelp;
577 eDEDXlayerstd[ilay] = xhelp;
578 }
579 } // ilay
580
581 } // if (pam_event->GetToFStoredTrack(-1) )...
582
583
584 };
585
586 //------------------------------------------------------------------------
587 void ToFPatch::Define_PMTsat()
588 {
589 Float_t sat[48] = {
590 3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
591 3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
592 3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
593 3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
594 3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
595 3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
596 PMTsat.Set(48,sat);
597 }
598
599
600 //------------------------------------------------------------------------
601 void ToFPatch::ReadParAtt(const char *pardir, const char *param)
602 {
603 char filename[100],filename1[100],filename2[100],name[50];
604 char str1[10];
605
606
607
608 //--------------------- flight ---------------------------
609
610 if (strcmp(param,"simu") != 0){
611 //cout<<"Flight calibrations "<<endl;
612
613 sprintf(filename,"%stime_intervals_10th_red.txt",pardir);
614
615 //cout<<"read Time Interval File "<<filename<<endl;
616
617 FILE *ftimein = fopen( filename , "r" );
618
619 for (int ii=0; ii<100; ii++) {
620 UInt_t tp[2];
621 if(fscanf(ftimein,"%s %u %u",
622 &str1[0], &tp[0], &tp[1])!=3) break;
623 //cout<<str1<<" "<<tp[0]<<" "<<tp[1]<<endl;
624 T_int_min[ii] = tp[0];
625 T_int_max[ii] = tp[1];
626
627
628 sprintf(name,"tofcalib_%s_atten.dat",str1);
629 // cout<<"test name "<<name<<endl;
630 sprintf(filename1,"%s%s",pardir,name);
631
632 //cout<<"read Attenuation file "<<filename1<<endl;
633 FILE *fattin = fopen( filename1, "r" );
634
635 Float_t A0help[48];
636 Float_t A1help[48];
637 Float_t A2help[48];
638 Float_t A3help[48];
639
640
641 for (int i=0; i<48; i++) {
642 int id=0;
643 float par[4];
644 if(fscanf(fattin,"%d %f %f %f %f",
645 &id, &par[0], &par[1], &par[2], &par[3] )!=5) break;
646 A0help[i] = par[0];
647 A1help[i] = par[1];
648 A2help[i] = par[2];
649 A3help[i] = par[3];
650 }
651
652 A0_array[ii].Set(48,A0help);
653 A1_array[ii].Set(48,A1help);
654 A2_array[ii].Set(48,A2help);
655 A3_array[ii].Set(48,A3help);
656
657
658 fclose(fattin);
659
660
661 if (T_int_max[ii] == 2000000000) break; // end of file
662
663 }
664 fclose(ftimein);
665
666 } // flight
667
668 //--------------------- simu ---------------------------
669
670 if ( strcmp(param,"simu") == 0){
671
672 //cout<<"Simulation calibrations "<<endl;
673 sprintf(name,"tofcalib_gp_atten.dat");
674 sprintf(filename1,"%s%s",pardir,name);
675
676 cout<<"read Attenuation file "<<filename1<<endl;
677 FILE *fattin = fopen( filename1, "r" );
678
679 Float_t A0help[48];
680 Float_t A1help[48];
681 Float_t A2help[48];
682 Float_t A3help[48];
683
684
685 for (int i=0; i<48; i++) {
686 int id=0;
687 float par[4];
688 if(fscanf(fattin,"%d %f %f %f %f",
689 &id, &par[0], &par[1], &par[2], &par[3] )!=5) break;
690 A0help[i] = par[0];
691 A1help[i] = par[1];
692 A2help[i] = par[2];
693 A3help[i] = par[3];
694 }
695
696 A0_array[0].Set(48,A0help);
697 A1_array[0].Set(48,A1help);
698 A2_array[0].Set(48,A2help);
699 A3_array[0].Set(48,A3help);
700
701
702 fclose(fattin);
703
704 T_int_min[0] = 1000000000;
705 T_int_max[0] = 2000000000;
706
707
708 } // simu
709
710
711
712
713 // set first time interval
714 ical_atten = 0;
715 tmin_atten = T_int_min[0];
716 tmax_atten = T_int_max[0];
717
718 // set first time interval attenuation
719
720 TArrayF &Apar0 = A0_array[0];
721 TArrayF &Apar1 = A1_array[0];
722 TArrayF &Apar2 = A2_array[0];
723 TArrayF &Apar3 = A3_array[0];
724
725 for (Int_t ii=0; ii<48;ii++) {
726 A0[ii]=Apar0[ii];
727 A1[ii]=Apar1[ii];
728 A2[ii]=Apar2[ii];
729 A3[ii]=Apar3[ii];
730 }
731
732 cout<<"First time interval attenuation: "<<tmin_atten<<" - "<<tmax_atten<<endl;
733
734
735 //-----------------------------------------------------------
736 // Here I read the dEdx_korr parameters 2nd order
737 //-----------------------------------------------------------
738
739 //Double_t t1,t2,tm;
740 UInt_t t1,t2,tm;
741 Float_t xmean1,xwidth1;
742 Float_t pmthelp[48];
743
744 //--------------------- flight ---------------------------
745 if ( strcmp(param,"simu") != 0){
746
747 sprintf(filename2,"%sdedx_2nd_corr_paddlehit.dat",pardir);
748
749 //cout<<"Filename 2nd order "<<filename2<<endl;
750 ifstream fin(filename2);
751
752 Int_t jj=0;
753 Int_t idummy;
754 while (! fin.eof()) {
755 fin>>t1>>tm>>t2;
756 //cout<<jj<<" "<<tm<<endl;
757 //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
758 //cout <<t1<<" "<<tm<<" "<<t2<<endl;
759 mtime[jj]=tm;
760 for (Int_t ii=0; ii<48;ii++) {
761 fin>>idummy>>xmean1>>xwidth1;
762 pmthelp[ii] = xmean1;
763 }
764 dedx_corr_m[jj].Set(48,pmthelp);
765 jj=jj+1;
766 }
767
768 /*
769 Int_t idummy;
770 for (Int_t jj=0; jj<2000; jj++) {
771 fin>>t1>>tm>>t2;
772 cout<<jj<<" "<<t1<<" "<<tm<<" "<<t2<<endl;
773 if (t2==2000000000) break;
774 cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
775 mtime[jj]=tm;
776 for (Int_t ii=0; ii<48;ii++) {
777 fin>>idummy>>xmean1>>xwidth1;
778 dedx_corr_m[jj][ii]=xmean1;
779 }
780 }
781 */
782
783 fin.close();
784
785 } // flight
786
787 //--------------------- simu ---------------------------
788
789 if ( strcmp(param,"simu") == 0){
790
791 sprintf(filename2,"%sdedx_2nd_corr_simu.dat",pardir);
792
793 cout<<"Filename 2nd order "<<filename2<<endl;
794 ifstream fin(filename2);
795
796 Int_t jj=0;
797 Int_t idummy;
798 while (! fin.eof()) {
799 fin>>t1>>tm>>t2;
800 //cout<<jj<<" "<<tm<<endl;
801 //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
802 cout <<t1<<" "<<tm<<" "<<t2<<endl;
803 mtime[jj]=tm;
804 for (Int_t ii=0; ii<48;ii++) {
805 fin>>idummy>>xmean1>>xwidth1;
806 pmthelp[ii] = xmean1;
807 }
808 dedx_corr_m[jj].Set(48,pmthelp);
809 jj=jj+1;
810 }
811
812 fin.close();
813
814 } // simu
815
816
817 // set first time interval
818 ical_2nd=0;
819 tmin_2nd = mtime[0];
820 tmax_2nd = mtime[1];
821
822 cout<<"First time interval 2nd order corr.: "<<tmin_2nd<<" - "<<tmax_2nd<<endl;
823
824
825 }
826
827
828 //----------------------------------------------------------------------------
829
830 double ToFPatch::atten(float C0, float C1, float C2, float C3, float x){
831 //
832 Float_t yval = C0*exp(x*C1)+C2*exp(x*C3) ;
833 // cout<<" in fl "<<C0<<" " <<C1<<" "<<C2<<" "<<C3<<" "<<C4<<" "<<x<<" "<<yval<<endl;
834 return yval ;
835 }
836
837
838 //------------------------------------------------------------------------
839
840 double ToFPatch::f_adcPC( float x )
841 {
842 // return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
843 return 28.0407 + 0.628929*x - 5.80901e-05*x*x + 3.14092e-08*x*x*x; // standard calibration curve since jan-07
844
845 }
846
847 //------------------------------------------------------------------------
848
849
850
851

  ViewVC Help
Powered by ViewVC 1.1.23