#include <ToFdEdx_patch.h>
#include <TGraph.h>
#include <TSpline.h>
#include <TH1F.h>
#include <TMVA/TSpline1.h>

//------------------------------------------------------------------------
void ToFdEdx_patch::Clear(Option_t *option)
{
  //
  L2    = 0;
  OBT   = 0;
  PKT   = 0;
  atime = 0;
  tr    = 0;
  //sntr  = 0;
  // Set arrays and initialize structure

  eDEDXpmt.Set(48);    eDEDXpmt.Reset(-1);   // Set array size  and reset structure
  eZpmt.Set(48);       eZpmt.Reset(-1);
  eDEDXpad.Set(24);    eDEDXpad.Reset(-1);
  eZpad.Set(24);       eZpad.Reset(-1); 
  eDEDXlayer.Set(6);   eDEDXlayer.Reset(-1);
  eZlayer.Set(6);      eZlayer.Reset(-1);
  INFOpmt.Set(48);     INFOpmt.Reset(0);
  INFOlayer.Set(6);    INFOlayer.Reset(0);
  //
  
  // ToF standalone 
  eDEDXpmtstd.Set(48);    eDEDXpmtstd.Reset(-1);   // Set array size  and reset structure
  eZpmtstd.Set(48);       eZpmtstd.Reset(-1);
  eDEDXpadstd.Set(24);    eDEDXpadstd.Reset(-1);
  eZpadstd.Set(24);       eZpadstd.Reset(-1); 
  eDEDXlayerstd.Set(6);   eDEDXlayerstd.Reset(-1);
  eZlayerstd.Set(6);      eZlayerstd.Reset(-1);
  INFOpmtstd.Set(48);     INFOpmtstd.Reset(0);
  INFOlayerstd.Set(6);    INFOlayerstd.Reset(0);

};

//------------------------------------------------------------------------
void ToFdEdx_patch::Print(Option_t *option)
{
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
};


//------------------------------------------------------------------------
void ToFdEdx_patch::InitPar(const char *pardir,  const char *param)
{
  // expensive function - call it once/run
  Define_PMTsat();
//  for(int i=0; i<48; i++) ReadParTD( i, Form("%s/TD_%d_200b.txt",pardir,i) );
//  ReadParAtt(            Form("%s/attenuation.txt"              , pardir) );

  ReadParAtt(pardir, param);

// Old ToFNaNuclei files
  ReadParPos(            Form("%s/desaturation_position.txt"    , pardir) );
  ReadParDesatBB(        Form("%s/desaturation_beta.txt"        , pardir) );
//  ReadParBBneg(          Form("%s/BetheBloch.txt"               , pardir) );
//  ReadParBBpos(          Form("%s/BetheBloch_betagt1.txt"       , pardir) );

// New Bethe-Bloch pol4 file 
  ReadParBBpol4(          Form("%s/BetheBloch_pol4_tri.txt"       , pardir) );


};


//------------------------------------------------------------------------
void ToFdEdx_patch::Process( PamLevel2 *l2p, const char* alg, const char  *tri_or_bi)
{
  // the parameters should be already initialised by InitPar()

  Float_t dedx_corrp[48], dedx_corrhe[48], dedx_corrc[48];  
  Float_t thelp1,thelp2;
  Float_t xv[6],yv[6]; //,zin[6];

  trkAlg = alg;
  TF1 *ftri = new TF1("ftri", "1++x++x*x", 0,100);

//  cout<<"In PROCESS "<<endl;

  Clear();
  L2 = l2p;
  if(!L2)                                   return;
  if(!L2) cout<<" L2  error !!!"<<endl;

//  cout<<" L2->GetNTracks(trkAlg) check  "<<L2->GetNTracks(trkAlg)<<endl;

  
//  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  TrkLevel2 *trkL2 = L2->GetTrkLevel2();
  if(!trkL2) cout<<" L2->GetTrkLevel2() error !!!"<<endl;
  if(!trkL2)                                return;
  ToFLevel2 *tofL2 = L2->GetToFLevel2();
  if(!tofL2) cout<<" L2->GetToFLevel2() error !!!"<<endl;
  if(!tofL2)                                return;
//  eNtr = trkL2->GetNTracks();
//  if(eNtr<1) cout<<" eNtr < 1 !!!"<<endl;
//  if(eNtr<1)                                return; 

  eNtr = L2->GetNTracks(trkAlg);
  if(eNtr<1) cout<<" eNtr < 1 !!!"<<endl;
  if(eNtr<1)                                return; 



//  cout<<"In PROCESS after first checks "<<endl;


//  OBT   = L2->GetOrbitalInfo()->OBT;
//  PKT   = L2->GetOrbitalInfo()->pkt_num;
//  atime = L2->GetOrbitalInfo()->absTime;

  if ( L2->IsORB() ){
  OBT   = L2->GetOrbitalInfo()->OBT;
  PKT   = L2->GetOrbitalInfo()->pkt_num;
  atime = L2->GetOrbitalInfo()->absTime; 
  } else {
  atime = 1153000000.;  // simulated data no OrbitalInfo  fixed date 15-jul-2006
  };

//----------------------------------------------------------------------
//----------------------------------------------------------------------
//  Check Attenuation correction file time interval
// tmin_atten  & tmax_atten  time interval limits for attenuation parameters
//----------------------------------------------------------------------

//cout<<atime<<" "<<tmin_atten<<" "<<tmax_atten<<endl;

if ((atime < tmin_atten) || (atime > tmax_atten))  {
ical_atten = 0;
for (Int_t ii=0; ii<100; ii++) {
//cout<<ii<<"  tmin "<<T_int_min[ii]<<"  tmax "<<T_int_max[ii]<<endl;
if ((atime>T_int_min[ii]) && (atime<=T_int_max[ii]))  {  // time interval found
ical_atten = ii;
tmin_atten = T_int_min[ii];
tmax_atten = T_int_max[ii];
cout<<"Time Interval found: "<<ii<<"  tmin "<<T_int_min[ii]<<"  tmax "<<T_int_max[ii]<<endl;
break;
}
                             } // ii

     TArrayF &Apar0 =  A0_array[ical_atten];
     TArrayF &Apar1 =  A1_array[ical_atten];
     TArrayF &Apar2 =  A2_array[ical_atten];
     TArrayF &Apar3 =  A3_array[ical_atten];
     TArrayF &Apar4 =  A4_array[ical_atten];
     TArrayF &Apar5 =  A5_array[ical_atten];

       for (Int_t ii=0; ii<48;ii++) {
     A0[ii]=Apar0[ii];
     A1[ii]=Apar1[ii];
     A2[ii]=Apar2[ii];
     A3[ii]=Apar3[ii];
     A4[ii]=Apar4[ii];
     A5[ii]=Apar5[ii];
//     cout<<A0[ii]<<" "<<A1[ii]<<" "<<A2[ii]<<" "<<A3[ii]<<endl;
                                    }

cout<<"attenuation corr: abs time "<<atime<<" interval "<<ical_atten<<" tmin "<<T_int_min[ical_atten]<<" tmax "<<T_int_max[ical_atten]<<endl;
                                                  }                               
                        
//----------------------------------------------------------------------
// tmin_2nd  & tmax_2nd  time interval limits for 2nd order correction
//----------------------------------------------------------------------

//cout<<atime<<" "<<tmin_2nd<<" "<<tmax_2nd<<endl;

if ((atime < tmin_2nd) || (atime > tmax_2nd))  {

ical_2nd=0;
for (Int_t ii=0; ii<3500; ii++) {
//cout<<ii<<" "<<atime<<"  tmin "<<mtime[ii]<<"  tmax "<<mtime[ii+1]<<endl;
if ((atime>mtime[ii]) && (atime<=mtime[ii+1]))  {  // time interval found
ical_2nd = ii;
tmin_2nd = mtime[ii];
tmax_2nd = mtime[ii+1];
//cout<<"Time Interval 2nd order corr  "<<ii<<"  tmin "<<mtime[ii]<<"  tmax "<<mtime[ii+1]<<endl;
break;
}
                                  };
cout<<"2nd order corr: abs time "<<atime<<" interval "<<ical_2nd<<" tmin "<<mtime[ical_2nd]<<" tmax "<<mtime[ical_2nd+1]<<endl;
                                            } 
//-----------------------------------------------------------------------

//== interpolate betwen time limits

       thelp1 = mtime[ical_2nd];
       thelp2 = mtime[ical_2nd+1];

     TArrayF &corr1 =  dedx_corr_mp[ical_2nd];
     TArrayF &corr2 =  dedx_corr_mp[ical_2nd+1];

     TArrayF &corr3 =  dedx_corr_mhe[ical_2nd];
     TArrayF &corr4 =  dedx_corr_mhe[ical_2nd+1];

     TArrayF &corr5 =  dedx_corr_mc[ical_2nd];
     TArrayF &corr6 =  dedx_corr_mc[ical_2nd+1];

       for (Int_t ii=0; ii<48;ii++) {
       Float_t yhelp1 = corr1[ii];
       Float_t yhelp2 = corr2[ii];
       Float_t slope  = (yhelp2-yhelp1)/(thelp2-thelp1);
       Float_t inter  = yhelp1 - slope*thelp1;
       dedx_corrp[ii] = slope*atime + inter;
//       if (ii==0) cout<<thelp1<<" "<<thelp2<<" "<<tabs<<" "<<yhelp1<<" "<<yhelp2<<" "<<dedx_corr[0]<<endl;
//       cout<<ii<<" "<<yhelp1<<" "<<yhelp2<<" => "<<dedx_corr[ii]<<endl;

        yhelp1 = corr3[ii];
        yhelp2 = corr4[ii];
        slope  = (yhelp2-yhelp1)/(thelp2-thelp1);
        inter  = yhelp1 - slope*thelp1;
       dedx_corrhe[ii] = slope*atime + inter;

        yhelp1 = corr5[ii];
        yhelp2 = corr6[ii];
        slope  = (yhelp2-yhelp1)/(thelp2-thelp1);
        inter  = yhelp1 - slope*thelp1;
       dedx_corrc[ii] = slope*atime + inter;

                       }


//--------------- TABLE OF PERIODS WITH HV PROBLEMS --------------------------

    int Aconn=0;    // PMT 0,20,22,24
    int Bconn=0;    // PMT 6,12,26,34
    int Cconn=0;    // PMT 4,14,28,32
    int Dconn=0;    // PMT 2,8,10,30
    int Econn=0;    // PMT 42,43,44,47
    int Fconn=0;    // PMT 7,19,23,27
    int Gconn=0;    // PMT 3,11,25,33
    int Hconn=0;    // PMT 1,9,13,21
    int Iconn=0;    // PMT 5,29,31,35
    int Lconn=0;    // PMT 37,40,45,46
    int Mconn=0;    // PMT 15,16,17,18
    int Nconn=0;    // PMT 36,38,39,41

    int standard=0;

    int S115B_ok=0;
    int S115B_break=0;

  if(atime>=1153660001 && atime<=1154375000)Dconn=1;
    else if(atime>=1155850001 && atime<=1156280000){
      Hconn=1;
      Nconn=1;
    }

 else if(atime>=1158160000 && atime<=1158220000)Cconn=1;  // new period found WM 11-jan-2018

 else if(atime>=1168490001 && atime<=1168940000)Dconn=1;
    else if(atime>=1168940001 && atime<=1169580000){
      Fconn=1;
      Mconn=1;
    }

 else if(atime>=1174665001 && atime<=1175000000)Bconn=1;
    else if(atime>=1176120001 && atime<=1176800000)Hconn=1;
    else if(atime>=1176800001 && atime<=1178330000)Econn=1;
    else if(atime>=1178330001 && atime<=1181322000)Hconn=1;
    else if(atime>=1182100001 && atime<=1183030000)Aconn=1;
    else if(atime>=1184000001 && atime<=1184570000)Hconn=1;
    else if(atime>=1185090001 && atime<=1185212000)Dconn=1;
    else if(atime>=1191100001 && atime<=1191940000)Dconn=1;
    else if(atime>=1196230001 && atime<=1196280000)Hconn=1;
    else if(atime>=1206100001 && atime<=1206375600)Cconn=1;
    else if(atime>=1217989201 && atime<=1218547800)Econn=1;
    else if(atime>=1225789201 && atime<=1226566800)Econn=1;
    else if(atime>=1229400901 && atime<=1229700000)Econn=1;
    else if(atime>=1230318001 && atime<=1230415200)Econn=1;

    else if(atime>=1320710401 && atime<=1322006400)Dconn=1;  // (0,"","",1320710401,1322006400,213," D connector ",NULL); 
    else if(atime>=1309910401 && atime<=1310083200)Hconn=1;  // (0,"","",1309910401,1310083200,217," H connector ",NULL);
    else if(atime>=1384700001 && atime<=1387400000)Hconn=1;  // (0,"","",1384700001,1387400000,217," H connector ",NULL);
    else if(atime>=1230318001 && atime<=1230410000)Hconn=1;  //wm period 38 HV_H  1230318001 - 1230410000    12-26 - 12-27

    else if(atime>=1426620001 && atime<=1427745000)Dconn=1;  // 15-03-17 - 15-03-30   HV_D problems
    else if(atime>=1433820001 && atime<=1433960000)Dconn=1;  // 15-06-09 - 15-06-10   HV_D problems
    else if(atime>=1434460001 && atime<=1435310000)Dconn=1;  // 15-06-17 - 15-06-26   HV_D problems


    else {
      standard=1;
    }
    
    if(atime<1158720000)S115B_ok=1;
    else S115B_break=1;


//================================================================

// cout<<" L2->GetNTracks(trkAlg) "<<L2->GetNTracks(trkAlg)<<endl;

if( L2->GetNTracks(trkAlg)>=1 ){

//  PamTrack *trackTRK = L2->GetTrack(itr);     // 9th
  PamTrack *trackTRK = L2->GetTrack(0,trkAlg);  // 10th

  Float_t betamean = fabs(trackTRK->GetToFTrack()->beta[12]);

  Float_t def =trackTRK->GetExtTrack()->al[4];
  Float_t chi =trackTRK->GetExtTrack()->chi2;
  Float_t rig = 1./def;
    
  if(betamean<0.05 || betamean>2){
    for(int i=0;i<48;i++)INFOpmt[i]=1;
    for(int i=0;i<6;i++)INFOlayer[i]=1;
  }

 // define angle:   

  double dx   = trackTRK->GetToFTrack()->xtr_tof[0] - trackTRK->GetToFTrack()->xtr_tof[5];
  double dy   = trackTRK->GetToFTrack()->ytr_tof[0] - trackTRK->GetToFTrack()->ytr_tof[5];
  double dr   = sqrt(dx*dx+dy*dy);
  double theta=atan(dr/77.31);

  //  TArrayF adc;
  Float_t adc[48];
  for(int i=0;i<48;i++){
    adc[i]=0;
  }


// Clear the dEdx arrays
  for(int i=0;i<48;i++){
  eDEDXpmt[i] = 1000;
                       }

  for(int i=0;i<24;i++) eDEDXpad[i]   = 1000.;                           
  for(int i=0;i<6;i++)  eDEDXlayer[i] = 1000.;    

//  Raw ADC in the ToFPMT class
  
  for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
    ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
    Int_t pmt_id   = tofpmt->pmt_id;
    adc[pmt_id]    = tofpmt->adc ;
  //cout<<ipmt<<" "<<pmt_id<<" "<<adc[pmt_id]<<endl;
  } //loop on 48 pmts
  
 Int_t iverbose=0;

  
for (Int_t jj=0; jj<6; jj++){
xv[jj]=trackTRK->GetToFTrack()->xtr_tof[jj];
yv[jj]=trackTRK->GetToFTrack()->ytr_tof[jj];
//zin[jj] = L2->GetToFLevel2()->GetZTOF(L2->GetToFLevel2()->GetToFPlaneID(jj));
//cout<<"xv "<<xv[jj]<<" yv "<<yv[jj]<<endl;
}


Int_t paddleidoftrack[6]= {-1} ;
for (Int_t  ilay=0; ilay<6; ilay++) {
paddleidoftrack[ilay] = L2->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
//cout<<ilay<<" "<<paddleidoftrack[ilay]<<endl;
                                    }

//-----------------------------------------------------------------------------------

// Old loop ToFTrack  class
//  for (Int_t ipmt=0; ipmt<trackTRK->GetToFTrack()->npmtadc; ipmt++){
//    Int_t ii = trackTRK->GetToFTrack()->pmtadc[ipmt];

// New loop over PMTs found with GetPaddleIdOfTrack
  for (Int_t ii=0; ii<48; ii++) {

  Int_t ilay;
  if(ii<16)            ilay=0;
  if((15<ii)&&(ii<28)) ilay=1;
  if((27<ii)&&(ii<32)) ilay=2;
  if((31<ii)&&(ii<36)) ilay=3;
  if((35<ii)&&(ii<42)) ilay=4;
  if((41<ii)&&(ii<48)) ilay=5;

  if (paddleidoftrack[ilay] > -1) {   // paddle hitted
   
   Float_t xpos=0;
  if(ii<16)            xpos = yv[0];
  if((15<ii)&&(ii<28)) xpos = xv[1];
  if((27<ii)&&(ii<32)) xpos = xv[2];
  if((31<ii)&&(ii<36)) xpos = yv[3];
  if((35<ii)&&(ii<42)) xpos = yv[4];
  if((41<ii)&&(ii<48)) xpos = xv[5];


    if(adc[ii]==4095)INFOpmt[ii]=2;
    else if(adc[ii]>=PMTsat[ii]-5)INFOpmt[ii]=3;

    
    if( adc[ii] >= PMTsat[ii]-5 )  continue;
    if( adc[ii] <= 0. )            continue;
    
    double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
    double adccorr = adcpC*fabs(cos(theta));    
    
    if(adccorr<=0)INFOpmt[ii]=4;
         if(adccorr<=0.)           continue;


double adcHe, adcnorm, adclin, dEdx, Zeta;

 adcHe=-2;
 adcnorm=-2;
 adclin=-2;
 dEdx=-2;
 Zeta=-2;


 float dEdxH=-2;
 double dEdxHe=-2;

 double day = (atime-1150000000)/84600;


    Double_t correction = 1.;

    if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
      correction = 1.675;
    }
    else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
      correction = 2.482;
    }
    else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
      correction = 1.464;
    }
    else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
      correction = 1.995;
    }
    else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
      correction = 1.273;
    }
    else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
      correction = 1.565;
    }
    else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
      correction = 1.565;
    }
//    else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
//      correction = 1.018;
//    }
//-- new WM   12-jan-2018    
    else if(Nconn==1 && (ii==36 || ii==38 || ii==39)){
      correction = 1.34;
    }
    else if(Nconn==1 && (ii==41)){
      correction = 1.0;
    }
//-- end new WM        

    else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
      correction = 1.84;
    }
    else if(S115B_break==1 && ii==9 && Hconn==1){
      correction = 1.64;
    }
    else if(S115B_break==1 && ii==9 && Hconn==0){
      correction = 1.0;
    }
    else correction = 1.;

//   adcHe is the value of the pol5 as a function of the incident position, function derived for helium 
       
    if( ii==9 && S115B_break==1 ){
      adcHe   = (f_att5B( xpos ))/correction;
    } else {
    adcHe   = Get_adc_he(ii, xpos) / correction;
    };

    if(adcHe<=0)   continue;

//  adcnorm is f_pos(adccorr), where the f_pos function describes the nonlinearity of the PMT

    if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
    else adcnorm = f_pos( (parPos[ii]), adccorr);
  
    if(adcnorm<=0) continue;
    
    // adclin  is the linear calculation of the dEdx using adcnorm and the attenuation function for He

    if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
    else  adclin  = 4.*adcnorm/adcHe;

    if(adclin<=0) INFOpmt[ii]=5;
    if(adclin<=0)  continue;

// dEdx = f_desatBB(adclin)  uses a nonlinear function which takes Birks' effect etc. into account

    if(ii==9 && S115B_break==1) {
     dEdx = f_desatBB5B( adclin );
                                 }
    else
    {                          
   dEdx = f_desatBB((parDesatBB[ii]), adclin );
// new linear behaviour for some PMTs...
   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)) {
   if (adclin<10.) {
   Float_t dEdx_help = f_desatBB((parDesatBB[ii]), 10. );
   Float_t mhelp = dEdx_help / 10.;
   dEdx = mhelp * adclin ;
                   }
   }
   } // else

    if(dEdx<0) INFOpmt[ii]=6;
    if(dEdx<=0)    continue;

//--------------------- 2nd order correction ----------------------------------

//------------------------  bi-scale  -----------------------------------------
if ( strcmp(tri_or_bi,"BI") == 0){

     Float_t slope,inter;

     if(dedx_corrhe[ii]>dedx_corrp[ii]) slope=3./(dedx_corrhe[ii]-dedx_corrp[ii]);
     else slope=1.;
     if(dedx_corrhe[ii]>dedx_corrp[ii]) inter=1. - (slope*dedx_corrp[ii]);
     else inter=0.;

//  For S115B_break dedx_corrp is NOT proton, but calibrated with carbon !!! 
//  dedx_corrhe  is calibrated with helium, but probably not very reliable
//  Do simple linear scaling:
    if(ii==9 && S115B_break==1) {
     if(dedx_corrp[ii]>10) slope=36./dedx_corrp[ii];
     else slope=1.;
     inter = 0.;
                                }
                                
     dEdxH = inter + slope*dEdx;

}  // bi-scale

//------------------------  tri-scale  -----------------------------------------

if ( strcmp(tri_or_bi,"TRI") == 0){

Float_t Xa[3],Ya[3];
Ya[0] = 1;
Ya[1] = 4;
Ya[2] = 36;
Xa[0] = dedx_corrp[ii];
Xa[1] = dedx_corrhe[ii];
Xa[2] = dedx_corrc[ii];
TGraph *g1 = new TGraph(3,Xa,Ya);
g1->Fit("ftri","q");
dEdxH = ftri->Eval(dEdx);
g1->Delete();
}

//---------------------------------------------------------------------------------

//  if (iverbose==1)   cout<<"dEdx "<<dEdx<<"  dEdxH "<<dEdxH<<endl;

     if(dEdxH!=-2)eDEDXpmt[ii]=(Float_t)dEdxH;
     else eDEDXpmt[ii]=(Float_t)dEdx;

//  The f_BB=f(beta) correction is only used to calculate the charge Zeta

    if(dEdxH!=-2)dEdx=(Float_t)dEdxH;

    dEdxHe=-2;

    if(betamean<100) {

    if(ii==9 && S115B_break==1){
      if( betamean <1. ) dEdxHe = f_BB5Bpol4( betamean );
      else               dEdxHe = f_BB5Bpol4( 1.0 );
    } else {
      if( betamean <1. ) dEdxHe = f_BBpol4( (parBBpol4[ii]), betamean );
      else               dEdxHe = f_BBpol4( (parBBpol4[ii]), 1.0 );
    }

    if(ii==9 && S115B_break==1)  Zeta = sqrt(36.*(dEdx/dEdxHe)); 
    else  Zeta = sqrt(4.*(dEdx/dEdxHe));

    if(Zeta<0) INFOpmt[ii]=7;
   
    eZpmt[ii]=(Float_t)Zeta;

                       }  // betamean < 100




// 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 );

// if (ii==9)  cout<<"eDEDXpmt[ii] "<<eDEDXpmt[ii]<<endl; //" eZpmt[ii] "<<eZpmt[ii]<<endl; 

  } // paddleidoftrack > -1
 }  //end loop on PMTs 0-47


//--------------------------  paddle + layer ----------------------
//-----------------------------------------------------------------
// The GetdEdx method "handmade"
//-----------------------------------------------------------------

Int_t ipad_a[6] = {0,8,14,16,18,21};
Int_t ipmt_a[6] = {0,16,28,32,36,42};

for (Int_t  ilay=0; ilay<6; ilay++) {
paddleidoftrack[ilay] = L2->GetToFLevel2()->GetPaddleIdOfTrack(xv[ilay], yv[ilay], ilay, 0.0) ;
                                    }


for (Int_t ilay=0; ilay<6; ilay++) {
Float_t  xhelp=1000.;
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
Int_t ipad  = ipad_a[ilay]+i1; 
Int_t ihelp = ipmt_a[ilay]+(i1*2); 
if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]==1000.)) xhelp = eDEDXpmt[ihelp];
if ((eDEDXpmt[ihelp+1]<1000.) && (eDEDXpmt[ihelp]==1000.)) xhelp = eDEDXpmt[ihelp+1];
if ((eDEDXpmt[ihelp]<1000.) && (eDEDXpmt[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmt[ihelp]+eDEDXpmt[ihelp+1]);
eDEDXpad[ipad] = xhelp;                           
eDEDXlayer[ilay] = xhelp;                           

if ((eZpmt[ihelp]>-1) && (eZpmt[ihelp+1]==-1)) xhelp = eZpmt[ihelp];
if ((eZpmt[ihelp+1]>-1) && (eZpmt[ihelp]==-1)) xhelp = eZpmt[ihelp+1];
if ((eZpmt[ihelp]>-1) && (eZpmt[ihelp+1]>-1)) xhelp = 0.5*(eZpmt[ihelp]+eZpmt[ihelp+1]);
eZpad[ipad] = xhelp;                           
eZlayer[ilay] = xhelp;                           


                        }                       
                                 } // ilay


  } // if( L2->GetNTracks(trkAlg)>=1 )...



//================================================================
//================================================================
//==============        standalone dEdx and beta     =============
//================================================================
//================================================================

if (L2->GetToFStoredTrack(-1)  ) {

// ToF track

  ToFTrkVar  *tof = L2->GetToFStoredTrack(-1);

  Float_t betamean = tof->CalcBeta(10., 10., 20.);

  if(betamean<0.05 || betamean>2){
    for(int i=0;i<48;i++)INFOpmtstd[i]=1;
    for(int i=0;i<6;i++)INFOlayerstd[i]=1;
  }

 // define angle:   

double theta =0.;
if (fabs((tof->xtofpos[0])<100) && (fabs(tof->xtofpos[2])<100) && (fabs(tof->ytofpos[0])<100) && (fabs(tof->ytofpos[2])<100)) {
  double dx   = tof->xtofpos[0] - tof->xtofpos[2];
  double dy   = tof->ytofpos[0] - tof->ytofpos[2];
  double dr   = sqrt(dx*dx+dy*dy);
  theta=atan(dr/77.31);
                                                               } // if...

  //  Clear ADC;
  Float_t adc[48];
  for(int i=0;i<48;i++){
//    adc[i]=0;
    adc[i]=4095;
  }
  
  
  for (Int_t ipmt=0; ipmt<tofL2->npmt() ; ipmt++){
    ToFPMT *tofpmt = tofL2->GetToFPMT(ipmt);
    Int_t pmt_id   = tofpmt->pmt_id;
    adc[pmt_id]    = tofpmt->adc ;  // raw ADC
  } //loop on 48 pmts
  
  
//====================================================================
//===========  Check ToF standalone using HitPaddle ==================
//====================================================================

Int_t paddleidoftrack[6] = {-1, -1, -1, -1, -1, -1};;

for(Int_t jj=0; jj<8; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(0,jj);
if (HitPad==1)  paddleidoftrack[0] = jj;
                           }
for(Int_t jj=0; jj<6; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(1,jj);
if (HitPad==1)  paddleidoftrack[1] = jj;
                           }

for(Int_t jj=0; jj<2; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(2,jj);
if (HitPad==1)   paddleidoftrack[2] = jj;
                           }
for(Int_t jj=0; jj<2; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(3,jj);
if (HitPad==1)   paddleidoftrack[3] = jj;
                           }

for(Int_t jj=0; jj<3; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(4,jj);
if (HitPad==1)   paddleidoftrack[4] = jj;
                           }
for(Int_t jj=0; jj<3; jj++){
Int_t HitPad = L2->GetToFLevel2()->HitPaddle(5,jj);
if (HitPad==1)   paddleidoftrack[5] = jj;
                           }

/*
cout<<"paddleidoftrack  std ";
for (Int_t  ilay=0; ilay<6; ilay++) {
cout<<" "<<paddleidoftrack[ilay];
                                    }
cout<<endl;
*/
//==================================================================
//==================================================================
//==================================================================

// Clear the dEdx arrays
  for(int i=0;i<48;i++){
  eDEDXpmtstd[i] = 1000;
                       }

  for(int i=0;i<24;i++) eDEDXpadstd[i]   = 1000.;                           
  for(int i=0;i<6;i++)  eDEDXlayerstd[i] = 1000.;                           


Int_t ipad_a[6] = {0,8,14,16,18,21};
Int_t ipmt_a[6] = {0,16,28,32,36,42};

Int_t pmthit[48];
for(int i=0;i<48;i++) pmthit[i]   = 0;                           

for(int i=0;i<6;i++) {
if (paddleidoftrack[i] > -1) {
Int_t ipad = ipad_a[i] + paddleidoftrack[i];
pmthit[2*ipad] = 1;
pmthit[2*ipad+1] = 1;
                             }
}



  for (Int_t ipmt=0; ipmt<48; ipmt++){
   if (pmthit[ipmt] == 1) {
  
  Int_t ii = ipmt;

  Float_t xpos=0;
  if(ii<16)            xpos = tof->ytofpos[0];
  if((15<ii)&&(ii<28)) xpos = tof->xtofpos[0];
  if((27<ii)&&(ii<32)) xpos = tof->xtofpos[1];
  if((31<ii)&&(ii<36)) xpos = tof->ytofpos[1];
  if((35<ii)&&(ii<42)) xpos = tof->ytofpos[2];
  if((41<ii)&&(ii<48)) xpos = tof->xtofpos[2];
     
    
    if(adc[ii]==4095)INFOpmtstd[ii]=2;
    else if(adc[ii]>=PMTsat[ii]-5)INFOpmtstd[ii]=3;

    //cout<<ii<<" "<<adc[ii]<<" "<<PMTsat[ii]<<endl;
    
    if( adc[ii] >= PMTsat[ii]-5 )  continue;
    if( adc[ii] <= 0. )            continue;
    
    double adcpC   = f_adcPC( adc[ii] );    // - adc conversion in pC
    double adccorr = adcpC*fabs(cos(theta));
    
     //cout<<ii<<" adc "<<adc[ii]<<" adcpC "<<adcpC<<" adccorr "<<adccorr<<endl;
 
    
    if(adccorr<=0)INFOpmt[ii]=4;
         if(adccorr<=0.)           continue;


double adcHe, adcnorm, adclin, dEdx, Zeta;

 adcHe=-2;
 adcnorm=-2;
 adclin=-2;
 dEdx=-2;
 Zeta=-2;

 float dEdxH=-2;
 double dEdxHe=-2;

// double day = (atime-1150000000)/84600;


    Double_t correction = 1.;

    if(Aconn==1 && (ii==0 || ii==20 || ii==22 || ii==24)){
      correction = 1.675;
    }
    else if(Bconn==1 && (ii==6 || ii==12 || ii==26 || ii==34)){
      correction = 2.482;
    }
    else if(Cconn==1 && (ii==4 || ii==14 || ii==28 || ii==32)){
      correction = 1.464;
    }
    else if(Dconn==1 && (ii==2 || ii==8 || ii==10 || ii==30)){
      correction = 1.995;
    }
    else if(Econn==1 && (ii==42 || ii==43 || ii==44 || ii==47)){
      correction = 1.273;
    }
    else if(Fconn==1 && (ii==7 || ii==19 || ii==23 || ii==27)){
      correction = 1.565;
    }
    else if(Mconn==1 && (ii==15 || ii==16 || ii==17 || ii==18)){
      correction = 1.565;
    }
//    else if(Nconn==1 && (ii==36 || ii==38 || ii==39 || ii==41)){
//      correction = 1.018;
//    }
//-- new WM   12-jan-2018    
    else if(Nconn==1 && (ii==36 || ii==38 || ii==39)){
      correction = 1.34;
    }
    else if(Nconn==1 && (ii==41)){
      correction = 1.0;
    }
//-- end new WM        

    else if(Hconn==1 && (ii==1 || ii==13 || ii==21 || (ii==9&&S115B_ok==1))){
      correction = 1.84;
    }
    else if(S115B_break==1 && ii==9 && Hconn==1){
      correction = 1.64;
    }
    else if(S115B_break==1 && ii==9 && Hconn==0){
      correction = 1.0;
    }
    else correction = 1.;

//   adcHe is the value of the pol5 as a function of the incident position, curve derived for helium 
    
    if( ii==9 && S115B_break==1 ){
      adcHe   = (f_att5B( xpos ))/correction;
    } else {
    adcHe   = Get_adc_he(ii, xpos) / correction;
    };

    if(adcHe<=0)   continue;

//  adcnorm is f_pos(adccorr), where the f_pos function describes the nonlinearity of the PMT

    if(ii==9 && S115B_break==1)  adcnorm = f_pos5B(adccorr);
    else adcnorm = f_pos( (parPos[ii]), adccorr);
  
    //cout<<"adcnorm "<<adcnorm<<endl;

    if(adcnorm<=0) continue;
    
    // adclin  is the linear calculation of the dEdx using adcnorm and the attenuation function for He

    if(ii==9 && S115B_break==1)  adclin  = 36.*adcnorm/adcHe;
    else  adclin  = 4.*adcnorm/adcHe;

    //cout<<"adclin "<<adclin<<endl;

    if(adclin<=0) INFOpmt[ii]=5;
    if(adclin<=0)  continue;

// dEdx = f_desatBB(adclin)  uses a nonlinear function which takes Birks' effect etc. into account


    if(ii==9 && S115B_break==1) {
     dEdx = f_desatBB5B( adclin );
                                 }
    else
    {                          
   dEdx = f_desatBB((parDesatBB[ii]), adclin );
// new linear behaviour for some PMTs...
   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)) {
   if (adclin<10.) {
   Float_t dEdx_help = f_desatBB((parDesatBB[ii]), 10. );
   Float_t mhelp = dEdx_help / 10.;
   dEdx = mhelp * adclin ;
                   }
//   if (iverbose==1) cout<<"dEdx after linear correction  f_desatBB  "<<dEdx<<endl;
   }
   } // else

    if(dEdx<0) INFOpmtstd[ii]=6;
    if(dEdx<=0)    continue;


//--------------------- 2nd order correction ----------------------------------

//------------------------  bi-scale  -----------------------------------------
if ( strcmp(tri_or_bi,"BI") == 0){

     Float_t slope,inter;

     if(dedx_corrhe[ii]>dedx_corrp[ii]) slope=3./(dedx_corrhe[ii]-dedx_corrp[ii]);
     else slope=1.;
     if(dedx_corrhe[ii]>dedx_corrp[ii]) inter=1. - (slope*dedx_corrp[ii]);
     else inter=0.;

//  For S115B_break dedx_corrp is NOT proton, but calibrated with carbon !!! 
//  dedx_corrhe  is calibrated with helium, but probably not very reliable
//  Do simple linear scaling:
    if(ii==9 && S115B_break==1) {
     if(dedx_corrp[ii]>10) slope=36./dedx_corrp[ii];
     else slope=1.;
     inter = 0.;
                                }
                                
     dEdxH = inter + slope*dEdx;

}  // bi-scale

//------------------------  tri-scale  -----------------------------------------

if ( strcmp(tri_or_bi,"TRI") == 0){

Float_t Xa[3],Ya[3];
Ya[0] = 1;
Ya[1] = 4;
Ya[2] = 36;
Xa[0] = dedx_corrp[ii];
Xa[1] = dedx_corrhe[ii];
Xa[2] = dedx_corrc[ii];
TGraph *g1 = new TGraph(3,Xa,Ya);
g1->Fit("ftri","q");
dEdxH = ftri->Eval(dEdx);
g1->Delete();
}

//---------------------------------------------------------------------------------

//  if (iverbose==1)   cout<<"dEdx "<<dEdx<<"  dEdxH "<<dEdxH<<endl;

     if(dEdxH!=-2)eDEDXpmtstd[ii]=(Float_t)dEdxH;
     else eDEDXpmtstd[ii]=(Float_t)dEdx;
     

//  The f_BB=f(beta) correction is only used in ToFNaNuclei to calculate the charge Zeta

    dEdxHe=-2;

    if(betamean<100) {

    if(ii==9 && S115B_break==1){
      if( betamean <1. ) dEdxHe = f_BB5Bpol4( betamean );
      else               dEdxHe = f_BB5Bpol4( 1.0 );
    } else {
      if( betamean <1. ) dEdxHe = f_BBpol4( (parBBpol4[ii]), betamean );
      else               dEdxHe = f_BBpol4( (parBBpol4[ii]), 1.0 );
    }

    if(ii==9 && S115B_break==1)  Zeta = sqrt(36.*(dEdx/dEdxHe)); 
    else  Zeta = sqrt(4.*(dEdx/dEdxHe));

    if(Zeta<0) INFOpmtstd[ii]=7;
   
    eZpmtstd[ii]=(Float_t)Zeta;

                       }  // betamean < 100



  //   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 );

     //cout<<"eDEDXpmtstd[ii] "<<eDEDXpmt[ii]<<" eZpmtstd[ii] "<<eZpmtstd[ii]<<endl; 

  } // hitted PMTs
 }  //end loop on PMTs along the track

//-----------------------------------------------------------------
// The GetdEdx method "handmade"
//-----------------------------------------------------------------


for (Int_t ilay=0; ilay<6; ilay++) {
Float_t  xhelp=1000.;
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
Int_t ipad  = ipad_a[ilay]+i1; 
Int_t ihelp = ipmt_a[ilay]+(i1*2); 
if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]==1000.)) xhelp = eDEDXpmtstd[ihelp];
if ((eDEDXpmtstd[ihelp+1]<1000.) && (eDEDXpmtstd[ihelp]==1000.)) xhelp = eDEDXpmtstd[ihelp+1];
if ((eDEDXpmtstd[ihelp]<1000.) && (eDEDXpmtstd[ihelp+1]<1000.)) xhelp = 0.5*(eDEDXpmtstd[ihelp]+eDEDXpmtstd[ihelp+1]);
eDEDXpadstd[ipad] = xhelp;                           
eDEDXlayerstd[ilay] = xhelp;                           

if ((eZpmtstd[ihelp]>-1) && (eZpmtstd[ihelp+1]==-1)) xhelp = eZpmtstd[ihelp];
if ((eZpmtstd[ihelp+1]>-1) && (eZpmtstd[ihelp]==-1)) xhelp = eZpmtstd[ihelp+1];
if ((eZpmtstd[ihelp]>-1) && (eZpmtstd[ihelp+1]>-1)) xhelp = 0.5*(eZpmtstd[ihelp]+eZpmtstd[ihelp+1]);
eZpadstd[ipad] = xhelp;                           
eZlayerstd[ilay] = xhelp;                           


                        }                       
                                 } // ilay

  } // if (pam_event->GetToFStoredTrack(-1)  )...

ftri->Delete();

};


//------------------------------------------------------------------------
void ToFdEdx_patch::PrintTD()
{
  for(int i=0; i<48; i++) {   
/*
    TArrayF &binx = TDx[i];
    TArrayF &biny = TDy[i];
    for(int k=0; k<200; k++) {  // bin temporali
      printf("%d %d %f %f", i,k, binx[k], biny[k]);
     
    }
*/
  }
}


//------------------------------------------------------------------------
void ToFdEdx_patch::Define_PMTsat()
{
  Float_t  sat[48] = {
    3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
    3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
    3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
    3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
    3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
    3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
  PMTsat.Set(48,sat);
}

//------------------------------------------------------------------------
/*
void ToFdEdx_patch::ReadParTD( Int_t ipmt, const char *fname )
{
  printf("read %s\n",fname);
  if(ipmt<0)  return;
  if(ipmt>47) return;
  FILE *fattin = fopen( fname , "r" );
  Float_t yTD[200],xTD[200];
  for(int j=0;j<200;j++){
    float x,y,ym,e;
    if(fscanf(fattin,"%f %f %f %f", 
	      &x, &y, &ym, &e )!=4) break;
    xTD[j]=x;
    if(ym>0&&fabs(y-ym)>1)  yTD[j]=ym;
    else                    yTD[j]=y;
  }
  TDx[ipmt].Set(200,xTD);
  TDy[ipmt].Set(200,yTD);
  fclose(fattin);
}
*/
//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParBBpos( const char *fname )
{
  printf("read %s\n",fname);
  parBBpos.Set(48);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp;
    if(fscanf(fattin,"%d %f", 
	      &tid, &tp )!=2) break;
    parBBpos[i]=tp;
  }
  fclose(fattin);
}

//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParDesatBB( const char *fname )
{
  printf("read %s\n",fname);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp[3];
    if(fscanf(fattin,"%d %f %f %f", 
	      &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
    parDesatBB[i].Set(3,tp);
  }
  fclose(fattin);
}


//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParBBneg( const char *fname )

{
  printf("read %s\n",fname);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp[3];
    if(fscanf(fattin,"%d %f %f %f", 
	      &tid, &tp[0], &tp[1], &tp[2] )!=4) break;
    parBBneg[i].Set(3,tp);
  }
  fclose(fattin);
}

//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParBBpol4( const char *fname )

{
  printf("read %s\n",fname);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp[5];
    if(fscanf(fattin,"%d %f %f %f %f %f", 
	      &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4] )!=6) break;
    parBBpol4[i].Set(5,tp);
  }
  fclose(fattin);
}
//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParPos( const char *fname )
{
  printf("read %s\n",fname);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp[4];
    if(fscanf(fattin,"%d %f %f %f %f", 
	      &tid, &tp[0], &tp[1], &tp[2], &tp[3])!=5) break;
    parPos[i].Set(4,tp);
  }
  fclose(fattin);
}

//------------------------------------------------------------------------
/*void ToFdEdx_patch::ReadParAtt( const char *fname )
{
  printf("read %s\n",fname);
  FILE *fattin = fopen( fname , "r" );
  for (int i=0; i<48; i++) {
    int   tid=0;
    float  tp[6];
    if(fscanf(fattin,"%d %f %f %f %f %f %f", 
	      &tid, &tp[0], &tp[1], &tp[2], &tp[3], &tp[4], &tp[5] )!=7) break;
    parAtt[i].Set(6,tp);
  }
  fclose(fattin);
}
*/

//---------------------------------------------------------------------------

double ToFdEdx_patch::f_att(float C0, float C1, float C2, float C3, float C4, float C5,  float x)
{

 //cout<<"in f_att: para "<<C0<<" "<<C1<<" "<<C2<<" "<<C3<<" "<<C4<<" "<<C5<<" x "<<x<<endl;

  return 
    C0 +
    C1*x +
    C2*x*x +
    C3*x*x*x +
    C4*x*x*x*x +
    C5*x*x*x*x*x;
}


/*
double ToFdEdx_patch::f_att( TArrayF &p, float x )
{

 //cout<<"in f_att: para "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" "<<p[3]<<" "<<p[4]<<" "<<p[5]<<" x "<<x<<endl;

  return 
    p[0] +
    p[1]*x +
    p[2]*x*x +
    p[3]*x*x*x +
    p[4]*x*x*x*x +
    p[5]*x*x*x*x*x;
}
*/

//------------------------------------------------------------------------
double ToFdEdx_patch::f_att5B( float x )
{
  return 
    101.9409 +
    6.643781*x +
    0.2765518*x*x +
    0.004617647*x*x*x +
    0.0006195132*x*x*x*x +
    0.00002813734*x*x*x*x*x;
} 


double ToFdEdx_patch::f_pos( TArrayF &p, float x )
{
  return 
    p[0] +
    p[1]*x +
    p[2]*x*x +
    p[3]*x*x*x;
}

double ToFdEdx_patch::f_pos5B( float x )
{
  return
    15.45132 +
    0.8369721*x +
    0.0005*x*x;
}



double ToFdEdx_patch::f_adcPC( float x )
{
  return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
}

/*
float ToFdEdx_patch::Get_adc_he( int id, float pl_x[6], float pl_y[6])
{
 
  //
  // input: id - pmt [0:47}
  //             pl_x - coord x of the tof plane
  //             pl_y - coord y



   adc_he = 0;
//  if( eGeom.GetXY(id)==1 )  adc_he = f_att( (parAtt[id]), pl_x[eGeom.GetPlane(id)] );
//  if( eGeom.GetXY(id)==2 )  adc_he = f_att( (parAtt[id]), pl_y[eGeom.GetPlane(id)] );
  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)] );
  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)] );

  //cout<<"in Get_adc_he "<<adc_he<<endl;

  return adc_he;
}
*/

//------------------------------------------------------------------------

float ToFdEdx_patch::Get_adc_he( int id, float x)
{
  // input: id - pmt [0:47] , position in the paddle

  adc_he = 0;
  adc_he = f_att( A0[id],A1[id],A2[id],A3[id], A4[id], A5[id], x );

  //cout<<"in Get_adc_he "<<adc_he<<endl;

  return adc_he;
}

//------------------------------------------------------------------------
double ToFdEdx_patch::f_BB( TArrayF &p, float x )
{
  return  p[0]/(x*x)*(log(x*x/(1-x*x)) - p[1]*x*x - p[2]);
}

//------------------------------------------------------------------------
double ToFdEdx_patch::f_BBpol4( TArrayF &p, float x )
{
  return  p[0] + p[1]*x + p[2]*x*x + p[3]*x*x*x + p[4 ]*x*x*x*x;
}

//------------------------------------------------------------------------
double ToFdEdx_patch::f_BB5B( float x )
{
  return  0.165797/(x*x)*(log(x*x/(1-x*x)) + 140.481*x*x + 52.9258);
}
//------------------------------------------------------------------------
double ToFdEdx_patch::f_BB5Bpol4( float x )
{
//  return  138.165 - 120.884*x - 122.502*x*x + 203.795*x*x*x - 64.1068*x*x*x*x ;
  return  67.2843 + 223.175*x - 721.068*x*x + 646.677*x*x*x - 181.868*x*x*x*x ;
 
}
//------------------------------------------------------------------------
double ToFdEdx_patch::f_desatBB( TArrayF &p, float x )
{
  return 
    p[0] +
    p[1]*x +
    p[2]*x*x;
}

//------------------------------------------------------------------------
double ToFdEdx_patch::f_desatBB5B( float x )
{
  return 
    -2.4 +
    0.75*x +
    0.009*x*x;
}

//------------------------------------------------------------------------
void ToFdEdx_patch::ReadParAtt(const char *pardir, const char *param)
{
  char filename[100],filename1[100],filename2[100],name[50];
  char str1[30];


//---------------------  flight ---------------------------

 if (strcmp(param,"simu") != 0){
 cout<<"Flight calibrations "<<endl;

  sprintf(filename,"%stime_intervals_atten_10th.txt",pardir);

  cout<<"read Time Interval File for Attenuation : "<<filename<<endl;

  FILE *ftimein = fopen( filename , "r" );

  for (int ii=0; ii<100; ii++) {
    UInt_t tp[2];
    if(fscanf(ftimein,"%s %u %u",
              &str1[0], &tp[0], &tp[1])!=3) break;
    cout<<str1<<" "<<tp[0]<<" "<<tp[1]<<endl;
    T_int_min[ii] = tp[0];
    T_int_max[ii] = tp[1];


  sprintf(name,"%s",str1);
//  cout<<"test name  "<<name<<endl;
  sprintf(filename1,"%s%s",pardir,name);

  cout<<"read Attenuation file "<<filename1<<endl;
  FILE *fattin = fopen( filename1, "r" );

  Float_t A0help[48];
  Float_t A1help[48];
  Float_t A2help[48];
  Float_t A3help[48];
  Float_t A4help[48];
  Float_t A5help[48];
  

  for (int i=0; i<48; i++) {
    int   id=0;
    float par[6];
    if(fscanf(fattin,"%d %f %f %f %f %f %f",
              &id, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5] )!=7) break;
    A0help[i] = par[0];
    A1help[i] = par[1];
    A2help[i] = par[2];
    A3help[i] = par[3];
    A4help[i] = par[4];
    A5help[i] = par[5];
  }
  
  A0_array[ii].Set(48,A0help);
  A1_array[ii].Set(48,A1help);
  A2_array[ii].Set(48,A2help);
  A3_array[ii].Set(48,A3help);
  A4_array[ii].Set(48,A4help);
  A5_array[ii].Set(48,A5help);

  
  fclose(fattin);


//  if (T_int_max[ii] == 2000000000)  break;  // end of file
  if (T_int_max[ii] > 2000000000)  break;  // end of file

  }
  fclose(ftimein);

} // flight

//---------------------  simu ---------------------------

 if ( strcmp(param,"simu") == 0){

  //cout<<"Simulation calibrations "<<endl;
  sprintf(name,"tofdedx_gp_atten.dat");
  sprintf(filename1,"%s%s",pardir,name);

  cout<<"read Attenuation file "<<filename1<<endl;
  FILE *fattin = fopen( filename1, "r" );

  Float_t A0help[48];
  Float_t A1help[48];
  Float_t A2help[48];
  Float_t A3help[48];
  Float_t A4help[48];
  Float_t A5help[48];
  

  for (int i=0; i<48; i++) {
    int   id=0;
    float par[6];
    if(fscanf(fattin,"%d %f %f %f %f %f %f",
              &id, &par[0], &par[1], &par[2], &par[3], &par[4], &par[5] )!=7) break;
    A0help[i] = par[0];
    A1help[i] = par[1];
    A2help[i] = par[2];
    A3help[i] = par[3];
    A4help[i] = par[4];
    A5help[i] = par[5];
  }
  
  A0_array[0].Set(48,A0help);
  A1_array[0].Set(48,A1help);
  A2_array[0].Set(48,A2help);
  A3_array[0].Set(48,A3help);
  A4_array[0].Set(48,A4help);
  A5_array[0].Set(48,A5help);

  
  fclose(fattin);

T_int_min[0] = 1000000000;
T_int_max[0] = 2000000000;


} // simu




// set first time interval
ical_atten = 0;
tmin_atten = T_int_min[0];
tmax_atten = T_int_max[0];

// set first time interval attenuation

     TArrayF &Apar0 =  A0_array[0];
     TArrayF &Apar1 =  A1_array[0];
     TArrayF &Apar2 =  A2_array[0];
     TArrayF &Apar3 =  A3_array[0];
     TArrayF &Apar4 =  A4_array[0];
     TArrayF &Apar5 =  A5_array[0];

       for (Int_t ii=0; ii<48;ii++) {
     A0[ii]=Apar0[ii];
     A1[ii]=Apar1[ii];
     A2[ii]=Apar2[ii];
     A3[ii]=Apar3[ii];
     A4[ii]=Apar4[ii];
     A5[ii]=Apar5[ii];
                                    }

cout<<"First time interval attenuation: "<<tmin_atten<<" - "<<tmax_atten<<endl;


//-----------------------------------------------------------
// Here I  read the "bi-scale" dEdx_korr parameters 2nd order
//-----------------------------------------------------------

//Double_t t1,t2,tm;
//UInt_t t1,t2,tm;
Double_t t1,t2,tm;

Float_t xm1,xm2,xm3;
Float_t pmthelpp[48];
Float_t pmthelphe[48];
Float_t pmthelpc[48];

//---------------------  flight ---------------------------
if ( strcmp(param,"simu") != 0){

sprintf(filename2,"%striscale_tofdedx_patch.txt",pardir);

cout<<"Filename 2nd order tri-scale "<<filename2<<endl;
ifstream fin(filename2);

Int_t jj=0;
Int_t idummy;
while (! fin.eof()) {
fin>>t1>>tm>>t2;
//cout<<jj<<" "<<tm<<endl;
cout << setiosflags(ios::fixed)  << setw(10) << setprecision(3) << tm << endl;
//cout <<t1<<" "<<tm<<" "<<t2<<endl;
mtime[jj]=tm;
for (Int_t ii=0; ii<48;ii++) {
fin>>xm1>>xm2>>xm3;  // He-corr, P-corr, C-corr
pmthelphe[ii] = xm1;
pmthelpp[ii]  = xm2;
pmthelpc[ii]  = xm3;
                       }
dedx_corr_mp[jj].Set(48,pmthelpp);
dedx_corr_mhe[jj].Set(48,pmthelphe);
dedx_corr_mc[jj].Set(48,pmthelpc);
jj=jj+1;
}

/*
Int_t idummy;
for (Int_t jj=0; jj<2000; jj++) {
fin>>t1>>tm>>t2;
cout<<jj<<" "<<t1<<" "<<tm<<" "<<t2<<endl;
if (t2==2000000000) break;
cout << setiosflags(ios::fixed)  << setw(10) << setprecision(3) << tm << endl;
mtime[jj]=tm;
for (Int_t ii=0; ii<48;ii++) {
fin>>idummy>>xmean1>>xwidth1;
dedx_corr_m[jj][ii]=xmean1;
                       }
}
*/

fin.close();

}  // flight

//---------------------  simu ---------------------------

if ( strcmp(param,"simu") == 0){

sprintf(filename2,"%striscale_tofdedx_patch_simu.txt",pardir);

cout<<"Filename 2nd order "<<filename2<<endl;
ifstream fin(filename2);

Int_t jj=0;
Int_t idummy;
while (! fin.eof()) {
fin>>t1>>tm>>t2;
//cout<<jj<<" "<<tm<<endl;
//cout << setiosflags(ios::fixed)  << setw(10) << setprecision(3) << tm << endl;
cout <<t1<<" "<<tm<<" "<<t2<<endl;
mtime[jj]=tm;
for (Int_t ii=0; ii<48;ii++) {
fin>>xm1>>xm2>>xm3;  // He-corr, P-corr, C-corr
pmthelphe[ii] = xm1;
pmthelpp[ii]  = xm2;
pmthelpc[ii]  = xm3;
                       }
dedx_corr_mp[jj].Set(48,pmthelpp);
dedx_corr_mhe[jj].Set(48,pmthelphe);
dedx_corr_mc[jj].Set(48,pmthelpc);
jj=jj+1;
}


fin.close();

}  // simu


// set first time interval
ical_2nd=0;
tmin_2nd = mtime[0];
tmax_2nd = mtime[1];

cout<<"First time interval 2nd order corr.: "<<tmin_2nd<<" - "<<tmax_2nd<<endl;


}


//----------------------------------------------------------------------------