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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide 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 mayorov 1.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