#include <ToFPatch.h>
#include <TGraph.h>
#include <TSpline.h>
#include <TMVA/TSpline1.h>

//------------------------------------------------------------------------
void ToFPatch::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
  eDEDXpmtraw.Set(48); eDEDXpmtraw.Reset(-1); // Set array size and reset structure
  eDEDXpad.Set(24);    eDEDXpad.Reset(-1);
  eDEDXlayer.Set(6);   eDEDXlayer.Reset(-1);
  INFOpmt.Set(48);     INFOpmt.Reset(0);
  INFOlayer.Set(6);    INFOlayer.Reset(0);
  
  eDEDXpmtstd.Set(48);    eDEDXpmtstd.Reset(-1);   // Set array size  and reset structure
  eDEDXpmtrawstd.Set(48); eDEDXpmtrawstd.Reset(-1); // Set array size and reset structure
  eDEDXpadstd.Set(24);    eDEDXpadstd.Reset(-1);
  eDEDXlayerstd.Set(6);   eDEDXlayerstd.Reset(-1);
  INFOpmtstd.Set(48);     INFOpmtstd.Reset(0);
  INFOlayerstd.Set(6);    INFOlayerstd.Reset(0);

  
  //
};

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


//------------------------------------------------------------------------
void ToFPatch::InitPar(const char *pardir, const char *param)
{
  // expensive function - call it once/run
  Define_PMTsat();
  ReadParAtt(pardir, param);
};


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

  Float_t thelp1,thelp2;
  Float_t dedx_corr[48];
  Float_t xv[6],yv[6],zin[6];

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


  trkAlg = alg;

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

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

// tmin_atten  & tmax_atten  time interval limits for attenuation parameters
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];

       for (Int_t ii=0; ii<48;ii++) {
     A0[ii]=Apar0[ii];
     A1[ii]=Apar1[ii];
     A2[ii]=Apar2[ii];
     A3[ii]=Apar3[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<1500; 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_m[ical_2nd];
     TArrayF &corr2 =  dedx_corr_m[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_corr[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;
                       }


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

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

  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:   

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

Float_t theta[6];

for (Int_t jj=0; jj<3; jj++){
      Float_t dx=0.;
      Float_t dy=0.;
      Float_t dr=0.;
      Int_t kk=2*jj;

      theta[kk]  =0;
      theta[kk+1]=0;
      Float_t dist = zin[kk] - zin[kk+1]; 

      if ((xv[kk])<100.) {
         dx   = xv[kk]  - xv[kk+1];
         dy   = yv[kk]  - yv[kk+1];
         dr   = sqrt(dx*dx+dy*dy);
        theta[kk]   = atan(dr/dist);
        theta[kk+1] = atan(dr/dist);
      }

                           }  // jj

  //  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
  
  
 
Int_t paddleidoftrack[6]= {0} ;

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

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


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

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

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

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


//  ratio helium to proton ca. 4
    Float_t  hepratio = 4.;

Int_t ipad,ihelp;
Float_t adclraw,adcrraw;
Float_t adcl,adcr;
Float_t dedxl,dedxr;
Float_t xkorr;

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++) {

Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {

ipad  = ipad_a[ilay]+i1; 
ihelp = ipmt_a[ilay]+(i1*2); 

/*    
    if(adc[ihelp]==4095)INFOpmt[ihelp]=2;
    else if(adc[iihelp]>=PMTsat[iihelp]-5)INFOpmt[iihelp]=3;
    
    if( adc[ihelp] >= PMTsat[iihelp]-5 )  continue;
    if( adc[ihelp] <= 0. )            continue;
    
    double adcpC   = f_adcPC( adc[iihelp] );    // - adc conversion in pC
    double adccorr = adcpC*fabs(cos(theta[ilay]));
    
    if(adccorr<=0)INFOpmt[ihelp]=4;
         if(adccorr<=0.)           continue;
*/

adcl= adc[ihelp];
adcr= adc[ihelp+1];
adclraw = adcl;
adcrraw = adcr;
adcl   = f_adcPC( adcl );    // - adc conversion in pC
adcr   = f_adcPC( adcr );    // - adc conversion in pC

Float_t x;
if (ilay==0) x=yv[0];
if (ilay==1) x=xv[1];
if (ilay==2) x=xv[2];
if (ilay==3) x=yv[3];
if (ilay==4) x=yv[4];
if (ilay==5) x=xv[5];


 if (adclraw<4095) {
 adcl  = adcl*cos(theta[ilay]);
 xkorr = atten(A0[ihelp],A1[ihelp],A2[ihelp],A3[ihelp],x);
 xkorr = xkorr/hepratio;
 dedxl  = adcl/xkorr;
 eDEDXpmtraw[ihelp] = dedxl;
 Float_t dedxl1  = dedxl*4./dedx_corr[ihelp];  // 2nd order corr
// cout<<A0[ihelp]<<" "<<A1[ihelp]<<" "<<A2[ihelp]<<" "<<A3[ihelp]<<" => "<<xkorr<<endl;
// cout<<"Layer "<<ilay<<" PMT "<<ihelp+1<<" ADClraw "<<adclraw<<" adcl "<<adcl<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxl<<" dedx_2nd "<<dedxl1<<endl;  
 eDEDXpmt[ihelp]    = dedxl1;
// cout<<eDEDXpmt[ihelp]<<endl;
}

 if (adcrraw<4095)  {
 adcr = adcr*cos(theta[ilay]);
 xkorr = atten(A0[ihelp+1],A1[ihelp+1],A2[ihelp+1],A3[ihelp+1],x);
// cout<<A0[ihelp+1]<<" "<<A1[ihelp+1]<<" "<<A2[ihelp+1]<<" "<<A3[ihelp+1]<<" => "<<xkorr<<endl;
 xkorr = xkorr/hepratio;
//cout<<adcrraw<<" "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr"<<xkorr<<endl;
 dedxr  = adcr/xkorr;
 eDEDXpmtraw[ihelp+1] = dedxr;
 Float_t dedxr1  = dedxr*4./dedx_corr[ihelp+1];  // 2nd order corr
// cout<<"Layer "<<ilay<<" PMT "<<ihelp+1+1<<" ADCrraw "<<adcrraw<<" adcr "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxr<<" dedx_2nd "<<dedxr1<<endl;  
 eDEDXpmt[ihelp+1]  = dedxr1;
// cout<<eDEDXpmt[ihelp+1]<<endl;
                                   }


} // i1>-1
} // ilay

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


for (Int_t ilay=0; ilay<6; ilay++) {
Float_t  xhelp=1000.;
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
ipad  = ipad_a[ilay]+i1; 
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;                           
                        }                       
                                 } // 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++)INFOpmt[i]=1;
    for(int i=0;i<6;i++)INFOlayer[i]=1;
  }



 // define angle (simple way...  all theta[i] are the same...)   

Float_t theta[6] = {0., 0., 0., 0., 0., 0. };

if (fabs((tof->xtofpos[0])<100) && (fabs(tof->xtofpos[2])<100) && (fabs(tof->ytofpos[0])<100) && (fabs(tof->ytofpos[2])<100)) {
  Float_t dx   = tof->xtofpos[0] - tof->xtofpos[2];
  Float_t dy   = tof->ytofpos[0] - tof->ytofpos[2];
  Float_t dr   = sqrt(dx*dx+dy*dy);
  for (Int_t jj=0; jj<6; jj++) theta[jj]=atan(dr/77.31);
// cout<<tof->xtofpos[0]<<" "<<tof->xtofpos[2]<<" => "<<dx<<endl;
// cout<<tof->ytofpos[0]<<" "<<tof->ytofpos[2]<<" => "<<dy<<endl;
                                                               } // if...

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

  //  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++){
  eDEDXpmtrawstd[i] = 1000;
  eDEDXpmtstd[i] = 1000;
                       }

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


//  ratio helium to proton ca. 4
    Float_t  hepratio = 4.;

Int_t ipad,ihelp;
Float_t adclraw,adcrraw;
Float_t adcl,adcr;
Float_t dedxl,dedxr;
Float_t xkorr;

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++) {

Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {

ipad  = ipad_a[ilay]+i1; 
ihelp = ipmt_a[ilay]+(i1*2); 

/*    
    if(adc[ihelp]==4095)INFOpmt[ihelp]=2;
    else if(adc[iihelp]>=PMTsat[iihelp]-5)INFOpmt[iihelp]=3;
    
    if( adc[ihelp] >= PMTsat[iihelp]-5 )  continue;
    if( adc[ihelp] <= 0. )            continue;
    
    double adcpC   = f_adcPC( adc[iihelp] );    // - adc conversion in pC
    double adccorr = adcpC*fabs(cos(theta[ilay]));
    
    if(adccorr<=0)INFOpmt[ihelp]=4;
         if(adccorr<=0.)           continue;
*/

adcl= adc[ihelp];
adcr= adc[ihelp+1];
adclraw = adcl;
adcrraw = adcr;
adcl   = f_adcPC( adcl );    // - adc conversion in pC
adcr   = f_adcPC( adcr );    // - adc conversion in pC

Float_t x;
if (ilay==0) x = tof->ytofpos[0];
if (ilay==1) x = tof->xtofpos[0];
if (ilay==2) x = tof->xtofpos[1];
if (ilay==3) x = tof->ytofpos[1];
if (ilay==4) x = tof->ytofpos[2];
if (ilay==5) x = tof->xtofpos[2];

/*
Float_t xx;
if (ilay==0) xx=yv[0];
if (ilay==1) xx=xv[1];
if (ilay==2) xx=xv[2];
if (ilay==3) xx=yv[3];
if (ilay==4) xx=yv[4];
if (ilay==5) xx=xv[5];

cout<<ilay<<" xv "<<xx<<" tof std "<<x<<endl;
*/

if  (fabs(x) < 100.) {

 if (adclraw<4095) {
 adcl  = adcl*cos(theta[ilay]);
 xkorr = atten(A0[ihelp],A1[ihelp],A2[ihelp],A3[ihelp],x);
 xkorr = xkorr/hepratio;
 dedxl  = adcl/xkorr;
 eDEDXpmtrawstd[ihelp] = dedxl;
 Float_t dedxl1  = dedxl*4./dedx_corr[ihelp];  // 2nd order corr
// cout<<A0[ihelp]<<" "<<A1[ihelp]<<" "<<A2[ihelp]<<" "<<A3[ihelp]<<" => "<<xkorr<<endl;
// cout<<"Layer "<<ilay<<" PMT "<<ihelp+1<<" ADClraw "<<adclraw<<" adcl "<<adcl<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxl<<" dedx_2nd "<<dedxl1<<endl;  
 eDEDXpmtstd[ihelp]    = dedxl1;
// cout<<"trk "<<eDEDXpmt[ihelp]<<" std "<<eDEDXpmtstd[ihelp]<<endl;
}

 if (adcrraw<4095)  {
 adcr = adcr*cos(theta[ilay]);
 xkorr = atten(A0[ihelp+1],A1[ihelp+1],A2[ihelp+1],A3[ihelp+1],x);
// cout<<A0[ihelp+1]<<" "<<A1[ihelp+1]<<" "<<A2[ihelp+1]<<" "<<A3[ihelp+1]<<" => "<<xkorr<<endl;
 xkorr = xkorr/hepratio;
//cout<<adcrraw<<" "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr"<<xkorr<<endl;
 dedxr  = adcr/xkorr;
 eDEDXpmtrawstd[ihelp+1] = dedxr;
 Float_t dedxr1  = dedxr*4./dedx_corr[ihelp+1];  // 2nd order corr
// cout<<"Layer "<<ilay<<" PMT "<<ihelp+1+1<<" ADCrraw "<<adcrraw<<" adcr "<<adcr<<" x "<<x<<" theta "<<theta[ilay]<<" xkorr "<<xkorr<<" dedx "<<dedxr<<" dedx_2nd "<<dedxr1<<endl;  
 eDEDXpmtstd[ihelp+1]  = dedxr1;
// cout<<eDEDXpmt[ihelp+1]<<endl;
// cout<<"trk "<<eDEDXpmt[ihelp+1]<<" std "<<eDEDXpmtstd[ihelp+1]<<endl;
                                   }
}  // tofpos < 100

} // i1>-1
} // ilay

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


for (Int_t ilay=0; ilay<6; ilay++) {
Float_t  xhelp=1000.;
Int_t i1 = paddleidoftrack[ilay];
if (i1 != -1) {
ipad  = ipad_a[ilay]+i1; 
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;                           
                        }                       
                                 } // ilay

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


};

//------------------------------------------------------------------------
void ToFPatch::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 ToFPatch::ReadParAtt(const char *pardir, const char *param)
{
  char filename[100],filename1[100],filename2[100],name[50];
  char str1[10];



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

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

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

  //cout<<"read Time Interval File "<<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,"tofcalib_%s_atten.dat",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];
  

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

  
  fclose(fattin);


  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,"tofcalib_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];
  

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

  
  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];

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

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


//-----------------------------------------------------------
// Here I  read the dEdx_korr parameters 2nd order
//-----------------------------------------------------------

//Double_t t1,t2,tm;
UInt_t t1,t2,tm;
Float_t xmean1,xwidth1;
Float_t pmthelp[48];

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

sprintf(filename2,"%sdedx_2nd_corr_paddlehit.dat",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>>idummy>>xmean1>>xwidth1;
pmthelp[ii] = xmean1;
                       }
dedx_corr_m[jj].Set(48,pmthelp);
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,"%sdedx_2nd_corr_simu.dat",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>>idummy>>xmean1>>xwidth1;
pmthelp[ii] = xmean1;
                       }
dedx_corr_m[jj].Set(48,pmthelp);
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;


}


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

double ToFPatch::atten(float C0, float C1, float C2, float C3,  float x){
//
   Float_t yval = C0*exp(x*C1)+C2*exp(x*C3) ;
//   cout<<" in fl "<<C0<<" " <<C1<<" "<<C2<<" "<<C3<<" "<<C4<<" "<<x<<" "<<yval<<endl;
   return yval ;
}


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

double ToFPatch::f_adcPC( float x )
{
//  return 28.12+0.6312*x-5.647e-05*x*x+3.064e-08*x*x*x;
  return 28.0407 + 0.628929*x  - 5.80901e-05*x*x + 3.14092e-08*x*x*x;  // standard calibration curve since jan-07

}

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