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

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

  ViewVC Help
Powered by ViewVC 1.1.23