//////////////////////////////////////
// TCHICsIEnergy implementation
// for all derived classes
//////////////////////////////////////

#include "TCHICsIEnergy.h"
#define tollup   1.03
#define tolldown 0.97

/////////////////////////////////////////////
// CsI Calibration for AsyEos experiment 
/////////////////////////////////////////////

//Constructor 
TCHICsIGSIEnergy::TCHICsIGSIEnergy(string filener, string filenerec, string calibdir) 
{
 ffilenerg = filener;
 ffilenergrec = filenerec;
 ffilecut=calibdir + "CHIcutsgood2.root";
 
 for(int i=0; i<MAXTEL; i++) {
  fa[i] = -1;
  fb[i] = -1;
 }

}

// parameters initialization for all telescopes
int TCHICsIGSIEnergy::Init()
{
 const int LEN=128;
 int row = 0;
 char line[LEN];
 int numtel;
 double a,b;
 ifstream freg(ffilenerg.c_str());
 if(!freg) {
  cout<<"TCHICsIGSIEnergy>> Calibration file "<<ffilenerg<<
        " can not be opened"<<endl;
  return -1;
 }
 
 istringstream is;
 while(freg.getline(line,LEN)) {
  line[strlen(line)]=' ';
   row++;
   switch(line[0]) {
   case '*' :
    break;
   default:
    is.str(line);
    is>>numtel>>a>>b;
    if(is.fail()) {
     cout<<"TCHICsIGSIEnergy>> Bad line "<<row<<" skipped"<<endl;
     is.clear();
     break;
    }
    if(numtel>=0 && numtel<MAXTEL) {
     fa[numtel] = a;
     fb[numtel] = b;
    }
   break;
  }
 }
 cout<<"TCHICsIGSIEnergy>> Read "<<row<<" lines from calibration file "
     <<ffilenerg<<endl;
 freg.close();
 fisdefined = true;
 //getchar();
//paolo 21/11/2011 
//  ifstream pippo("../CALFILES/de-e-bs.txt"); 
  ifstream pippo(ffilenergrec.c_str()); 
  if(!pippo) cout << "error: file " << ffilenergrec << " not found"<< endl;
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<31;iii++){
   pippo >> DE11[iii] >> E11[iii];
//   cout << DE11[iii] <<" "<< E11[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<45;iii++){
   pippo >> DE12[iii] >> E12[iii];
//   cout << DE12[iii] <<" "<< E12[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<41;iii++){
   pippo >> DE13[iii] >> E13[iii];
//   cout << DE13[iii] <<" "<< E13[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<45;iii++){
   pippo >> DE24[iii] >> E24[iii];
//   cout << DE24[iii] <<" "<< E24[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<41;iii++){
   pippo >> DE37[iii] >> E37[iii];
//   cout << DE37[iii] <<" "<< E37[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
 //getchar();
  for (int iii=0;iii<37;iii++){
   pippo >> DE49[iii] >> E49[iii];
//   cout << DE49[iii] <<" "<< E49[iii]<<endl;
  }
  pippo >> ZZZ >> AAA;
  cout <<"reading energy loss tabke for Z="<<  ZZZ <<" A="<< AAA<< endl;
  //getchar();
  for (int iii=0;iii<12;iii++){
   pippo >> DE511[iii] >> E511[iii];
//   cout << DE511[iii] <<" "<< E511[iii]<<endl;
  }
  pippo.close();

  TFile *fcut=new TFile(ffilecut.c_str());
  cut3Hebis=  (TCutG *)fcut->Get("bad3HeStop_bis");  
  cut3He=     (TCutG *)fcut->Get("bad3HeStop");  
  cut7Li =    (TCutG *)fcut->Get("bad7LiStop");  
  cut4He=     (TCutG *)fcut->Get("bad4HeStop");  
  fcut->Close();  
//paolo 21/11/2011 
 return 0;
}

// List data for debug 
void TCHICsIGSIEnergy::PrintData()
{
 if(!fisdefined)return;
 for(int i=0; i<MAXTEL; i++) {
  if(fa[i]!=-1 && fb[i]!= -1) {
   cout<<"Telescope: "<<setw(4)<<i<<" "<<fa[i]<<" "<<fb[i]<<endl;
  }
 }
}

// Energy calibration routine. 
// Da completare per Paolo con la formula per il calcolo dell'energia 
// totale per le particelle che traversano
void TCHICsIGSIEnergy::EvalEnergy(int numtel, int fast, int slow, TCHIResult *idr )
{
 double erec;
 int stopped=1;
 fDE=-1.; 
 fEnergy=-1.;
 fecod = kENER2;
 idr->SetStopped();
 fZ = idr->GetZ();             //charge
 fA = idr->GetA();             //mass
 fPID=idr->GetPID();            // PID 
 ECode cod = idr->Getcod();    //identification code
                  
 double de = -1.0;
 if(fa[numtel]!=-1 && fb[numtel]!=-1) {
  if(fast>0) {
   de = fast*fa[numtel] + fb[numtel];
  }
 }
//paolo 21/11/2011
int iii;
int myflag=0;
 if(de>0){
 if(fZ==2 && fA==2 && de<330){
  idr->SetA(0);fA=0;
  idr->SetZ(1);fZ=1;
  idr->Setcod(kICODE4);
 }  
 if(fZ==1){
  if(fA==1){
   if(de>222.1*tollup){
    idr->Setcod(kICODE3);
    idr->SetA(2);fA=2;
   }
   if(de>292.8*tollup){
    idr->Setcod(kICODE3);
    idr->SetA(3);fA=3;
   } 
   if(de>347.6*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(4);fA=4;
    idr->SetZ(2);fZ=2;
   }
   if(fA==1 && (fPID<0.5 || fPID>1.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }
  if(fA==2){
   if(de>292.8*tollup){
    idr->Setcod(kICODE3);
    idr->SetA(3);fA=3;
   } 
   if(de>347.6*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(4);fA=4;
    idr->SetZ(2);fZ=2;
   }  
   if(fA==2 && (fPID<1.5 || fPID>2.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }  
  if(fA==3){
   if(de>347.6*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(4);fA=4;
    idr->SetZ(2);fZ=2;
   }
   if(fA==3 && (fPID<2.5 || fPID>3.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
   if(fA==3 && (de>292.8*tollup) && (fPID!=3) && (idr->Getcod()<3))idr->Setcod(kICODE4);
   if(fA==3 && (de>292.8*tollup) && (fPID==3) && (idr->Getcod()<3))idr->Setcod(kICODE2);
  }  
  if(fA==0){
   if(de < 222.1 && de>0){
    for (iii=0;iii<31;iii++){
     if(de>= DE11[iii] && de < DE11[iii+1]){
      erec=E11[iii]+(de-DE11[iii])/(DE11[iii+1]-DE11[iii])*(E11[iii+1]-E11[iii]);
     }
     if(de<DE11[0]){erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE11[30])erec=de+1;
    }
    idr->SetA(1);fA=1;
    stopped=0; 
    idr->SetUnStopped();
    if(idr->Getcod()<=3)idr->Setcod(kICODE3);
   }
   if(de>222.1 && de<292.8){
    for (iii=0;iii<45;iii++){
     if(de>= DE12[iii] && de < DE12[iii+1]){;
      erec=E12[iii]+(de-DE12[iii])/(DE12[iii+1]-DE12[iii])*(E12[iii+1]-E12[iii]);
     }
     if(de<DE12[0]){
      erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE12[44])erec=de+1;
    }
    idr->SetA(2);fA=2;
    stopped=0;
    idr->SetUnStopped();
    if(idr->Getcod()<=3)idr->Setcod(kICODE3);
   }  
   if(de>292.8 && de<347.6*tollup){
    for (iii=0;iii<41;iii++){
     if(de>= DE13[iii] && de < DE13[iii+1]){;
      erec=E13[iii]+(de-DE13[iii])/(DE13[iii+1]-DE13[iii])*(E13[iii+1]-E13[iii]);
     }
     if(de<DE13[0]){
      erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE13[40])erec=de+1;
    }
    idr->SetA(3);fA=3;
    stopped=0;
    idr->SetUnStopped();
    if(idr->Getcod()<=3)idr->Setcod(kICODE3);
   }  
   if(de>347.6*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(2);fA=2;
    idr->SetZ(2);fZ=2;
   }  
  }
 }
 if(fZ==2){
  myflag==0;
  if(fA==3){
   if(de>792.4*tollup){
    idr->Setcod(kICODE3);
    idr->SetA(4);fA=4;
   }
   if(cut3Hebis->IsInside(fPID,de)){
//    idr->Setcod(kICODE4);
    idr->SetA(2);fA=2;
    myflag=1;
   }
   if(fA==3 && (fPID<2.5 || fPID>3.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }
  if(fA==4){
   if(de>888.1*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(7);fA=7;
    idr->SetZ(3);fA=3;
   }
   if(cut4He->IsInside(fPID,de))idr->Setcod(kICODE5);
   if(fA==4 && (fPID<3.5 || fPID>4.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }
  if(fA==2){
   if(de>888.1*tollup && myflag==0){
    idr->Setcod(kICODE4);
    idr->SetA(3);fA=3;
    idr->SetZ(3);fZ=3;
    stopped=0;
    idr->SetUnStopped();
   }  
   if(de<888.1*tollup){
    for (iii=0;iii<45;iii++){
     if(de>= DE24[iii] && de < DE24[iii+1]){;
      erec=E24[iii]+(de-DE24[iii])/(DE24[iii+1]-DE24[iii])*(E24[iii+1]-E24[iii]);
     }
     if(de<DE24[0]){
      erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE24[44])erec=de+1;
    }
    idr->SetA(4);fA=4;
    stopped=0;
    idr->SetUnStopped();
    if(idr->Getcod()!=5)idr->Setcod(kICODE3);
    if(myflag==1)idr->Setcod(kICODE5);
   }
  }
 }
 if(fZ==3){
  if(fA==7){
   if(de>1808.1*tollup){
    idr->Setcod(kICODE4);
    idr->SetZ(4);fA=4;
    idr->SetA(9);fA=9;
   } 
   if(cut7Li->IsInside(fPID,de))idr->Setcod(kICODE5);
   if(fA==7 && (fPID<5.5 || fPID>7.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }
  if(fA==3){
   if(de>1808.1*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(4);fA=4;
    idr->SetZ(4);fZ=4;
   }  
   if(de < 1808.1*tollup){
    for (iii=0;iii<41;iii++){
     if(de>= DE37[iii] && de < DE37[iii+1]){;
      erec=E37[iii]+(de-DE37[iii])/(DE37[iii+1]-DE37[iii])*(E37[iii+1]-E37[iii]);
     }
     if(de<DE37[0]){
      erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE37[40])erec=de+1;
    }
    idr->SetA(7);fA=7;
    stopped=0;
    idr->SetUnStopped();
    if(idr->Getcod()!=5)idr->Setcod(kICODE3);
   }  
  }  
 }  
 if(fZ==4){
  if(fA==9){
   if(de>2836.8*tollup){
    idr->Setcod(kICODE4);
    idr->SetA(11);fA=11;
    idr->SetZ(5);fZ=5;
   }  
   if(fA==9 && (fPID<6.5 || fPID>10.5) && (idr->Getcod()<3))idr->Setcod(kICODE5);
  }
  if(fA==4){
   if(de>2863.9*tollup){
   idr->Setcod(kICODE4);
   idr->SetA(5);fA=5;
   idr->SetZ(5);fZ=5;
  }
  if(de<2863.9*tollup){
   for (iii=0;iii<37;iii++){
     if(de>= DE49[iii] && de < DE49[iii+1]){;
      erec=E49[iii]+(de-DE49[iii])/(DE49[iii+1]-DE49[iii])*(E49[iii+1]-E49[iii]);
     }
     if(de<DE49[0]){
      erec=-100;
      idr->Setcod(kICODE5);
     }
     if(de>DE49[36])erec=de+1;
    }
    idr->SetA(9);fA=9;
    stopped=0;
    idr->SetUnStopped();
    if(idr->Getcod()!=5)idr->Setcod(kICODE3);
   }  
  }  
 }  
 if(fZ==5 && fA==5){
  for (iii=0;iii<12;iii++){
   if(de>= DE511[iii] && de < DE511[iii+1]){;
    erec=E511[iii]+(de-DE511[iii])/(DE511[iii+1]-DE511[iii])*(E511[iii+1]-E511[iii]);
   }
   if(de<DE511[0]){
    erec=-100;
    idr->Setcod(kICODE5);
   }
   if(de>DE511[11])erec=de+1;
  }
  idr->SetA(11);fA=11;
  stopped=0;
  idr->SetUnStopped();
  if(idr->Getcod()!=5)idr->Setcod(kICODE3);
 }  
 if((idr->Getcod()==8) && (idr->GetZ()==-1) && (slow>2500)){
  if(de<888.1){
   idr->SetA(4);fA=4;
   idr->SetZ(2);fZ=2;
  }
  if(de> 888.1 && de<1808.1){
   idr->SetA(7);fA=7;
   idr->SetZ(3);fZ=3;
  }
  if(de> 1808.1 && de<2863.9){
   idr->SetA(9);fA=9;
   idr->SetZ(4);fZ=4;
  }
  if(de> 2863.9){
   idr->SetA(11);fA=11;
   idr->SetZ(5);fZ=5;
  }
  idr->Setcod(kICODE6);
  idr->SetStopped();
  stopped==1;
 }

 if(fZ==2 && fA==2 && stopped==1 && (idr->Getcod()<8)){
  //getchar();
  cout <<" 22 " << de << endl;
 } 
 if(fZ==1 && fA==0 && stopped==1 && (idr->Getcod()<8)&& de>0){
  //getchar();
  cout <<" 10 " << de << endl;
 } 
 if(fZ==3 && fA==3 && stopped==1 && (idr->Getcod()<8)){
  //getchar();
  cout <<" 33 " << de << endl;
 } 
 if(fZ==4 && fA==4 && stopped==1 && (idr->Getcod()<8)){
  //getchar();
  cout <<" 44 " << de << endl;
 } 
 if(fZ==5 && fA==5 && stopped==1 && (idr->Getcod()<8)){
  //getchar();
  cout <<" 55 " << de << endl;
 } 
 }
 if(stopped)fDE = de;
 if(!stopped)fDE = erec;
 fecod =  kENER0;
}