#include <iostream>
#include <libgen.h>
#include <math.h>

#include "TChain.h"
#include "TFile.h"
#include "TRandom3.h"
#include "TObjArray.h"
#include <Compression.h>

#include "anyoption.h"

#include "ECALboxcox.h"

using namespace std;

bool VERBOSE;
bool QUICK;
bool DEBUG;
bool PRINTENE;

int ENERGY;
int SPAN;

int NENE;

unsigned int Run;
unsigned int Event;
unsigned int Time;
unsigned int RunTag; 

float trkQ;
float trkInnQ;
float tofQ;
float trkDefRig;
float ecalEnergyE;
float ecalCorrEneE_E;
float ecalCorrEneP_E;
float ecalGetCorrEne2015;
float ecalGetCorrEne2017;
bool goodTrdKCalib;
float trkDefChiY;
float trkDefChiX;
float partCutOff;
int passRTI; 

short nTrTracks;    
float trkFullRig;
float trkInnRig;
float trkHalfUpRig;
float trkHalfDwnRig;
float ecalCoG[3];
float ecalTrkCoG[3];
float ecalLayerFrac[18];
short trdKHitsOnTrkTrack;
short trdKHitsOnTrdTrack;
short trkSpan;
short nACC;
float trdKlhrEP_ene;
float trdKlhE_ene;
float trdKlhE_ene_trdstd;
float trdKlhrEH_ene;
float ecalBDT;
float ecalBDTchi2;
float ecalEnergyD;
short runAnalysisTag;
int _runAnalysisTag;

float rtiNtrig;
float rtiNev;
float rtiNerr;
float rtiNpart;
float rtiLifetime;
float rtiZenith;
int rtiGood;

float accTheta[21];
float accPhi[21];
float accP[21][3];

float eop;
float trdclass;
float trdclass_trdstd;
float thetaontop;
float phiontop;
float bcnlay[6];
float bcnlay_ns[6];
float bcnslay[5];

float scmean,    scmean_5lay;
float scdummylh, scdummylh_5lay;
float sclhv0,    sclhv5,   sclhv6;
float sclhdv0,   sclhdv5,  sclhdv6;
float scbdtv0,   scbdtv5,  scbdtv6;
float scbdtgv0,  scbdtgv5, scbdtgv6;

float ecalEneEP;

//Additional variables for ECAL BDT training
// float ecalS1tot[3];
// float ecalS3tot[3];
// float ecalRadiusFrac[3];
// float ecalLayerS13[18];
// float ecalShowerHits; //this is float, should be int
// short ecalTotalHits;

Int_t ecal_k_hadflag, ecal_k_N_Shwr;
Float_t ecal_k_EleEne, ecal_k_EmBDT, ecal_k_EmLkhd[2];

Int_t ecal_k_ShwrStat;
Float_t ecal_k_ShwrE0, ecal_k_ShwrA0, ecal_k_ShwrX0,
  ecal_k_ShwrY0, ecal_k_ShwrZ0, ecal_k_ShwrKX, ecal_k_ShwrKY;

//#define ALSOPRESCALED
#ifdef ALSOPRESCALED
#define NSIGN 3
const char *sig_name[NSIGN] = { "Negative", "Positive", "Prescale" };
#else
#define NSIGN 2
const char *sig_name[NSIGN] = { "Negative", "Positive" };
#endif
#define NTAG 2
const char *tag_name[NTAG] = { "ForAnal", "NotForAnal" };

string OUTPUTDIR;
string FILELIST;
int ENECLASS;
int NTUPLEVERSION;
string DIRFRIEND;
bool ADDFRIEND;

bool PassRTI(){

  if( rtiGood!=0 ) return false; //{ cout<<"Fail 0"<<endl; }
  if( !(rtiNtrig/rtiNev>0.98) ) return false;//{ cout<<"Fail 1"<<endl; }
  if( !(rtiNpart/rtiNtrig>0.07/1600*rtiNtrig&&rtiNpart/rtiNtrig<0.25) ) return false;// { cout<<"Fail 2"<<endl; }
  if( !(rtiLifetime>0.5) ) return false;//{ cout<<"Fail 3"<<endl; }
  if( !(rtiZenith<40) ) return false;//{ cout<<"Fail 4"<<endl; }
  if( !(rtiNerr>=0&&rtiNerr/rtiNev<0.1) ) return false;//{ cout<<"Fail 5"<<endl; }
  if( !(rtiNpart>0&&rtiNev<1800) ) return false;// cout<<"Fail 6"<<endl; }
  return true;

}


bool matchX;
bool matchY;

void TrkEcalMatch( float ene, float rig, float ecalX, float trkX, float ecalY, float trkY){

  float Dx = (trkX-ecalX)/(sqrt(pow(-0.05118,2)+pow((0.3724/ene),2))+pow(10.,-0.6946-0.407*log10(ene)+0.1675*pow(log10(ene),2)-0.254*pow(log10(ene),3)));
  float Dy = (trkY-ecalY)/(sqrt(pow(0.05618,2)+pow(0.5049/ene,2))+pow(10,-1.223+0.01248*log10(ene)+0.4523*pow(log10(ene),2)-0.5666*pow(log10(ene),3)));
  
  matchX = false;
  if(Dx>=-0.21-3*0.93 && Dx<=-0.21+3*0.93) matchX=true;
  
  matchY = false;
  if( rig<0 )
    { if(Dy>=(0.5821-3*1.028) && Dy<=(0.5821+3*2.333)) matchY=true; }
  if( rig>0 )
    { if(Dy<=(0.5821+3*1.028) && Dy>=(0.5821-3*2.333)) matchY=true; }

  return;

}

Double_t *MinEne;
Double_t *MaxEne;

int GetEnergyIndex( float ene ){
  
  ene = TMath::Abs( ene );
  for(int iene=0; iene<NENE; iene++)
    {
      if( ene>MinEne[iene] && ene<=MaxEne[iene] ) return iene;
    }
  if( ene > MaxEne[NENE-1] ) return 999;
  return -1;
  
}


char listname[1024];           //Filelist to be processed   
int list_first_line;           //First line in the list to be processed                                                                      
int list_last_line;            //Last line in the list to be processed     

//Was used as a "patch" for pass4 data. Should be unused now
float EneCorrT( TSpline3 *spcorr, const float ene, const unsigned int Time ){
  float corr = spcorr->Eval(Time);
  corr = TMath::Abs(corr-1.)>0.05 ? 1.0055 : 1.0055/corr;
  return ene*corr;
}

int Build_superskimmed(){

  // string fileecalcorr;
//   if( !getenv("AMSACOMMONSW") ) { printf("Cannot find AMSACOMMONSW env var\n"); return 76; }
//   else { fileecalcorr.assign( Form( "%s/TreeBuilder/include/Result_Kalman.root", getenv("AMSACOMMONSW") ) ); 
//     cout<<"Reading energy time-dependent corrections from "<<fileecalcorr.c_str()<<endl; }
  
  TChain *chain = new TChain("tree");
  TChain *friendchain = NULL;
  if (ADDFRIEND) {
    friendchain = new TChain("tree");
  }

  sprintf( listname, "%s", FILELIST.c_str() );
  printf("Adding files from list %s\n", listname);
  printf("First line to add: %d\n",list_first_line);
  printf("Last  line to add: %d\n",list_last_line);
  FILE* listfile = fopen(listname,"r");
  char file_to_add[1024];
  int ifile=0;
  ofstream  rejfile;
  while (!feof(listfile)){
    if( ifile>list_last_line ) break;
    fscanf(listfile,"%s\n", file_to_add );
    if(ifile>=list_first_line && ifile<=list_last_line){
      string sname; sname.assign(file_to_add);
      string bname; bname.assign(basename(file_to_add));
      size_t root_pos=sname.find(".root");
      bool addthis = true;
      if(root_pos == string::npos){
	cout<<Form("%s is not a rootfile. Skipping this",file_to_add )<<endl;
	cerr<<Form("%s is not a rootfile. Skipping this",file_to_add )<<endl;
	addthis=false;
      }
      else{
	chain->Add(file_to_add);
	//	printf("Adding %s\n", file_to_add);
	if (ADDFRIEND) {
	  string dirfriend = DIRFRIEND;
	  string fullfriendname = dirfriend+"/"+bname;
	  //	  printf("Adding %s\n", fullfriendname.c_str());
	  friendchain->Add(fullfriendname.c_str());
	}
      }
    }
    ifile++;
  }

  if (ADDFRIEND){
    chain->AddFriend(friendchain);
  }
  
  cout<<Form("\n **** Files in the chain ****\n");
  chain->GetListOfFiles()->Print();
  cout<<Form("\n");
  cout<<Form("%d Entries in the chain\n\n",(int)chain->GetEntries());

  printf("Building the index...\n\n");
  chain->BuildIndex("Run", "Event");
  if (ADDFRIEND) {
    friendchain->BuildIndex("Run", "Event");
  }
  chain->BuildIndex("Run", "Event");
  
  chain->SetBranchAddress("Run",                  &Run);
  chain->SetBranchAddress("RunTag",               &RunTag);
  chain->SetBranchAddress("Event",                &Event);
  chain->SetBranchAddress("Time",                 &Time);
  chain->SetBranchAddress("trkInnQ",              &trkInnQ);
  chain->SetBranchAddress("ecalEnergyE",          &ecalEnergyE);
  chain->SetBranchAddress("trkDefRig",            &trkDefRig);
  chain->SetBranchAddress("ecalCorrEneE_E",       &ecalCorrEneE_E);
  chain->SetBranchAddress("ecalCorrEneP_E",       &ecalCorrEneP_E);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchAddress("ecalGetCorrEne2015",   &ecalGetCorrEne2015);
  }
  if (NTUPLEVERSION>=2) {
    chain->SetBranchAddress("ecalGetCorrEne2017",   &ecalGetCorrEne2017);
  }
  if (NTUPLEVERSION==0) {
    chain->SetBranchAddress("trdKLikeEP_ene",        &trdKlhrEP_ene );
    chain->SetBranchAddress("trdKLikeEH_ene",        &trdKlhrEH_ene );
    chain->SetBranchAddress("trdKLikeEP_ene_prob",   &trdKlhE_ene );
  }
  else {
    chain->SetBranchAddress("trdKlhrEP_ene",        &trdKlhrEP_ene );
    chain->SetBranchAddress("trdKlhrEH_ene",        &trdKlhrEH_ene );
    chain->SetBranchAddress("trdKlhE_ene",          &trdKlhE_ene );
    chain->SetBranchAddress("trdKlhE_ene_trdstd",   &trdKlhE_ene_trdstd );
  }
  
  chain->SetBranchAddress("passRTI",              &passRTI);
  chain->SetBranchAddress("nTrTracks",            &nTrTracks);
  chain->SetBranchAddress("goodTrdKCalib",        &goodTrdKCalib);
  chain->SetBranchAddress("trdKHitsOnTrkTrack",   &trdKHitsOnTrkTrack);
  if (NTUPLEVERSION>=1){
    chain->SetBranchAddress("trdKHitsOnTrdTrack",   &trdKHitsOnTrdTrack);
  }
  chain->SetBranchAddress("ecalEnergyD",          &ecalEnergyD );
  chain->SetBranchAddress("ecalLayerFrac",        &ecalLayerFrac );
  chain->SetBranchAddress("trkSpan",              &trkSpan);
  if (NTUPLEVERSION==0) {
    chain->SetBranchAddress("runAnalysisTag",       &_runAnalysisTag);
  }
  else {
    chain->SetBranchAddress("runAnalysisTag",       &runAnalysisTag);
  }
  chain->SetBranchAddress("rtiLifetime",          &rtiLifetime);

  chain->SetBranchAddress("accTheta",             &accTheta[0]);
  chain->SetBranchAddress("accPhi",               &accPhi[0]);
  chain->SetBranchAddress("accP",                  accP);
  
  if (NTUPLEVERSION>=2) {
    chain->SetBranchAddress("ecalKHadFlag",  &ecal_k_hadflag );
    chain->SetBranchAddress("ecalKNShwr",    &ecal_k_N_Shwr  );
    chain->SetBranchAddress("ecalKEleEne",   &ecal_k_EleEne  );
    chain->SetBranchAddress("ecalKEmBDT",    &ecal_k_EmBDT   );
    chain->SetBranchAddress("ecalKEmLkhd",    ecal_k_EmLkhd  );
    chain->SetBranchAddress("ecalKShwrStat", &ecal_k_ShwrStat);
    chain->SetBranchAddress("ecalKShwrE0",   &ecal_k_ShwrE0  );
    chain->SetBranchAddress("ecalKShwrA0",   &ecal_k_ShwrA0  );
    chain->SetBranchAddress("ecalKShwrX0",   &ecal_k_ShwrX0  );
    chain->SetBranchAddress("ecalKShwrY0",   &ecal_k_ShwrY0  );
    chain->SetBranchAddress("ecalKShwrZ0",   &ecal_k_ShwrZ0  );
    chain->SetBranchAddress("ecalKShwrKX",   &ecal_k_ShwrKX  );
    chain->SetBranchAddress("ecalKShwrKY",   &ecal_k_ShwrKY  );
  }

  chain->SetBranchStatus("*",0);
  
  //****
  //Re-activate old branches
  //****
  chain->SetBranchStatus("runAnalysisTag",             1);
  chain->SetBranchStatus("Run",                        1);
  chain->SetBranchStatus("RunTag",                     1);
  chain->SetBranchStatus("whichBadRun",                1);
  chain->SetBranchStatus("Event",                      1);
  chain->SetBranchStatus("Time",                       1);
  chain->SetBranchStatus("IssTheta",                   1);
  chain->SetBranchStatus("IssPhi",                     1);
  chain->SetBranchStatus("magTheta",                   1);
  chain->SetBranchStatus("magPhi",                     1);
  chain->SetBranchStatus("nACC",                       1);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("IssInSAA",                   1);
  }
  
  chain->SetBranchStatus("trkQ",                       1);
  chain->SetBranchStatus("trkInnQ",                    1);
  chain->SetBranchStatus("trkQlay",                    1);
  if (NTUPLEVERSION>=2) {
    chain->SetBranchStatus("trkQStatus",                 1);
    chain->SetBranchStatus("trkMaxQL1",                  1);
    chain->SetBranchStatus("trkMaxQL1Status",            1);
  }
  
  chain->SetBranchStatus("tofQ",                       1);
  chain->SetBranchStatus("tofQlay",                    1);
  // chain->SetBranchStatus("tofLayerEdep",               1);
  chain->SetBranchStatus("tofUpperQ",                  1);
  chain->SetBranchStatus("tofLowerQ",                  1);
  chain->SetBranchStatus("tofBeta",                    1);
  
  chain->SetBranchStatus("trigJMembPatt",              1);
  chain->SetBranchStatus("trigPhysBPatt",              1);
  chain->SetBranchStatus("trigPhysBPatt_notrestored",   1);
  
  chain->SetBranchStatus("nTrdTracks",                 1);
  chain->SetBranchStatus("trdTotalHits",               1);
  chain->SetBranchStatus("trdTrackHits",               1);
  chain->SetBranchStatus("goodTrdKCalib",              1);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("trdKlhrEP",                 1);
    chain->SetBranchStatus("trdKlhrHP",                 1);
    chain->SetBranchStatus("trdKlhrEH",                 1);
    chain->SetBranchStatus("trdKlhrEP_ene",             1);
    chain->SetBranchStatus("trdKlhrHP_ene",             1);
    chain->SetBranchStatus("trdKlhrEH_ene",             1);
    chain->SetBranchStatus("trdKlhE_ene",        1);
    chain->SetBranchStatus("trdKlhP_ene",        1);
    chain->SetBranchStatus("trdKlhH_ene",        1);
    chain->SetBranchStatus("trdKlhrEP_ene_trdstd",             1);
    chain->SetBranchStatus("trdKlhrHP_ene_trdstd",             1);
    chain->SetBranchStatus("trdKlhrEH_ene_trdstd",             1);
    chain->SetBranchStatus("trdKlhE_ene_trdstd",        1);
    chain->SetBranchStatus("trdKlhP_ene_trdstd",        1);
    chain->SetBranchStatus("trdKlhH_ene_trdstd",        1);
  }
  else {
    chain->SetBranchStatus("trdKLikeEP",                 1);
    chain->SetBranchStatus("trdKLikeHP",                 1);
    chain->SetBranchStatus("trdKLikeEH",                 1);
    chain->SetBranchStatus("trdKLikeEP_ene",             1);
    chain->SetBranchStatus("trdKLikeHP_ene",             1);
    chain->SetBranchStatus("trdKLikeEH_ene",             1);
    chain->SetBranchStatus("trdKLikeEP_ene_prob",        1);
    chain->SetBranchStatus("trdKLikeHP_ene_prob",        1);
    chain->SetBranchStatus("trdKLikeEH_ene_prob",        1);
  }
  chain->SetBranchStatus("trkDefRig",                  1);
  chain->SetBranchStatus("trkInnRig",                  1);
  chain->SetBranchStatus("trkHalfDwnRig",              1);  
  chain->SetBranchStatus("trkHalfUpRig",               1);  
  chain->SetBranchStatus("trkFullRig",                 1); 
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("trkEnhDefRig",                  1);
    chain->SetBranchStatus("trkEnhInnRig",                  1);
    chain->SetBranchStatus("trkEnhHalfDwnRig",              1);  
    chain->SetBranchStatus("trkEnhHalfUpRig",               1);  
    chain->SetBranchStatus("trkEnhFullRig",                 1); 
  }
  chain->SetBranchStatus("ecalEnergyE",                1);  
  chain->SetBranchStatus("ecalEnergyA",                1);  
  chain->SetBranchStatus("ecalCorrEneE_E",             1);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("ecalGetCorrEne2015",         1);
  }
  if (NTUPLEVERSION>=2) {
    chain->SetBranchStatus("ecalGetCorrEne2017",         1);
  }
  chain->SetBranchStatus("ecalCorrEneE_P",             1);
  chain->SetBranchStatus("ecalCorrEneP_E",             1);
  chain->SetBranchStatus("ecalCorrEneP_P",             1);
  chain->SetBranchStatus("rtiEvent",                   1);  
  chain->SetBranchStatus("rtiMaxCutoff",               1);
  chain->SetBranchStatus("rtiMaxIGRFCutoff",           1);
  chain->SetBranchStatus("rtiLifetime",                1);  
  chain->SetBranchStatus("rtiRad",                     1);  
  chain->SetBranchStatus("rtiNev",                     1);  
  chain->SetBranchStatus("rtiNerr",                    1);  
  chain->SetBranchStatus("rtiNtrig",                   1);  
  chain->SetBranchStatus("rtiNpart",                   1);  
  chain->SetBranchStatus("rtiGood",                    1);  
  chain->SetBranchStatus("rtiTheta",                   1);  
  chain->SetBranchStatus("rtiPhi",                     1);  
  chain->SetBranchStatus("rtiZenith",                  1);
  chain->SetBranchStatus("passRTI",                    1);
  //chain->SetBranchStatus("trkMinoise",               1);
  //  chain->SetBranchStatus("GeoCutOff_Pos",          1);
  //  chain->SetBranchStatus("GeoCutOff_Neg",          1);
  // chain->SetBranchStatus("GeoCutOff",               1);
  
  chain->SetBranchStatus("trdKHitsOnTrkTrack",         1);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("trdKHitsOnTrdTrack",         1);
  }

  chain->SetBranchStatus("trkSpan",                    1);
  chain->SetBranchStatus("trkHasLay2",                 1);
  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("trkEnhSpan",                 1);
    chain->SetBranchStatus("trkEnhHasLay2",              1);
  }

  chain->SetBranchStatus("nTrTracks",                  1);
  chain->SetBranchStatus("trkDefErrInvRig",            1);
  chain->SetBranchStatus("trkDefChiY",                 1);
  chain->SetBranchStatus("trkDefChiX",                 1);
  chain->SetBranchStatus("trkDefTheta",                1);
  chain->SetBranchStatus("trkDefPhi",                  1);
  chain->SetBranchStatus("trkMaxRes",                  1);
  chain->SetBranchStatus("trkMaxResLay",               1);
  // chain->SetBranchStatus("trkFixBadX",                 1);
  // chain->SetBranchStatus("trkSimpleChi2X",             1);
  

  chain->SetBranchStatus("nEcalShowers",               1);  
  chain->SetBranchStatus("ecalEnergyD",                1);  
  
 //  chain->SetBranchStatus("ecalBDTv5",                  1);  
//   chain->SetBranchStatus("ecalBDTv5s",                 1);  
//   chain->SetBranchStatus("ecalBDTv5E",                 1);  
  chain->SetBranchStatus("ecalBDTv5Es",                1);  
  chain->SetBranchStatus("ecalBDTv6Es",                1);  
  chain->SetBranchStatus("ecalBDTv7Es",                1);  
 //  chain->SetBranchStatus("ecalBDTv4",                  1);  
//   chain->SetBranchStatus("ecalBDTv104",                1);  
//   chain->SetBranchStatus("ecalBDTchi2v3",              1);  
//   chain->SetBranchStatus("ecalBDTchi2v3s",             1);  
//    chain->SetBranchStatus("ecalBDTchi2v2",              1);  
//   chain->SetBranchStatus("ecalBDTchi2v3E",             1);  
//   chain->SetBranchStatus("ecalBDTchi2v3Es",            1);  
//   chain->SetBranchStatus("ecalBDTchi2v102",            1);  
//   chain->SetBranchStatus("ecalLappEst",                1);  
//   chain->SetBranchStatus("ecalaxisChi2lf",             1);
//   chain->SetBranchStatus("ecalaxisChi2cm",             1);

  chain->SetBranchStatus("ecalLayerFrac",              1);

  chain->SetBranchStatus("ecalTrkEntry",               1);  
  chain->SetBranchStatus("ecalTrkExit",                1);  
  chain->SetBranchStatus("ecalCoG",                    1);
  chain->SetBranchStatus("ecalTrkCoG",                 1);
  chain->SetBranchStatus("ecalDir",                    1);
  chain->SetBranchStatus("ecalTrkDir",                 1);
  chain->SetBranchStatus("ecalStatus",                 1);

  chain->SetBranchStatus("ecalS1tot",                  1);
  chain->SetBranchStatus("ecalS3tot",                  1);
  chain->SetBranchStatus("ecalRadiusFrac",             1);
  chain->SetBranchStatus("ecalLayerS13",               1);
  chain->SetBranchStatus("ecalShowerHits",             1);
  chain->SetBranchStatus("ecalTotalHits",              1);
  
  chain->SetBranchStatus("nParticles",                 1);  
  chain->SetBranchStatus("partMomentum",               1);  
  chain->SetBranchStatus("partBeta",                   1);  

  chain->SetBranchStatus("accPhi",                     1);
  chain->SetBranchStatus("accTheta",                   1);
  chain->SetBranchStatus("accP",                       1);

  if (NTUPLEVERSION>=1) {
    chain->SetBranchStatus("trkccEval",                  1);
    chain->SetBranchStatus("trkccVar",                   1);
  }

  if (NTUPLEVERSION>=2) {
    chain->SetBranchStatus("ecalKHadFlag",               1);
    chain->SetBranchStatus("ecalKNShwr",                 1);
    chain->SetBranchStatus("ecalKEleEne",                1);
    chain->SetBranchStatus("ecalKEmBDT",                 1);
    chain->SetBranchStatus("ecalKEmLkhd",                1);
    chain->SetBranchStatus("ecalKShwrStat",              1);
    chain->SetBranchStatus("ecalKShwrE0",                1);
    chain->SetBranchStatus("ecalKShwrA0",                1);
    chain->SetBranchStatus("ecalKShwrX0",                1);
    chain->SetBranchStatus("ecalKShwrY0",                1);
    chain->SetBranchStatus("ecalKShwrZ0",                1);
    chain->SetBranchStatus("ecalKShwrKX",                1);
    chain->SetBranchStatus("ecalKShwrKY",                1);
  }
  
  TFile *f[NSIGN][NTAG][NENE+1];
  TTree *t[NSIGN][NTAG][NENE+1];

  char dir[1024];
  sprintf( dir, "%s", OUTPUTDIR.c_str() );

  printf("Creating the output files and trees...\n");
  
  int complevel=ROOT::CompressionSettings(ROOT::kLZMA, 2);
  printf("The choosen compression level is %d\n", complevel);

  for(int itag=0; itag<NTAG; itag++)
    {
      
      for(int isign=0; isign<NSIGN; isign++)
	for(int iene=0; iene<NENE+1; iene++)
	  {
	    char fname[1024];
	    sprintf( fname, "%s_%s%s_%.2d_%.6d_%.6d.root", sig_name[isign], tag_name[itag], (SPAN<=5&&SPAN>=0)?Form("_sp%d",SPAN):"", iene, list_first_line, list_last_line);
	    f[isign][itag][iene] = new TFile( Form("%s/%s", dir, fname), "RECREATE", "File with Superskimmed Hadded trees", complevel);
	    f[isign][itag][iene]->cd();
	    t[isign][itag][iene] = chain->CloneTree(0);
	    
	    // t[isign][itag][iene]->SetAutoSave(1000000);//1 MB
	    //	    t[isign][itag][iene]->SetAutoFlush(-1000000);//1 MB
	    t[isign][itag][iene]->SetAutoFlush(-1000000);//1 MB
	    //	    t[isign][itag][iene]->SetAutoSave(100000);//0.1 MB
	    //	    t[isign][itag][iene]->SetAutoFlush(-100000);//0.1 MB
	    t[isign][itag][iene]->SetDirectory(f[isign][itag][iene]);

	    //****
	    //New Branches
	    //****
	    t[isign][itag][iene]->Branch("eop",           &eop,           "eop/F");
	    t[isign][itag][iene]->Branch("trdclass",      &trdclass,      "trdclass/F");
	    t[isign][itag][iene]->Branch("trdclass_trdstd",      &trdclass_trdstd,      "trdclass_trdstd/F");
	    t[isign][itag][iene]->Branch("thetaontop",    &thetaontop,    "thetaontop/F");
	    t[isign][itag][iene]->Branch("phiontop",      &phiontop,      "phiontop/F");
	    //
	    t[isign][itag][iene]->Branch("bcnlay",        bcnlay,         "bcnlay[6]/F");
	    t[isign][itag][iene]->Branch("bcnlay_ns",     bcnlay_ns,      "bcnlay_ns[6]/F");
	    //
	    t[isign][itag][iene]->Branch("scmean",        &scmean,        "scmean/F");
	    t[isign][itag][iene]->Branch("scmean_5lay",   &scmean_5lay,   "scmean_5lay/F");
	    t[isign][itag][iene]->Branch("scdummylh",     &scdummylh,     "scdummylh/F");
	    t[isign][itag][iene]->Branch("scdummylh_5lay",&scdummylh_5lay,"scdummylh_5lay/F");
	    t[isign][itag][iene]->Branch("sclhv0",        &sclhv0,        "sclhv0/F");
	    t[isign][itag][iene]->Branch("sclhv5",        &sclhv5,        "sclhv5/F");
	    t[isign][itag][iene]->Branch("sclhv6",        &sclhv6,        "sclhv6/F");
	    t[isign][itag][iene]->Branch("sclhdv0",       &sclhdv0,       "sclhdv0/F");
	    t[isign][itag][iene]->Branch("sclhdv5",       &sclhdv5,       "sclhdv5/F");
	    t[isign][itag][iene]->Branch("sclhdv6",       &sclhdv6,       "sclhdv6/F");
	    t[isign][itag][iene]->Branch("scbdtv0",       &scbdtv0,       "scbdtv0/F");
	    t[isign][itag][iene]->Branch("scbdtv5",       &scbdtv5,       "scbdtv5/F");
	    t[isign][itag][iene]->Branch("scbdtv6",       &scbdtv6,       "scbdtv6/F");
	    t[isign][itag][iene]->Branch("scbdtgv0",      &scbdtgv0,      "scbdtgv0/F");
	    t[isign][itag][iene]->Branch("scbdtgv5",      &scbdtgv5,      "scbdtgv5/F");
	    t[isign][itag][iene]->Branch("scbdtgv6",      &scbdtgv6,      "scbdtgv6/F");
	    //	    
	    t[isign][itag][iene]->Branch("ecalEneEP",        &ecalEneEP,         "ecalEneEP/F");
	   
	    TObjArray* obj = t[isign][itag][iene]->GetListOfBranches();
	    for (int ii=0; ii<obj->GetEntries(); ii++) {
	      TBranch* branch = (TBranch*)(obj->At(ii));
	      branch->SetCompressionLevel(6);
	    }
 
	  }
    }

  printf("Output files and trees created...\n");
  
  int nentries = (int)chain->GetEntries();
  if(QUICK) nentries = TMath::Min(nentries,10000);

  cout<<Form("Analyzing %d entries\n", nentries);

  TRandom3 *rgen = new TRandom3(567);
  const int prescale=20;

  int cont=1;
  for(int ientry=0; ientry<nentries; ientry++)
    {
      if( (ientry/(float)nentries) >= ((1./20)*cont)) {
	printf("%.2d%%\n",cont*100/20);
	cont++;
	for(int itag=0; itag<NTAG; itag++){
	  for(int isign=0; isign<NSIGN; isign++) {
	    for(int iene=0; iene<NENE+1; iene++) {
	      t[isign][itag][iene]->FlushBaskets();
	    }
	  }
	}
      }
      chain->GetEntry(ientry);

      if( SPAN>=0&&SPAN<=5 )   //if SPAN<6 select each span. otherwise put all spans together 
	{ if( trkSpan!=SPAN ) continue; }
      
      bool good_for_anal = true;
      if( trkInnQ>1.5 )          good_for_anal=false;
      if( (trdKHitsOnTrkTrack<8 && trdKHitsOnTrdTrack<8) ) good_for_anal=false;
      if( trdKlhrEH_ene>0.8 )   good_for_anal=false;
      if( trdKlhrEP_ene<0 )     good_for_anal=false;
      if( nTrTracks != 1 )       good_for_anal=false;
      if( rtiLifetime<0 )        good_for_anal=false;
      if( passRTI!=1 )           good_for_anal=false;
      if( runAnalysisTag<0 )     good_for_anal=false;
      if( RunTag<0xF000 )        good_for_anal=false;

      ecalCorrEneE_E = fabs(ecalCorrEneE_E);
      ecalCorrEneP_E = fabs(ecalCorrEneP_E);
      
      ecalGetCorrEne2015 = fabs(ecalGetCorrEne2015);//why this 'fabs'?

      eop = ecalEnergyD/TMath::Abs(trkDefRig);

      ecalEneEP = (ecalCorrEneE_E+ecalCorrEneP_E)/2;

      float class_ene = -9999;
      if( ENECLASS==0 ) class_ene = ecalCorrEneE_E;
      else if( ENECLASS==1 ) class_ene = ecalCorrEneP_E;
      else if( ENECLASS==2 ) class_ene = ecalEneEP;
      else if( ENECLASS==3 ) class_ene = trkDefRig;
      else if( ENECLASS==4 ) class_ene = ecalEnergyE;
      else if( ENECLASS==5 ) class_ene = ecalGetCorrEne2015;
      else if( ENECLASS==6 ) class_ene = ecalEnergyD;
      else if( ENECLASS==7 ) class_ene = 0.988*ecalCorrEneE_E;
      else if( ENECLASS==8 ) class_ene = 0.9745*ecalCorrEneE_E;
      else if( ENECLASS==9 ) class_ene = ecalGetCorrEne2017;
      //      if( class_ene<0 ) continue; //this should be removed! Otherwise in the case ENECLASS==3 no negatives are produced...
      
      int ene_index = GetEnergyIndex( class_ene );
      if( ene_index<0 ) continue;
      if( ene_index==999 ) ene_index = NENE; //Fill tree for energies greater than max energy
      
      trdclass         = -( log10(trdKlhE_ene)+2);
      trdclass_trdstd  = -( log10(trdKlhE_ene_trdstd)+2);
      
      thetaontop = accTheta[0];
      phiontop = accPhi[0];


      //BoxCoxNormalize
      double bcnlay_d[6];
      double bcnlay_ns_d[6];
      double bcnslay_d[5];
      bool binned = kFALSE;
      float layfrac[6]; for(int ilay=0; ilay<6; ilay++) layfrac[ilay] = ecalLayerFrac[ilay];
      int goodboxcox = GetBoxCoxNorm( ecalCorrEneE_E, ecalEnergyD, layfrac, bcnlay_d, bcnlay_ns_d, bcnslay_d, binned );
      if( goodboxcox!=0 ) { printf( "GetBoxCoxNorm failed\n"); return false; }
      for(int ilay=0; ilay<6; ilay++) bcnlay[ilay]     = (float)bcnlay_d[ilay];
      for(int ilay=0; ilay<6; ilay++) bcnlay_ns[ilay]  = (float)bcnlay_ns_d[ilay];
      for(int ilay=0; ilay<5; ilay++) bcnslay[ilay]    = (float)bcnslay_d[ilay];

      //TMVA simple cuts
      sclhv0 = (float)GetSCBDT(bcnlay_ns_d, "Likelihood", 0);
      sclhv5 = (float)GetSCBDT(bcnlay_d,    "Likelihood", 5);
      sclhv6 = (float)GetSCBDT(bcnlay_d,    "Likelihood", 6);
    
      sclhdv0 = (float)GetSCBDT(bcnlay_ns_d, "LikelihoodD", 0);
      sclhdv5 = (float)GetSCBDT(bcnlay_d,    "LikelihoodD", 5);
      sclhdv6 = (float)GetSCBDT(bcnlay_d,    "LikelihoodD", 6);

      scbdtv0 = (float)GetSCBDT(bcnlay_ns_d, "BDT", 0);
      scbdtv5 = (float)GetSCBDT(bcnlay_d,    "BDT", 5);
      scbdtv6 = (float)GetSCBDT(bcnlay_d,    "BDT", 6);

      scbdtgv0 = (float)GetSCBDT(bcnlay_ns_d, "BDTG", 0);
      scbdtgv5 = (float)GetSCBDT(bcnlay_d,    "BDTG", 5);
      scbdtgv6 = (float)GetSCBDT(bcnlay_d,    "BDTG", 6);

      
      //Stupid Simple cuts
      float mean=-999; 
      mean=0; for(int ilay=0; ilay<6; ilay++) mean+=bcnlay[ilay]; mean/=6; scmean      = (float)mean;
      mean=0; for(int ilay=0; ilay<5; ilay++) mean+=bcnlay[ilay]; mean/=5; scmean_5lay = (float)mean;
       
      int minlay, maxlay;
      int dummylhversion=1;
      minlay=0; maxlay=5;
      scdummylh      = (float)GetBoxCoxDummyLikelihood(ecalCorrEneE_E, bcnlay_d, minlay, maxlay, dummylhversion);
      minlay=0; maxlay=4;
      scdummylh_5lay = (float)GetBoxCoxDummyLikelihood(ecalCorrEneE_E, bcnlay_d, minlay, maxlay, dummylhversion);
       
      int tree_sign=-1;
      if( trkDefRig>0 ) tree_sign=1; else tree_sign=0;

      int tree_analtag=0;
      if( good_for_anal==false ) tree_analtag=1;    
      
      //      gDebug=1;
      t[tree_sign][tree_analtag][ene_index]->Fill();
#ifdef ALSOPRESCALED
      if( rgen->Uniform() < 1./prescale ) t[2][tree_analtag][ene_index]->Fill();
#endif
      
    }

  for(int isign=0; isign<NSIGN; isign++)
    for(int itag=0; itag<NTAG; itag++)
      for(int iene=0; iene<NENE+1; iene++)
	{
	  cout<<Form("\n%d/%d in file %s\n", (int)t[isign][itag][iene]->GetEntries(),nentries, f[isign][itag][iene]->GetName());
	  cout<<Form("Saving file %s ...",f[isign][itag][iene]->GetName());
	  f[isign][itag][iene]->cd();  t[isign][itag][iene]->Write(); f[isign][itag][iene]->Close();
	  cout<<Form("Done\n");
	}
 
  return 1;
}

bool PrintSingleEne(int __eneclass__){
  
  if( __eneclass__==0 ) printf("Energy Classifier %d: CorrEne(2,2) aka OldEne                    \n", __eneclass__);
  else if( __eneclass__==1 ) printf("Energy Classifier %d: CorrEne(0,2) aka EneP                 \n", __eneclass__);
  else if( __eneclass__==2 ) printf("Energy Classifier %d: <EneP,EneE>                           \n", __eneclass__);
  else if( __eneclass__==3 ) printf("Energy Classifier %d: Rigidity                              \n", __eneclass__);
  else if( __eneclass__==4 ) printf("Energy Classifier %d: EnergyE                               \n", __eneclass__);
  else if( __eneclass__==5 ) printf("Energy Classifier %d: CorrEne2015 aka NewEne                \n", __eneclass__);
  else if( __eneclass__==6 ) printf("Energy Classifier %d: EnergyD                               \n", __eneclass__);
  else if( __eneclass__==7 ) printf("Energy Classifier %d: 0.988*CorrEne(2,2) aka 0.988*OldEne   \n", __eneclass__);
  else if( __eneclass__==8 ) printf("Energy Classifier %d: 0.9745*CorrEne(2,2) aka 0.9745*OldEne \n", __eneclass__);
  else if( __eneclass__==9 ) printf("Energy Classifier %d: CorrEne2017 by N.Zimmermann (May2017) \n", __eneclass__);
  else {
    printf("EnergyClassifier (%d) not valid\n", __eneclass__ );
    return false;
  }

  return true;
}

void PrintVariousEne(){

  for (int ii=0; ii<10; ii++) {
    PrintSingleEne(ii);
  }

  return;
}


AnyOption *opt;

bool GetOptions(int argc, char** argv){

  opt = new AnyOption();

  //**************
  //*****Add Usage
  //**************
  opt->addUsage( Form("Usage: %s [options] [arguments] --filelist <filelist> <index of the first file to process> <index of the last file to process>\n", argv[0]) );
  opt->addUsage( "" );
  opt->addUsage( "Options: " );
  opt->addUsage( "  -h, --help ........................... Print this help " );
  opt->addUsage( "  -v, --verbose ........................ Verbose Mode" );
  opt->addUsage( "  -d, --debug .......................... Debug ON" );
  opt->addUsage( "  " );
  opt->addUsage( "  --span <N> ........................... To process just events with a particular Tracker span, N in [0,5]" );
  opt->addUsage( "  --outputdir <path> ................... To specify where to write the output files (\"./\" is the default)" );
  opt->addUsage( "  --filelist <filelist> ................ To specify the filelist to be processed (MANDATORY)" );
  opt->addUsage( "  --eneclass <N> ....................... The energy classifier (by code) used to split the output files (0 is the default)" );
  opt->addUsage( "  --ntupleversion <N> .................. The ntuple format version (2 is the default)" );
  opt->addUsage( "  " );
  opt->addUsage( "  -q, --quick .......................... Quick mode (useful for debug) ON" );
  opt->addUsage( "  --printene ........................... Print the energies code number" );
  opt->addUsage( "  " );
  opt->addUsage( "  --addfriend <path to Friend files> ... For each file to be processed, in the filelist, add a friend (with same name) from the dir passed" );
  opt->addUsage( "  " );
  
  //***********
  //set Flags
  //***********
  opt->setFlag( "help",   'h' );
  opt->setFlag( "verbose",'v' );
  opt->setFlag( "debug",  'd' );
  opt->setFlag( "quick",  'q' );
  opt->setFlag( "printene" );
  opt->setOption( "span" );
  opt->setOption( "outputdir" );
  opt->setOption( "filelist" );
  opt->setOption( "eneclass" );
  opt->setOption( "ntupleversion" );
  opt->setOption( "addfriend" );

  //****************
  //Get Line Command
  //****************
  opt->processCommandArgs( argc, argv );

  //*************
  //Get Options
  //*************
  if( opt->getFlag( "help" ) || opt->getFlag( 'h' ) ){ opt->printUsage(); exit(2); }

  if( opt->getFlag( "printene" ) ) {
    PrintVariousEne();
    return false;
  }

  VERBOSE=false;
  if( opt->getFlag( "verbose" ) || opt->getFlag( 'v' ) ) VERBOSE = true;
  printf("Verbose set: %s\n", VERBOSE ? "true" : "false" );

  DEBUG = false;
  if( opt->getFlag( "debug" ) || opt->getFlag( 'd' ) ) DEBUG = true;
  printf("Debug set: %s\n", DEBUG ? "true" : "false" );

  QUICK = false;
  if( opt->getFlag( "quick" ) || opt->getFlag( 'd' ) ) QUICK = true;
  printf("QUICK set: %s\n", QUICK ? "true" : "false" );
 
  OUTPUTDIR.assign(".");
  if( opt->getValue("outputdir") )
    { OUTPUTDIR.assign( opt->getValue("outputdir") ); }
  printf("Outputdir: %s\n", OUTPUTDIR.c_str() );

  ENECLASS=0;
  if( opt->getValue("eneclass") ) {
    ENECLASS = atoi( opt->getValue("eneclass") );
  }
  if (!PrintSingleEne(ENECLASS)) return false;

  if( opt->getValue("filelist") ) {
    FILELIST.assign( opt->getValue("filelist") );
    printf("Filelist: %s\n", FILELIST.c_str() );
  }
  else {
    printf("The \"--filelist\" option is MANDATORY\n");
    return false;
  }
     
  SPAN=-1;
  if( !opt->getValue("span") )
    { printf("No span specified. No selection on span\n"); }
  else 
    {
      SPAN = atoi( opt->getValue("span") );
      if( SPAN<0 || SPAN>5 ) {printf ("Please specify a tracker span in the range [0,5]\n"); return false; }
      else printf("Tracker Span: %d\n", SPAN );
    }
  
  NTUPLEVERSION = 2;//default (current)
  if( opt->getValue("ntupleversion") )
    { NTUPLEVERSION=atoi(opt->getValue("ntupleversion")); }
  printf("Ntuple version: %d\n", NTUPLEVERSION );

  ADDFRIEND=false;
  DIRFRIEND.assign(".");
  if(opt->getValue("addfriend")){
    ADDFRIEND=true;
    DIRFRIEND.assign(opt->getValue("addfriend"));
    printf("We'll add friend files from: %s\n", DIRFRIEND.c_str());
  }

  // printf("%d\n", opt->getArgc());
  // for (int ii=0; ii<opt->getArgc(); ii++) {
  //   printf("%d) %s\n", ii, opt->getArgv(ii));
  // }
  
  //*************                                                                                                                            
  //Get Arguments                                                                                                                            
  //************                                                                                                                             
  switch (opt->getArgc()){
    
  case 2:
    list_first_line = atoi( opt->getArgv(0) );
    list_last_line = atoi( opt->getArgv(1) );
    break;
    
  default:
    opt->printUsage();
    exit(-1);
    break;
    
  }
  
  return true;
}

int main(int argc, char** argv){
  
  if( !GetOptions(argc,argv) ) return 1;

  const int nbins=77;
  double PaperoFourBinning[nbins] = {
    0.50, 0.65, 0.82, 1.01, 1.22,
    1.46, 1.72, 2.00, 2.31, 2.65,
    3.00, 3.36, 3.73, 4.12, 4.54,
    5.00, 5.49, 6.00, 6.54, 7.10,
    7.69, 8.30, 8.95, 9.62, 10.32,
    11.04, 11.80, 12.59, 13.41, 14.25,
    15.14, 16.05, 17.00, 17.98, 18.99,
    20.04, 21.13, 22.25, 23.42, 24.62,
    25.90, 27.25, 28.68, 30.21, 31.82,
    33.53, 35.36, 37.31, 39.39, 41.61,
    44.00, 46.57, 49.33, 52.33, 55.58,
    59.13, 63.02, 67.30, 72.05, 77.37,
    83.36, 90.19, 98.08, 107.34, 118.42,
    132.11, 148.81, 169.86, 197.69, 237.16,
    290.0, 370.0, 500.0, 700.0, 1000.0,
    1500.00, 2200.00
  };

  const int nbins_antip = 31;
  double AntiProtonsBinning[nbins_antip] = {
    1, 1.33, 1.71, 2.15, 2.67, 3.29, 4.02, 4.88, 5.9, 7.09,
    8.48, 10.1, 12, 14.1, 16.6, 19.5, 22.8, 26.7, 31.1, 36.1,
    41.9, 48.5, 56.1, 64.8, 74.9, 86.5, 100, 125, 175, 259,
    450
  };
  
  NENE = (nbins-1);
  if (ENECLASS==3) NENE = (nbins_antip-1);
  MinEne = new Double_t[NENE];
  MaxEne = new Double_t[NENE];
  for (int ii=0; ii<NENE; ii++) {
    if (ENECLASS==3) {
      MinEne[ii]=AntiProtonsBinning[ii];
      MaxEne[ii]=AntiProtonsBinning[ii+1];
    }
    else {
      MinEne[ii]=PaperoFourBinning[ii];
      MaxEne[ii]=PaperoFourBinning[ii+1];
    }
  }
  
  for(int iene=0; iene<NENE; iene++)
    cout<<Form(" iene:%.2d) [%.2f - %.2f]", iene, MinEne[iene], MaxEne[iene] )<<endl;
  cout<<endl;
  
  Build_superskimmed();
  
  return 0;
}