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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Mon Dec 10 18:47:02 2018 UTC (6 years ago) by mayorov
Branch: MAIN
CVS Tags: HEAD
ToFdEdx_patch from Wolfgang. Follows Rita Carbone's advanced method to
derive the ToF dEdx (including PMT saturation, Bethe-Bloch effect, etc.)

1 mayorov 1.1 #include <ToFdEdx_patch.h>
2     #include <TGraph.h>
3     #include <TSpline.h>
4     #include <TH1F.h>
5     #include <TMVA/TSpline1.h>
6    
7     //------------------------------------------------------------------------
8     void ToFdEdx_patch::Clear(Option_t *option)
9     {
10     //
11     L2 = 0;
12     OBT = 0;
13     PKT = 0;
14     atime = 0;
15     tr = 0;
16     //sntr = 0;
17     // Set arrays and initialize structure
18    
19     eDEDXpmt.Set(48); eDEDXpmt.Reset(-1); // Set array size and reset structure
20     eZpmt.Set(48); eZpmt.Reset(-1);
21     eDEDXpad.Set(24); eDEDXpad.Reset(-1);
22     eZpad.Set(24); eZpad.Reset(-1);
23     eDEDXlayer.Set(6); eDEDXlayer.Reset(-1);
24     eZlayer.Set(6); eZlayer.Reset(-1);
25     INFOpmt.Set(48); INFOpmt.Reset(0);
26     INFOlayer.Set(6); INFOlayer.Reset(0);
27     //
28    
29     // ToF standalone
30     eDEDXpmtstd.Set(48); eDEDXpmtstd.Reset(-1); // Set array size and reset structure
31     eZpmtstd.Set(48); eZpmtstd.Reset(-1);
32     eDEDXpadstd.Set(24); eDEDXpadstd.Reset(-1);
33     eZpadstd.Set(24); eZpadstd.Reset(-1);
34     eDEDXlayerstd.Set(6); eDEDXlayerstd.Reset(-1);
35     eZlayerstd.Set(6); eZlayerstd.Reset(-1);
36     INFOpmtstd.Set(48); INFOpmtstd.Reset(0);
37     INFOlayerstd.Set(6); INFOlayerstd.Reset(0);
38    
39     };
40    
41     //------------------------------------------------------------------------
42     void ToFdEdx_patch::Print(Option_t *option)
43     {
44     //
45     printf("========================================================================\n");
46     printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
47     };
48    
49    
50     //------------------------------------------------------------------------
51     void ToFdEdx_patch::InitPar(const char *pardir, const char *param)
52     {
53     // expensive function - call it once/run
54     Define_PMTsat();
55     // for(int i=0; i<48; i++) ReadParTD( i, Form("%s/TD_%d_200b.txt",pardir,i) );
56     // ReadParAtt( Form("%s/attenuation.txt" , pardir) );
57    
58     ReadParAtt(pardir, param);
59    
60     // Old ToFNaNuclei files
61     ReadParPos( Form("%s/desaturation_position.txt" , pardir) );
62     ReadParDesatBB( Form("%s/desaturation_beta.txt" , pardir) );
63     // ReadParBBneg( Form("%s/BetheBloch.txt" , pardir) );
64     // ReadParBBpos( Form("%s/BetheBloch_betagt1.txt" , pardir) );
65    
66     // New Bethe-Bloch pol4 file
67     ReadParBBpol4( Form("%s/BetheBloch_pol4_tri.txt" , pardir) );
68    
69    
70     };
71    
72    
73     //------------------------------------------------------------------------
74     void ToFdEdx_patch::Process( PamLevel2 *l2p, const char* alg, const char *tri_or_bi)
75     {
76     // the parameters should be already initialised by InitPar()
77    
78     Float_t dedx_corrp[48], dedx_corrhe[48], dedx_corrc[48];
79     Float_t thelp1,thelp2;
80     Float_t xv[6],yv[6]; //,zin[6];
81    
82     trkAlg = alg;
83     TF1 *ftri = new TF1("ftri", "1++x++x*x", 0,100);
84    
85     // cout<<"In PROCESS "<<endl;
86    
87     Clear();
88     L2 = l2p;
89     if(!L2) return;
90     if(!L2) cout<<" L2 error !!!"<<endl;
91    
92     // cout<<" L2->GetNTracks(trkAlg) check "<<L2->GetNTracks(trkAlg)<<endl;
93    
94    
95     // if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
96     TrkLevel2 *trkL2 = L2->GetTrkLevel2();
97     if(!trkL2) cout<<" L2->GetTrkLevel2() error !!!"<<endl;
98     if(!trkL2) return;
99     ToFLevel2 *tofL2 = L2->GetToFLevel2();
100     if(!tofL2) cout<<" L2->GetToFLevel2() error !!!"<<endl;
101     if(!tofL2) return;
102     // eNtr = trkL2->GetNTracks();
103     // if(eNtr<1) cout<<" eNtr < 1 !!!"<<endl;
104     // if(eNtr<1) return;
105    
106     eNtr = L2->GetNTracks(trkAlg);
107     if(eNtr<1) cout<<" eNtr < 1 !!!"<<endl;
108     if(eNtr<1) return;
109    
110    
111    
112     // cout<<"In PROCESS after first checks "<<endl;
113    
114    
115     // OBT = L2->GetOrbitalInfo()->OBT;
116     // PKT = L2->GetOrbitalInfo()->pkt_num;
117     // atime = L2->GetOrbitalInfo()->absTime;
118    
119     if ( L2->IsORB() ){
120     OBT = L2->GetOrbitalInfo()->OBT;
121     PKT = L2->GetOrbitalInfo()->pkt_num;
122     atime = L2->GetOrbitalInfo()->absTime;
123     } else {
124     atime = 1153000000.; // simulated data no OrbitalInfo fixed date 15-jul-2006
125     };
126    
127     //----------------------------------------------------------------------
128     //----------------------------------------------------------------------
129     // Check Attenuation correction file time interval
130     // tmin_atten & tmax_atten time interval limits for attenuation parameters
131     //----------------------------------------------------------------------
132    
133     //cout<<atime<<" "<<tmin_atten<<" "<<tmax_atten<<endl;
134    
135     if ((atime < tmin_atten) || (atime > tmax_atten)) {
136     ical_atten = 0;
137     for (Int_t ii=0; ii<100; ii++) {
138     //cout<<ii<<" tmin "<<T_int_min[ii]<<" tmax "<<T_int_max[ii]<<endl;
139     if ((atime>T_int_min[ii]) && (atime<=T_int_max[ii])) { // time interval found
140     ical_atten = ii;
141     tmin_atten = T_int_min[ii];
142     tmax_atten = T_int_max[ii];
143     cout<<"Time Interval found: "<<ii<<" tmin "<<T_int_min[ii]<<" tmax "<<T_int_max[ii]<<endl;
144     break;
145     }
146     } // ii
147    
148     TArrayF &Apar0 = A0_array[ical_atten];
149     TArrayF &Apar1 = A1_array[ical_atten];
150     TArrayF &Apar2 = A2_array[ical_atten];
151     TArrayF &Apar3 = A3_array[ical_atten];
152     TArrayF &Apar4 = A4_array[ical_atten];
153     TArrayF &Apar5 = A5_array[ical_atten];
154    
155     for (Int_t ii=0; ii<48;ii++) {
156     A0[ii]=Apar0[ii];
157     A1[ii]=Apar1[ii];
158     A2[ii]=Apar2[ii];
159     A3[ii]=Apar3[ii];
160     A4[ii]=Apar4[ii];
161     A5[ii]=Apar5[ii];
162     // cout<<A0[ii]<<" "<<A1[ii]<<" "<<A2[ii]<<" "<<A3[ii]<<endl;
163     }
164    
165     cout<<"attenuation corr: abs time "<<atime<<" interval "<<ical_atten<<" tmin "<<T_int_min[ical_atten]<<" tmax "<<T_int_max[ical_atten]<<endl;
166     }
167    
168     //----------------------------------------------------------------------
169     // tmin_2nd & tmax_2nd time interval limits for 2nd order correction
170     //----------------------------------------------------------------------
171    
172     //cout<<atime<<" "<<tmin_2nd<<" "<<tmax_2nd<<endl;
173    
174     if ((atime < tmin_2nd) || (atime > tmax_2nd)) {
175    
176     ical_2nd=0;
177     for (Int_t ii=0; ii<3500; ii++) {
178     //cout<<ii<<" "<<atime<<" tmin "<<mtime[ii]<<" tmax "<<mtime[ii+1]<<endl;
179     if ((atime>mtime[ii]) && (atime<=mtime[ii+1])) { // time interval found
180     ical_2nd = ii;
181     tmin_2nd = mtime[ii];
182     tmax_2nd = mtime[ii+1];
183     //cout<<"Time Interval 2nd order corr "<<ii<<" tmin "<<mtime[ii]<<" tmax "<<mtime[ii+1]<<endl;
184     break;
185     }
186     };
187     cout<<"2nd order corr: abs time "<<atime<<" interval "<<ical_2nd<<" tmin "<<mtime[ical_2nd]<<" tmax "<<mtime[ical_2nd+1]<<endl;
188     }
189     //-----------------------------------------------------------------------
190    
191     //== interpolate betwen time limits
192    
193     thelp1 = mtime[ical_2nd];
194     thelp2 = mtime[ical_2nd+1];
195    
196     TArrayF &corr1 = dedx_corr_mp[ical_2nd];
197     TArrayF &corr2 = dedx_corr_mp[ical_2nd+1];
198    
199     TArrayF &corr3 = dedx_corr_mhe[ical_2nd];
200     TArrayF &corr4 = dedx_corr_mhe[ical_2nd+1];
201    
202     TArrayF &corr5 = dedx_corr_mc[ical_2nd];
203     TArrayF &corr6 = dedx_corr_mc[ical_2nd+1];
204    
205     for (Int_t ii=0; ii<48;ii++) {
206     Float_t yhelp1 = corr1[ii];
207     Float_t yhelp2 = corr2[ii];
208     Float_t slope = (yhelp2-yhelp1)/(thelp2-thelp1);
209     Float_t inter = yhelp1 - slope*thelp1;
210     dedx_corrp[ii] = slope*atime + inter;
211     // if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
212     // cout<<ii<<" "<<yhelp1<<" "<<yhelp2<<" => "<<dedx_corr[ii]<<endl;
213    
214     yhelp1 = corr3[ii];
215     yhelp2 = corr4[ii];
216     slope = (yhelp2-yhelp1)/(thelp2-thelp1);
217     inter = yhelp1 - slope*thelp1;
218     dedx_corrhe[ii] = slope*atime + inter;
219    
220     yhelp1 = corr5[ii];
221     yhelp2 = corr6[ii];
222     slope = (yhelp2-yhelp1)/(thelp2-thelp1);
223     inter = yhelp1 - slope*thelp1;
224     dedx_corrc[ii] = slope*atime + inter;
225    
226     }
227    
228    
229     //--------------- TABLE OF PERIODS WITH HV PROBLEMS --------------------------
230    
231     int Aconn=0; // PMT 0,20,22,24
232     int Bconn=0; // PMT 6,12,26,34
233     int Cconn=0; // PMT 4,14,28,32
234     int Dconn=0; // PMT 2,8,10,30
235     int Econn=0; // PMT 42,43,44,47
236     int Fconn=0; // PMT 7,19,23,27
237     int Gconn=0; // PMT 3,11,25,33
238     int Hconn=0; // PMT 1,9,13,21
239     int Iconn=0; // PMT 5,29,31,35
240     int Lconn=0; // PMT 37,40,45,46
241     int Mconn=0; // PMT 15,16,17,18
242     int Nconn=0; // PMT 36,38,39,41
243    
244     int standard=0;
245    
246     int S115B_ok=0;
247     int S115B_break=0;
248    
249     if(atime>=1153660001 && atime<=1154375000)Dconn=1;
250     else if(atime>=1155850001 && atime<=1156280000){
251     Hconn=1;
252     Nconn=1;
253     }
254    
255     else if(atime>=1158160000 && atime<=1158220000)Cconn=1; // new period found WM 11-jan-2018
256    
257     else if(atime>=1168490001 && atime<=1168940000)Dconn=1;
258     else if(atime>=1168940001 && atime<=1169580000){
259     Fconn=1;
260     Mconn=1;
261     }
262    
263     else if(atime>=1174665001 && atime<=1175000000)Bconn=1;
264     else if(atime>=1176120001 && atime<=1176800000)Hconn=1;
265     else if(atime>=1176800001 && atime<=1178330000)Econn=1;
266     else if(atime>=1178330001 && atime<=1181322000)Hconn=1;
267     else if(atime>=1182100001 && atime<=1183030000)Aconn=1;
268     else if(atime>=1184000001 && atime<=1184570000)Hconn=1;
269     else if(atime>=1185090001 && atime<=1185212000)Dconn=1;
270     else if(atime>=1191100001 && atime<=1191940000)Dconn=1;
271     else if(atime>=1196230001 && atime<=1196280000)Hconn=1;
272     else if(atime>=1206100001 && atime<=1206375600)Cconn=1;
273     else if(atime>=1217989201 && atime<=1218547800)Econn=1;
274     else if(atime>=1225789201 && atime<=1226566800)Econn=1;
275     else if(atime>=1229400901 && atime<=1229700000)Econn=1;
276     else if(atime>=1230318001 && atime<=1230415200)Econn=1;
277    
278     else if(atime>=1320710401 && atime<=1322006400)Dconn=1; // (0,"","",1320710401,1322006400,213," D connector ",NULL);
279     else if(atime>=1309910401 && atime<=1310083200)Hconn=1; // (0,"","",1309910401,1310083200,217," H connector ",NULL);
280     else if(atime>=1384700001 && atime<=1387400000)Hconn=1; // (0,"","",1384700001,1387400000,217," H connector ",NULL);
281     else if(atime>=1230318001 && atime<=1230410000)Hconn=1; //wm period 38 HV_H 1230318001 - 1230410000 12-26 - 12-27
282    
283     else if(atime>=1426620001 && atime<=1427745000)Dconn=1; // 15-03-17 - 15-03-30 HV_D problems
284     else if(atime>=1433820001 && atime<=1433960000)Dconn=1; // 15-06-09 - 15-06-10 HV_D problems
285     else if(atime>=1434460001 && atime<=1435310000)Dconn=1; // 15-06-17 - 15-06-26 HV_D problems
286    
287    
288     else {
289     standard=1;
290     }
291    
292     if(atime<1158720000)S115B_ok=1;
293     else S115B_break=1;
294    
295    
296     //================================================================
297    
298     // cout<<" L2->GetNTracks(trkAlg) "<<L2->GetNTracks(trkAlg)<<endl;
299    
300     if( L2->GetNTracks(trkAlg)>=1 ){
301    
302     // PamTrack *trackTRK = L2->GetTrack(itr); // 9th
303     PamTrack *trackTRK = L2->GetTrack(0,trkAlg); // 10th
304    
305     Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]);
306    
307     Float_t def =trackTRK->GetExtTrack()->al[4];
308     Float_t chi =trackTRK->GetExtTrack()->chi2;
309     Float_t rig = 1./def;
310    
311     if(betamean<0.05 || betamean>2){
312     for(int i=0;i<48;i++)INFOpmt[i]=1;
313     for(int i=0;i<6;i++)INFOlayer[i]=1;
314     }
315    
316     // define angle:
317    
318     double dx = trackTRK->GetToFTrack()->xtr_tof[0] - trackTRK->GetToFTrack()->xtr_tof[5];
319     double dy = trackTRK->GetToFTrack()->ytr_tof[0] - trackTRK->GetToFTrack()->ytr_tof[5];
320     double dr = sqrt(dx*dx+dy*dy);
321     double theta=atan(dr/77.31);
322    
323     // TArrayF adc;
324     Float_t adc[48];
325     for(int i=0;i<48;i++){
326     adc[i]=0;
327     }
328    
329    
330     // Clear the dEdx arrays
331     for(int i=0;i<48;i++){
332     eDEDXpmt[i] = 1000;
333     }
334    
335     for(int i=0;i<24;i++) eDEDXpad[i] = 1000.;
336     for(int i=0;i<6;i++) eDEDXlayer[i] = 1000.;
337    
338     // Raw ADC in the ToFPMT class
339    
340     for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
341     ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
342     Int_t pmt_id = tofpmt->pmt_id;
343     adc[pmt_id] = tofpmt->adc ;
344     //cout<<ipmt<<" "<<pmt_id<<" "<<adc[pmt_id]<<endl;
345     } //loop on 48 pmts
346    
347     Int_t iverbose=0;
348    
349    
350     for (Int_t jj=0; jj<6; jj++){
351     xv[jj]=trackTRK->GetToFTrack()->xtr_tof[jj];
352     yv[jj]=trackTRK->GetToFTrack()->ytr_tof[jj];
353     //zin[jj] = L2->GetToFLevel2()->GetZTOF(L2->GetToFLevel2()->GetToFPlaneID(jj));
354     //cout<<"xv "<<xv[jj]<<" yv "<<yv[jj]<<endl;
355     }
356    
357    
358     Int_t paddleidoftrack[6]= {-1} ;
359     for (Int_t ilay=0; ilay<6; ilay++) {
360     paddleidoftrack[ilay] = L2->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
361     //cout<<ilay<<" "<<paddleidoftrack[ilay]<<endl;
362     }
363    
364     //-----------------------------------------------------------------------------------
365    
366     // Old loop ToFTrack class
367     // for (Int_t ipmt=0; ipmt<trackTRK->GetToFTrack()->npmtadc; ipmt++){
368     // Int_t ii = trackTRK->GetToFTrack()->pmtadc[ipmt];
369    
370     // New loop over PMTs found with GetPaddleIdOfTrack
371     for (Int_t ii=0; ii<48; ii++) {
372    
373     Int_t ilay;
374     if(ii<16) ilay=0;
375     if((15<ii)&&(ii<28)) ilay=1;
376     if((27<ii)&&(ii<32)) ilay=2;
377     if((31<ii)&&(ii<36)) ilay=3;
378     if((35<ii)&&(ii<42)) ilay=4;
379     if((41<ii)&&(ii<48)) ilay=5;
380    
381     if (paddleidoftrack[ilay] > -1) { // paddle hitted
382    
383     Float_t xpos=0;
384     if(ii<16) xpos = yv[0];
385     if((15<ii)&&(ii<28)) xpos = xv[1];
386     if((27<ii)&&(ii<32)) xpos = xv[2];
387     if((31<ii)&&(ii<36)) xpos = yv[3];
388     if((35<ii)&&(ii<42)) xpos = yv[4];
389     if((41<ii)&&(ii<48)) xpos = xv[5];
390    
391    
392     if(adc[ii]==4095)INFOpmt[ii]=2;
393     else if(adc[ii]>=PMTsat[ii]-5)INFOpmt[ii]=3;
394    
395    
396     if( adc[ii] >= PMTsat[ii]-5 ) continue;
397     if( adc[ii] <= 0. ) continue;
398    
399     double adcpC = f_adcPC( adc[ii] ); // - adc conversion in pC
400     double adccorr = adcpC*fabs(cos(theta));
401    
402     if(adccorr<=0)INFOpmt[ii]=4;
403     if(adccorr<=0.) continue;
404    
405    
406     double adcHe, adcnorm, adclin, dEdx, Zeta;
407    
408     adcHe=-2;
409     adcnorm=-2;
410     adclin=-2;
411     dEdx=-2;
412     Zeta=-2;
413    
414    
415     float dEdxH=-2;
416     double dEdxHe=-2;
417    
418     double day = (atime-1150000000)/84600;
419    
420    
421     Double_t correction = 1.;
422    
423     if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
424     correction = 1.675;
425     }
426     else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
427     correction = 2.482;
428     }
429     else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
430     correction = 1.464;
431     }
432     else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
433     correction = 1.995;
434     }
435     else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
436     correction = 1.273;
437     }
438     else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
439     correction = 1.565;
440     }
441     else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
442     correction = 1.565;
443     }
444     // else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
445     // correction = 1.018;
446     // }
447     //-- new WM 12-jan-2018
448     else if(Nconn==1 && (ii==36 || ii==38 || ii==39)){
449     correction = 1.34;
450     }
451     else if(Nconn==1 && (ii==41)){
452     correction = 1.0;
453     }
454     //-- end new WM
455    
456     else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
457     correction = 1.84;
458     }
459     else if(S115B_break==1 && ii==9 && Hconn==1){
460     correction = 1.64;
461     }
462     else if(S115B_break==1 && ii==9 && Hconn==0){
463     correction = 1.0;
464     }
465     else correction = 1.;
466    
467     // adcHe is the value of the pol5 as a function of the incident position, function derived for helium
468    
469     if( ii==9 && S115B_break==1 ){
470     adcHe = (f_att5B( xpos ))/correction;
471     } else {
472     adcHe = Get_adc_he(ii, xpos) / correction;
473     };
474    
475     if(adcHe<=0) continue;
476    
477     // adcnorm is f_pos(adccorr), where the f_pos function describes the nonlinearity of the PMT
478    
479     if(ii==9 && S115B_break==1) adcnorm = f_pos5B(adccorr);
480     else adcnorm = f_pos( (parPos[ii]), adccorr);
481    
482     if(adcnorm<=0) continue;
483    
484     // adclin is the linear calculation of the dEdx using adcnorm and the attenuation function for He
485    
486     if(ii==9 && S115B_break==1) adclin = 36.*adcnorm/adcHe;
487     else adclin = 4.*adcnorm/adcHe;
488    
489     if(adclin<=0) INFOpmt[ii]=5;
490     if(adclin<=0) continue;
491    
492     // dEdx = f_desatBB(adclin) uses a nonlinear function which takes Birks' effect etc. into account
493    
494     if(ii==9 && S115B_break==1) {
495     dEdx = f_desatBB5B( adclin );
496     }
497     else
498     {
499     dEdx = f_desatBB((parDesatBB[ii]), adclin );
500     // new linear behaviour for some PMTs...
501     if ((ii==8) || (ii==9) || (ii==14) || (ii==17) || (ii==28) || (ii==29) || (ii==31) || (ii==32) || (ii==33) || (ii==35) || (ii==40) || (ii==45)) {
502     if (adclin<10.) {
503     Float_t dEdx_help = f_desatBB((parDesatBB[ii]), 10. );
504     Float_t mhelp = dEdx_help / 10.;
505     dEdx = mhelp * adclin ;
506     }
507     }
508     } // else
509    
510     if(dEdx<0) INFOpmt[ii]=6;
511     if(dEdx<=0) continue;
512    
513     //--------------------- 2nd order correction ----------------------------------
514    
515     //------------------------ bi-scale -----------------------------------------
516     if ( strcmp(tri_or_bi,"BI") == 0){
517    
518     Float_t slope,inter;
519    
520     if(dedx_corrhe[ii]>dedx_corrp[ii]) slope=3./(dedx_corrhe[ii]-dedx_corrp[ii]);
521     else slope=1.;
522     if(dedx_corrhe[ii]>dedx_corrp[ii]) inter=1. - (slope*dedx_corrp[ii]);
523     else inter=0.;
524    
525     // For S115B_break dedx_corrp is NOT proton, but calibrated with carbon !!!
526     // dedx_corrhe is calibrated with helium, but probably not very reliable
527     // Do simple linear scaling:
528     if(ii==9 && S115B_break==1) {
529     if(dedx_corrp[ii]>10) slope=36./dedx_corrp[ii];
530     else slope=1.;
531     inter = 0.;
532     }
533    
534     dEdxH = inter + slope*dEdx;
535    
536     } // bi-scale
537    
538     //------------------------ tri-scale -----------------------------------------
539    
540     if ( strcmp(tri_or_bi,"TRI") == 0){
541    
542     Float_t Xa[3],Ya[3];
543     Ya[0] = 1;
544     Ya[1] = 4;
545     Ya[2] = 36;
546     Xa[0] = dedx_corrp[ii];
547     Xa[1] = dedx_corrhe[ii];
548     Xa[2] = dedx_corrc[ii];
549     TGraph *g1 = new TGraph(3,Xa,Ya);
550     g1->Fit("ftri","q");
551     dEdxH = ftri->Eval(dEdx);
552     g1->Delete();
553     }
554    
555     //---------------------------------------------------------------------------------
556    
557     // if (iverbose==1) cout<<"dEdx "<<dEdx<<" dEdxH "<<dEdxH<<endl;
558    
559     if(dEdxH!=-2)eDEDXpmt[ii]=(Float_t)dEdxH;
560     else eDEDXpmt[ii]=(Float_t)dEdx;
561    
562     // The f_BB=f(beta) correction is only used to calculate the charge Zeta
563    
564     if(dEdxH!=-2)dEdx=(Float_t)dEdxH;
565    
566     dEdxHe=-2;
567    
568     if(betamean<100) {
569    
570     if(ii==9 && S115B_break==1){
571     if( betamean <1. ) dEdxHe = f_BB5Bpol4( betamean );
572     else dEdxHe = f_BB5Bpol4( 1.0 );
573     } else {
574     if( betamean <1. ) dEdxHe = f_BBpol4( (parBBpol4[ii]), betamean );
575     else dEdxHe = f_BBpol4( (parBBpol4[ii]), 1.0 );
576     }
577    
578     if(ii==9 && S115B_break==1) Zeta = sqrt(36.*(dEdx/dEdxHe));
579     else Zeta = sqrt(4.*(dEdx/dEdxHe));
580    
581     if(Zeta<0) INFOpmt[ii]=7;
582    
583     eZpmt[ii]=(Float_t)Zeta;
584    
585     } // betamean < 100
586    
587    
588    
589    
590     // if (ii==9) printf("%5d %5.0f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %5.4f \n", ii, adc[ii], adcpC, adccorr, adcHe, dEdxHe, dEdx, Zeta, betamean );
591    
592     // if (ii==9) cout<<"eDEDXpmt[ii] "<<eDEDXpmt[ii]<<endl; //" eZpmt[ii] "<<eZpmt[ii]<<endl;
593    
594     } // paddleidoftrack > -1
595     } //end loop on PMTs 0-47
596    
597    
598     //-------------------------- paddle + layer ----------------------
599     //-----------------------------------------------------------------
600     // The GetdEdx method "handmade"
601     //-----------------------------------------------------------------
602    
603     Int_t ipad_a[6] = {0,8,14,16,18,21};
604     Int_t ipmt_a[6] = {0,16,28,32,36,42};
605    
606     for (Int_t ilay=0; ilay<6; ilay++) {
607     paddleidoftrack[ilay] = L2->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
608     }
609    
610    
611     for (Int_t ilay=0; ilay<6; ilay++) {
612     Float_t xhelp=1000.;
613     Int_t i1 = paddleidoftrack[ilay];
614     if (i1 != -1) {
615     Int_t ipad = ipad_a[ilay]+i1;
616     Int_t ihelp = ipmt_a[ilay]+(i1*2);
617     if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]==1000.)) xhelp = eDEDXpmt[ihelp];
618     if ((eDEDXpmt[ihelp+1]<1000.) && (eDEDXpmt[ihelp]==1000.)) xhelp = eDEDXpmt[ihelp+1];
619     if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmt[ihelp]+eDEDXpmt[ihelp+1]);
620     eDEDXpad[ipad] = xhelp;
621     eDEDXlayer[ilay] = xhelp;
622    
623     if ((eZpmt[ihelp]>-1) && (eZpmt[ihelp+1]==-1)) xhelp = eZpmt[ihelp];
624     if ((eZpmt[ihelp+1]>-1) && (eZpmt[ihelp]==-1)) xhelp = eZpmt[ihelp+1];
625     if ((eZpmt[ihelp]>-1) && (eZpmt[ihelp+1]>-1)) xhelp = 0.5*(eZpmt[ihelp]+eZpmt[ihelp+1]);
626     eZpad[ipad] = xhelp;
627     eZlayer[ilay] = xhelp;
628    
629    
630     }
631     } // ilay
632    
633    
634     } // if( L2->GetNTracks(trkAlg)>=1 )...
635    
636    
637    
638     //================================================================
639     //================================================================
640     //============== standalone dEdx and beta =============
641     //================================================================
642     //================================================================
643    
644     if (L2->GetToFStoredTrack(-1) ) {
645    
646     // ToF track
647    
648     ToFTrkVar *tof = L2->GetToFStoredTrack(-1);
649    
650     Float_t betamean = tof->CalcBeta(10., 10., 20.);
651    
652     if(betamean<0.05 || betamean>2){
653     for(int i=0;i<48;i++)INFOpmtstd[i]=1;
654     for(int i=0;i<6;i++)INFOlayerstd[i]=1;
655     }
656    
657     // define angle:
658    
659     double theta =0.;
660     if (fabs((tof->xtofpos[0])<100) && (fabs(tof->xtofpos[2])<100) && (fabs(tof->ytofpos[0])<100) && (fabs(tof->ytofpos[2])<100)) {
661     double dx = tof->xtofpos[0] - tof->xtofpos[2];
662     double dy = tof->ytofpos[0] - tof->ytofpos[2];
663     double dr = sqrt(dx*dx+dy*dy);
664     theta=atan(dr/77.31);
665     } // if...
666    
667     // Clear ADC;
668     Float_t adc[48];
669     for(int i=0;i<48;i++){
670     // adc[i]=0;
671     adc[i]=4095;
672     }
673    
674    
675     for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
676     ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
677     Int_t pmt_id = tofpmt->pmt_id;
678     adc[pmt_id] = tofpmt->adc ; // raw ADC
679     } //loop on 48 pmts
680    
681    
682     //====================================================================
683     //=========== Check ToF standalone using HitPaddle ==================
684     //====================================================================
685    
686     Int_t paddleidoftrack[6] = {-1, -1, -1, -1, -1, -1};;
687    
688     for(Int_t jj=0; jj<8; jj++){
689     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(0,jj);
690     if (HitPad==1) paddleidoftrack[0] = jj;
691     }
692     for(Int_t jj=0; jj<6; jj++){
693     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(1,jj);
694     if (HitPad==1) paddleidoftrack[1] = jj;
695     }
696    
697     for(Int_t jj=0; jj<2; jj++){
698     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(2,jj);
699     if (HitPad==1) paddleidoftrack[2] = jj;
700     }
701     for(Int_t jj=0; jj<2; jj++){
702     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(3,jj);
703     if (HitPad==1) paddleidoftrack[3] = jj;
704     }
705    
706     for(Int_t jj=0; jj<3; jj++){
707     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(4,jj);
708     if (HitPad==1) paddleidoftrack[4] = jj;
709     }
710     for(Int_t jj=0; jj<3; jj++){
711     Int_t HitPad = L2->GetToFLevel2()->HitPaddle(5,jj);
712     if (HitPad==1) paddleidoftrack[5] = jj;
713     }
714    
715     /*
716     cout<<"paddleidoftrack std ";
717     for (Int_t ilay=0; ilay<6; ilay++) {
718     cout<<" "<<paddleidoftrack[ilay];
719     }
720     cout<<endl;
721     */
722     //==================================================================
723     //==================================================================
724     //==================================================================
725    
726     // Clear the dEdx arrays
727     for(int i=0;i<48;i++){
728     eDEDXpmtstd[i] = 1000;
729     }
730    
731     for(int i=0;i<24;i++) eDEDXpadstd[i] = 1000.;
732     for(int i=0;i<6;i++) eDEDXlayerstd[i] = 1000.;
733    
734    
735     Int_t ipad_a[6] = {0,8,14,16,18,21};
736     Int_t ipmt_a[6] = {0,16,28,32,36,42};
737    
738     Int_t pmthit[48];
739     for(int i=0;i<48;i++) pmthit[i] = 0;
740    
741     for(int i=0;i<6;i++) {
742     if (paddleidoftrack[i] > -1) {
743     Int_t ipad = ipad_a[i] + paddleidoftrack[i];
744     pmthit[2*ipad] = 1;
745     pmthit[2*ipad+1] = 1;
746     }
747     }
748    
749    
750    
751     for (Int_t ipmt=0; ipmt<48; ipmt++){
752     if (pmthit[ipmt] == 1) {
753    
754     Int_t ii = ipmt;
755    
756     Float_t xpos=0;
757     if(ii<16) xpos = tof->ytofpos[0];
758     if((15<ii)&&(ii<28)) xpos = tof->xtofpos[0];
759     if((27<ii)&&(ii<32)) xpos = tof->xtofpos[1];
760     if((31<ii)&&(ii<36)) xpos = tof->ytofpos[1];
761     if((35<ii)&&(ii<42)) xpos = tof->ytofpos[2];
762     if((41<ii)&&(ii<48)) xpos = tof->xtofpos[2];
763    
764    
765     if(adc[ii]==4095)INFOpmtstd[ii]=2;
766     else if(adc[ii]>=PMTsat[ii]-5)INFOpmtstd[ii]=3;
767    
768     //cout<<ii<<" "<<adc[ii]<<" "<<PMTsat[ii]<<endl;
769    
770     if( adc[ii] >= PMTsat[ii]-5 ) continue;
771     if( adc[ii] <= 0. ) continue;
772    
773     double adcpC = f_adcPC( adc[ii] ); // - adc conversion in pC
774     double adccorr = adcpC*fabs(cos(theta));
775    
776     //cout<<ii<<" adc "<<adc[ii]<<" adcpC "<<adcpC<<" adccorr "<<adccorr<<endl;
777    
778    
779     if(adccorr<=0)INFOpmt[ii]=4;
780     if(adccorr<=0.) continue;
781    
782    
783     double adcHe, adcnorm, adclin, dEdx, Zeta;
784    
785     adcHe=-2;
786     adcnorm=-2;
787     adclin=-2;
788     dEdx=-2;
789     Zeta=-2;
790    
791     float dEdxH=-2;
792     double dEdxHe=-2;
793    
794     // double day = (atime-1150000000)/84600;
795    
796    
797     Double_t correction = 1.;
798    
799     if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
800     correction = 1.675;
801     }
802     else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
803     correction = 2.482;
804     }
805     else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
806     correction = 1.464;
807     }
808     else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
809     correction = 1.995;
810     }
811     else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
812     correction = 1.273;
813     }
814     else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
815     correction = 1.565;
816     }
817     else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
818     correction = 1.565;
819     }
820     // else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
821     // correction = 1.018;
822     // }
823     //-- new WM 12-jan-2018
824     else if(Nconn==1 && (ii==36 || ii==38 || ii==39)){
825     correction = 1.34;
826     }
827     else if(Nconn==1 && (ii==41)){
828     correction = 1.0;
829     }
830     //-- end new WM
831    
832     else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
833     correction = 1.84;
834     }
835     else if(S115B_break==1 && ii==9 && Hconn==1){
836     correction = 1.64;
837     }
838     else if(S115B_break==1 && ii==9 && Hconn==0){
839     correction = 1.0;
840     }
841     else correction = 1.;
842    
843     // adcHe is the value of the pol5 as a function of the incident position, curve derived for helium
844    
845     if( ii==9 && S115B_break==1 ){
846     adcHe = (f_att5B( xpos ))/correction;
847     } else {
848     adcHe = Get_adc_he(ii, xpos) / correction;
849     };
850    
851     if(adcHe<=0) continue;
852    
853     // adcnorm is f_pos(adccorr), where the f_pos function describes the nonlinearity of the PMT
854    
855     if(ii==9 && S115B_break==1) adcnorm = f_pos5B(adccorr);
856     else adcnorm = f_pos( (parPos[ii]), adccorr);
857    
858     //cout<<"adcnorm "<<adcnorm<<endl;
859    
860     if(adcnorm<=0) continue;
861    
862     // adclin is the linear calculation of the dEdx using adcnorm and the attenuation function for He
863    
864     if(ii==9 && S115B_break==1) adclin = 36.*adcnorm/adcHe;
865     else adclin = 4.*adcnorm/adcHe;
866    
867     //cout<<"adclin "<<adclin<<endl;
868    
869     if(adclin<=0) INFOpmt[ii]=5;
870     if(adclin<=0) continue;
871    
872     // dEdx = f_desatBB(adclin) uses a nonlinear function which takes Birks' effect etc. into account
873    
874    
875     if(ii==9 && S115B_break==1) {
876     dEdx = f_desatBB5B( adclin );
877     }
878     else
879     {
880     dEdx = f_desatBB((parDesatBB[ii]), adclin );
881     // new linear behaviour for some PMTs...
882     if ((ii==8) || (ii==9) || (ii==14) || (ii==17) || (ii==28) || (ii==29) || (ii==31) || (ii==32) || (ii==33) || (ii==35) || (ii==40) || (ii==45)) {
883     if (adclin<10.) {
884     Float_t dEdx_help = f_desatBB((parDesatBB[ii]), 10. );
885     Float_t mhelp = dEdx_help / 10.;
886     dEdx = mhelp * adclin ;
887     }
888     // if (iverbose==1) cout<<"dEdx after linear correction f_desatBB "<<dEdx<<endl;
889     }
890     } // else
891    
892     if(dEdx<0) INFOpmtstd[ii]=6;
893     if(dEdx<=0) continue;
894    
895    
896     //--------------------- 2nd order correction ----------------------------------
897    
898     //------------------------ bi-scale -----------------------------------------
899     if ( strcmp(tri_or_bi,"BI") == 0){
900    
901     Float_t slope,inter;
902    
903     if(dedx_corrhe[ii]>dedx_corrp[ii]) slope=3./(dedx_corrhe[ii]-dedx_corrp[ii]);
904     else slope=1.;
905     if(dedx_corrhe[ii]>dedx_corrp[ii]) inter=1. - (slope*dedx_corrp[ii]);
906     else inter=0.;
907    
908     // For S115B_break dedx_corrp is NOT proton, but calibrated with carbon !!!
909     // dedx_corrhe is calibrated with helium, but probably not very reliable
910     // Do simple linear scaling:
911     if(ii==9 && S115B_break==1) {
912     if(dedx_corrp[ii]>10) slope=36./dedx_corrp[ii];
913     else slope=1.;
914     inter = 0.;
915     }
916    
917     dEdxH = inter + slope*dEdx;
918    
919     } // bi-scale
920    
921     //------------------------ tri-scale -----------------------------------------
922    
923     if ( strcmp(tri_or_bi,"TRI") == 0){
924    
925     Float_t Xa[3],Ya[3];
926     Ya[0] = 1;
927     Ya[1] = 4;
928     Ya[2] = 36;
929     Xa[0] = dedx_corrp[ii];
930     Xa[1] = dedx_corrhe[ii];
931     Xa[2] = dedx_corrc[ii];
932     TGraph *g1 = new TGraph(3,Xa,Ya);
933     g1->Fit("ftri","q");
934     dEdxH = ftri->Eval(dEdx);
935     g1->Delete();
936     }
937    
938     //---------------------------------------------------------------------------------
939    
940     // if (iverbose==1) cout<<"dEdx "<<dEdx<<" dEdxH "<<dEdxH<<endl;
941    
942     if(dEdxH!=-2)eDEDXpmtstd[ii]=(Float_t)dEdxH;
943     else eDEDXpmtstd[ii]=(Float_t)dEdx;
944    
945    
946     // The f_BB=f(beta) correction is only used in ToFNaNuclei to calculate the charge Zeta
947    
948     dEdxHe=-2;
949    
950     if(betamean<100) {
951    
952     if(ii==9 && S115B_break==1){
953     if( betamean <1. ) dEdxHe = f_BB5Bpol4( betamean );
954     else dEdxHe = f_BB5Bpol4( 1.0 );
955     } else {
956     if( betamean <1. ) dEdxHe = f_BBpol4( (parBBpol4[ii]), betamean );
957     else dEdxHe = f_BBpol4( (parBBpol4[ii]), 1.0 );
958     }
959    
960     if(ii==9 && S115B_break==1) Zeta = sqrt(36.*(dEdx/dEdxHe));
961     else Zeta = sqrt(4.*(dEdx/dEdxHe));
962    
963     if(Zeta<0) INFOpmtstd[ii]=7;
964    
965     eZpmtstd[ii]=(Float_t)Zeta;
966    
967     } // betamean < 100
968    
969    
970    
971     // printf("%5d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %5.4f \n", ii, adcpC, adccorr, adcHe, dEdxHe, dEdx, Zeta, betamean );
972    
973     //cout<<"eDEDXpmtstd[ii] "<<eDEDXpmt[ii]<<" eZpmtstd[ii] "<<eZpmtstd[ii]<<endl;
974    
975     } // hitted PMTs
976     } //end loop on PMTs along the track
977    
978     //-----------------------------------------------------------------
979     // The GetdEdx method "handmade"
980     //-----------------------------------------------------------------
981    
982    
983     for (Int_t ilay=0; ilay<6; ilay++) {
984     Float_t xhelp=1000.;
985     Int_t i1 = paddleidoftrack[ilay];
986     if (i1 != -1) {
987     Int_t ipad = ipad_a[ilay]+i1;
988     Int_t ihelp = ipmt_a[ilay]+(i1*2);
989     if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]==1000.)) xhelp = eDEDXpmtstd[ihelp];
990     if ((eDEDXpmtstd[ihelp+1]<1000.) && (eDEDXpmtstd[ihelp]==1000.)) xhelp = eDEDXpmtstd[ihelp+1];
991     if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmtstd[ihelp]+eDEDXpmtstd[ihelp+1]);
992     eDEDXpadstd[ipad] = xhelp;
993     eDEDXlayerstd[ilay] = xhelp;
994    
995     if ((eZpmtstd[ihelp]>-1) && (eZpmtstd[ihelp+1]==-1)) xhelp = eZpmtstd[ihelp];
996     if ((eZpmtstd[ihelp+1]>-1) && (eZpmtstd[ihelp]==-1)) xhelp = eZpmtstd[ihelp+1];
997     if ((eZpmtstd[ihelp]>-1) && (eZpmtstd[ihelp+1]>-1)) xhelp = 0.5*(eZpmtstd[ihelp]+eZpmtstd[ihelp+1]);
998     eZpadstd[ipad] = xhelp;
999     eZlayerstd[ilay] = xhelp;
1000    
1001    
1002     }
1003     } // ilay
1004    
1005     } // if (pam_event->GetToFStoredTrack(-1) )...
1006    
1007     ftri->Delete();
1008    
1009     };
1010    
1011    
1012     //------------------------------------------------------------------------
1013     void ToFdEdx_patch::PrintTD()
1014     {
1015     for(int i=0; i<48; i++) {
1016     /*
1017     TArrayF &binx = TDx[i];
1018     TArrayF &biny = TDy[i];
1019     for(int k=0; k<200; k++) { // bin temporali
1020     printf("%d %d %f %f", i,k, binx[k], biny[k]);
1021    
1022     }
1023     */
1024     }
1025     }
1026    
1027    
1028     //------------------------------------------------------------------------
1029     void ToFdEdx_patch::Define_PMTsat()
1030     {
1031     Float_t sat[48] = {
1032     3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
1033     3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
1034     3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
1035     3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
1036     3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
1037     3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
1038     PMTsat.Set(48,sat);
1039     }
1040    
1041     //------------------------------------------------------------------------
1042     /*
1043     void ToFdEdx_patch::ReadParTD( Int_t ipmt, const char *fname )
1044     {
1045     printf("read %s\n",fname);
1046     if(ipmt<0) return;
1047     if(ipmt>47) return;
1048     FILE *fattin = fopen( fname , "r" );
1049     Float_t yTD[200],xTD[200];
1050     for(int j=0;j<200;j++){
1051     float x,y,ym,e;
1052     if(fscanf(fattin,"%f %f %f %f",
1053     &x, &y, &ym, &e )!=4) break;
1054     xTD[j]=x;
1055     if(ym>0&&fabs(y-ym)>1) yTD[j]=ym;
1056     else yTD[j]=y;
1057     }
1058     TDx[ipmt].Set(200,xTD);
1059     TDy[ipmt].Set(200,yTD);
1060     fclose(fattin);
1061     }
1062     */
1063     //------------------------------------------------------------------------
1064     void ToFdEdx_patch::ReadParBBpos( const char *fname )
1065     {
1066     printf("read %s\n",fname);
1067     parBBpos.Set(48);
1068     FILE *fattin = fopen( fname , "r" );
1069     for (int i=0; i<48; i++) {
1070     int tid=0;
1071     float tp;
1072     if(fscanf(fattin,"%d %f",
1073     &tid, &tp )!=2) break;
1074     parBBpos[i]=tp;
1075     }
1076     fclose(fattin);
1077     }
1078    
1079     //------------------------------------------------------------------------
1080     void ToFdEdx_patch::ReadParDesatBB( const char *fname )
1081     {
1082     printf("read %s\n",fname);
1083     FILE *fattin = fopen( fname , "r" );
1084     for (int i=0; i<48; i++) {
1085     int tid=0;
1086     float tp[3];
1087     if(fscanf(fattin,"%d %f %f %f",
1088     &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1089     parDesatBB[i].Set(3,tp);
1090     }
1091     fclose(fattin);
1092     }
1093    
1094    
1095     //------------------------------------------------------------------------
1096     void ToFdEdx_patch::ReadParBBneg( const char *fname )
1097    
1098     {
1099     printf("read %s\n",fname);
1100     FILE *fattin = fopen( fname , "r" );
1101     for (int i=0; i<48; i++) {
1102     int tid=0;
1103     float tp[3];
1104     if(fscanf(fattin,"%d %f %f %f",
1105     &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
1106     parBBneg[i].Set(3,tp);
1107     }
1108     fclose(fattin);
1109     }
1110    
1111     //------------------------------------------------------------------------
1112     void ToFdEdx_patch::ReadParBBpol4( const char *fname )
1113    
1114     {
1115     printf("read %s\n",fname);
1116     FILE *fattin = fopen( fname , "r" );
1117     for (int i=0; i<48; i++) {
1118     int tid=0;
1119     float tp[5];
1120     if(fscanf(fattin,"%d %f %f %f %f %f",
1121     &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4] )!=6) break;
1122     parBBpol4[i].Set(5,tp);
1123     }
1124     fclose(fattin);
1125     }
1126     //------------------------------------------------------------------------
1127     void ToFdEdx_patch::ReadParPos( const char *fname )
1128     {
1129     printf("read %s\n",fname);
1130     FILE *fattin = fopen( fname , "r" );
1131     for (int i=0; i<48; i++) {
1132     int tid=0;
1133     float tp[4];
1134     if(fscanf(fattin,"%d %f %f %f %f",
1135     &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
1136     parPos[i].Set(4,tp);
1137     }
1138     fclose(fattin);
1139     }
1140    
1141     //------------------------------------------------------------------------
1142     /*void ToFdEdx_patch::ReadParAtt( const char *fname )
1143     {
1144     printf("read %s\n",fname);
1145     FILE *fattin = fopen( fname , "r" );
1146     for (int i=0; i<48; i++) {
1147     int tid=0;
1148     float tp[6];
1149     if(fscanf(fattin,"%d %f %f %f %f %f %f",
1150     &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
1151     parAtt[i].Set(6,tp);
1152     }
1153     fclose(fattin);
1154     }
1155     */
1156    
1157     //---------------------------------------------------------------------------
1158    
1159     double ToFdEdx_patch::f_att(float C0, float C1, float C2, float C3, float C4, float C5, float x)
1160     {
1161    
1162     //cout<<"in f_att: para "<<C0<<" "<<C1<<" "<<C2<<" "<<C3<<" "<<C4<<" "<<C5<<" x "<<x<<endl;
1163    
1164     return
1165     C0 +
1166     C1*x +
1167     C2*x*x +
1168     C3*x*x*x +
1169     C4*x*x*x*x +
1170     C5*x*x*x*x*x;
1171     }
1172    
1173    
1174     /*
1175     double ToFdEdx_patch::f_att( TArrayF &p, float x )
1176     {
1177    
1178     //cout<<"in f_att: para "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" "<<p[3]<<" "<<p[4]<<" "<<p[5]<<" x "<<x<<endl;
1179    
1180     return
1181     p[0] +
1182     p[1]*x +
1183     p[2]*x*x +
1184     p[3]*x*x*x +
1185     p[4]*x*x*x*x +
1186     p[5]*x*x*x*x*x;
1187     }
1188     */
1189    
1190     //------------------------------------------------------------------------
1191     double ToFdEdx_patch::f_att5B( float x )
1192     {
1193     return
1194     101.9409 +
1195     6.643781*x +
1196     0.2765518*x*x +
1197     0.004617647*x*x*x +
1198     0.0006195132*x*x*x*x +
1199     0.00002813734*x*x*x*x*x;
1200     }
1201    
1202    
1203     double ToFdEdx_patch::f_pos( TArrayF &p, float x )
1204     {
1205     return
1206     p[0] +
1207     p[1]*x +
1208     p[2]*x*x +
1209     p[3]*x*x*x;
1210     }
1211    
1212     double ToFdEdx_patch::f_pos5B( float x )
1213     {
1214     return
1215     15.45132 +
1216     0.8369721*x +
1217     0.0005*x*x;
1218     }
1219    
1220    
1221    
1222     double ToFdEdx_patch::f_adcPC( float x )
1223     {
1224     return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
1225     }
1226    
1227     /*
1228     float ToFdEdx_patch::Get_adc_he( int id, float pl_x[6], float pl_y[6])
1229     {
1230    
1231     //
1232     // input: id - pmt [0:47}
1233     // pl_x - coord x of the tof plane
1234     // pl_y - coord y
1235    
1236    
1237    
1238     adc_he = 0;
1239     // if( eGeom.GetXY(id)==1 ) adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
1240     // if( eGeom.GetXY(id)==2 ) adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
1241     if( eGeom.GetXY(id)==1 ) adc_he = f_att( A0[id],A1[id],A2[id],A3[id], A4[id], A5[id], pl_x[eGeom.GetPlane(id)] );
1242     if( eGeom.GetXY(id)==2 ) adc_he = f_att( A0[id],A1[id],A2[id],A3[id], A4[id], A5[id], pl_y[eGeom.GetPlane(id)] );
1243    
1244     //cout<<"in Get_adc_he "<<adc_he<<endl;
1245    
1246     return adc_he;
1247     }
1248     */
1249    
1250     //------------------------------------------------------------------------
1251    
1252     float ToFdEdx_patch::Get_adc_he( int id, float x)
1253     {
1254     // input: id - pmt [0:47] , position in the paddle
1255    
1256     adc_he = 0;
1257     adc_he = f_att( A0[id],A1[id],A2[id],A3[id], A4[id], A5[id], x );
1258    
1259     //cout<<"in Get_adc_he "<<adc_he<<endl;
1260    
1261     return adc_he;
1262     }
1263    
1264     //------------------------------------------------------------------------
1265     double ToFdEdx_patch::f_BB( TArrayF &p, float x )
1266     {
1267     return p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
1268     }
1269    
1270     //------------------------------------------------------------------------
1271     double ToFdEdx_patch::f_BBpol4( TArrayF &p, float x )
1272     {
1273     return p[0] + p[1]*x + p[2]*x*x + p[3]*x*x*x + p[4 ]*x*x*x*x;
1274     }
1275    
1276     //------------------------------------------------------------------------
1277     double ToFdEdx_patch::f_BB5B( float x )
1278     {
1279     return 0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
1280     }
1281     //------------------------------------------------------------------------
1282     double ToFdEdx_patch::f_BB5Bpol4( float x )
1283     {
1284     // return 138.165 - 120.884*x - 122.502*x*x + 203.795*x*x*x - 64.1068*x*x*x*x ;
1285     return 67.2843 + 223.175*x - 721.068*x*x + 646.677*x*x*x - 181.868*x*x*x*x ;
1286    
1287     }
1288     //------------------------------------------------------------------------
1289     double ToFdEdx_patch::f_desatBB( TArrayF &p, float x )
1290     {
1291     return
1292     p[0] +
1293     p[1]*x +
1294     p[2]*x*x;
1295     }
1296    
1297     //------------------------------------------------------------------------
1298     double ToFdEdx_patch::f_desatBB5B( float x )
1299     {
1300     return
1301     -2.4 +
1302     0.75*x +
1303     0.009*x*x;
1304     }
1305    
1306     //------------------------------------------------------------------------
1307     void ToFdEdx_patch::ReadParAtt(const char *pardir, const char *param)
1308     {
1309     char filename[100],filename1[100],filename2[100],name[50];
1310     char str1[30];
1311    
1312    
1313     //--------------------- flight ---------------------------
1314    
1315     if (strcmp(param,"simu") != 0){
1316     cout<<"Flight calibrations "<<endl;
1317    
1318     sprintf(filename,"%stime_intervals_atten_10th.txt",pardir);
1319    
1320     cout<<"read Time Interval File for Attenuation : "<<filename<<endl;
1321    
1322     FILE *ftimein = fopen( filename , "r" );
1323    
1324     for (int ii=0; ii<100; ii++) {
1325     UInt_t tp[2];
1326     if(fscanf(ftimein,"%s %u %u",
1327     &str1[0], &tp[0], &tp[1])!=3) break;
1328     cout<<str1<<" "<<tp[0]<<" "<<tp[1]<<endl;
1329     T_int_min[ii] = tp[0];
1330     T_int_max[ii] = tp[1];
1331    
1332    
1333     sprintf(name,"%s",str1);
1334     // cout<<"test name "<<name<<endl;
1335     sprintf(filename1,"%s%s",pardir,name);
1336    
1337     cout<<"read Attenuation file "<<filename1<<endl;
1338     FILE *fattin = fopen( filename1, "r" );
1339    
1340     Float_t A0help[48];
1341     Float_t A1help[48];
1342     Float_t A2help[48];
1343     Float_t A3help[48];
1344     Float_t A4help[48];
1345     Float_t A5help[48];
1346    
1347    
1348     for (int i=0; i<48; i++) {
1349     int id=0;
1350     float par[6];
1351     if(fscanf(fattin,"%d %f %f %f %f %f %f",
1352     &id, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5] )!=7) break;
1353     A0help[i] = par[0];
1354     A1help[i] = par[1];
1355     A2help[i] = par[2];
1356     A3help[i] = par[3];
1357     A4help[i] = par[4];
1358     A5help[i] = par[5];
1359     }
1360    
1361     A0_array[ii].Set(48,A0help);
1362     A1_array[ii].Set(48,A1help);
1363     A2_array[ii].Set(48,A2help);
1364     A3_array[ii].Set(48,A3help);
1365     A4_array[ii].Set(48,A4help);
1366     A5_array[ii].Set(48,A5help);
1367    
1368    
1369     fclose(fattin);
1370    
1371    
1372     // if (T_int_max[ii] == 2000000000) break; // end of file
1373     if (T_int_max[ii] > 2000000000) break; // end of file
1374    
1375     }
1376     fclose(ftimein);
1377    
1378     } // flight
1379    
1380     //--------------------- simu ---------------------------
1381    
1382     if ( strcmp(param,"simu") == 0){
1383    
1384     //cout<<"Simulation calibrations "<<endl;
1385     sprintf(name,"tofdedx_gp_atten.dat");
1386     sprintf(filename1,"%s%s",pardir,name);
1387    
1388     cout<<"read Attenuation file "<<filename1<<endl;
1389     FILE *fattin = fopen( filename1, "r" );
1390    
1391     Float_t A0help[48];
1392     Float_t A1help[48];
1393     Float_t A2help[48];
1394     Float_t A3help[48];
1395     Float_t A4help[48];
1396     Float_t A5help[48];
1397    
1398    
1399     for (int i=0; i<48; i++) {
1400     int id=0;
1401     float par[6];
1402     if(fscanf(fattin,"%d %f %f %f %f %f %f",
1403     &id, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5] )!=7) break;
1404     A0help[i] = par[0];
1405     A1help[i] = par[1];
1406     A2help[i] = par[2];
1407     A3help[i] = par[3];
1408     A4help[i] = par[4];
1409     A5help[i] = par[5];
1410     }
1411    
1412     A0_array[0].Set(48,A0help);
1413     A1_array[0].Set(48,A1help);
1414     A2_array[0].Set(48,A2help);
1415     A3_array[0].Set(48,A3help);
1416     A4_array[0].Set(48,A4help);
1417     A5_array[0].Set(48,A5help);
1418    
1419    
1420     fclose(fattin);
1421    
1422     T_int_min[0] = 1000000000;
1423     T_int_max[0] = 2000000000;
1424    
1425    
1426     } // simu
1427    
1428    
1429    
1430    
1431     // set first time interval
1432     ical_atten = 0;
1433     tmin_atten = T_int_min[0];
1434     tmax_atten = T_int_max[0];
1435    
1436     // set first time interval attenuation
1437    
1438     TArrayF &Apar0 = A0_array[0];
1439     TArrayF &Apar1 = A1_array[0];
1440     TArrayF &Apar2 = A2_array[0];
1441     TArrayF &Apar3 = A3_array[0];
1442     TArrayF &Apar4 = A4_array[0];
1443     TArrayF &Apar5 = A5_array[0];
1444    
1445     for (Int_t ii=0; ii<48;ii++) {
1446     A0[ii]=Apar0[ii];
1447     A1[ii]=Apar1[ii];
1448     A2[ii]=Apar2[ii];
1449     A3[ii]=Apar3[ii];
1450     A4[ii]=Apar4[ii];
1451     A5[ii]=Apar5[ii];
1452     }
1453    
1454     cout<<"First time interval attenuation: "<<tmin_atten<<" - "<<tmax_atten<<endl;
1455    
1456    
1457     //-----------------------------------------------------------
1458     // Here I read the "bi-scale" dEdx_korr parameters 2nd order
1459     //-----------------------------------------------------------
1460    
1461     //Double_t t1,t2,tm;
1462     //UInt_t t1,t2,tm;
1463     Double_t t1,t2,tm;
1464    
1465     Float_t xm1,xm2,xm3;
1466     Float_t pmthelpp[48];
1467     Float_t pmthelphe[48];
1468     Float_t pmthelpc[48];
1469    
1470     //--------------------- flight ---------------------------
1471     if ( strcmp(param,"simu") != 0){
1472    
1473     sprintf(filename2,"%striscale_tofdedx_patch.txt",pardir);
1474    
1475     cout<<"Filename 2nd order tri-scale "<<filename2<<endl;
1476     ifstream fin(filename2);
1477    
1478     Int_t jj=0;
1479     Int_t idummy;
1480     while (! fin.eof()) {
1481     fin>>t1>>tm>>t2;
1482     //cout<<jj<<" "<<tm<<endl;
1483     cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
1484     //cout <<t1<<" "<<tm<<" "<<t2<<endl;
1485     mtime[jj]=tm;
1486     for (Int_t ii=0; ii<48;ii++) {
1487     fin>>xm1>>xm2>>xm3; // He-corr, P-corr, C-corr
1488     pmthelphe[ii] = xm1;
1489     pmthelpp[ii] = xm2;
1490     pmthelpc[ii] = xm3;
1491     }
1492     dedx_corr_mp[jj].Set(48,pmthelpp);
1493     dedx_corr_mhe[jj].Set(48,pmthelphe);
1494     dedx_corr_mc[jj].Set(48,pmthelpc);
1495     jj=jj+1;
1496     }
1497    
1498     /*
1499     Int_t idummy;
1500     for (Int_t jj=0; jj<2000; jj++) {
1501     fin>>t1>>tm>>t2;
1502     cout<<jj<<" "<<t1<<" "<<tm<<" "<<t2<<endl;
1503     if (t2==2000000000) break;
1504     cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
1505     mtime[jj]=tm;
1506     for (Int_t ii=0; ii<48;ii++) {
1507     fin>>idummy>>xmean1>>xwidth1;
1508     dedx_corr_m[jj][ii]=xmean1;
1509     }
1510     }
1511     */
1512    
1513     fin.close();
1514    
1515     } // flight
1516    
1517     //--------------------- simu ---------------------------
1518    
1519     if ( strcmp(param,"simu") == 0){
1520    
1521     sprintf(filename2,"%striscale_tofdedx_patch_simu.txt",pardir);
1522    
1523     cout<<"Filename 2nd order "<<filename2<<endl;
1524     ifstream fin(filename2);
1525    
1526     Int_t jj=0;
1527     Int_t idummy;
1528     while (! fin.eof()) {
1529     fin>>t1>>tm>>t2;
1530     //cout<<jj<<" "<<tm<<endl;
1531     //cout << setiosflags(ios::fixed) << setw(10) << setprecision(3) << tm << endl;
1532     cout <<t1<<" "<<tm<<" "<<t2<<endl;
1533     mtime[jj]=tm;
1534     for (Int_t ii=0; ii<48;ii++) {
1535     fin>>xm1>>xm2>>xm3; // He-corr, P-corr, C-corr
1536     pmthelphe[ii] = xm1;
1537     pmthelpp[ii] = xm2;
1538     pmthelpc[ii] = xm3;
1539     }
1540     dedx_corr_mp[jj].Set(48,pmthelpp);
1541     dedx_corr_mhe[jj].Set(48,pmthelphe);
1542     dedx_corr_mc[jj].Set(48,pmthelpc);
1543     jj=jj+1;
1544     }
1545    
1546    
1547     fin.close();
1548    
1549     } // simu
1550    
1551    
1552     // set first time interval
1553     ical_2nd=0;
1554     tmin_2nd = mtime[0];
1555     tmax_2nd = mtime[1];
1556    
1557     cout<<"First time interval 2nd order corr.: "<<tmin_2nd<<" - "<<tmax_2nd<<endl;
1558    
1559    
1560     }
1561    
1562    
1563     //----------------------------------------------------------------------------
1564    
1565    
1566    
1567    
1568    
1569    
1570    
1571    
1572    
1573    
1574    
1575    
1576    

  ViewVC Help
Powered by ViewVC 1.1.23