////////////////////////////////////////////////////////////////////////////////////
//
//             _____ Optimize cuts _____
//
//  Take as input root files with hillas parameters and compute the optimal cuts by
//  mapping the width/length plane (only upper cuts so far)
//
//                 Jose Flix    <jflix@ifae.es>
//                 Javier Rico  <jrico@ifae.es>
////////////////////////////////////////////////////////////////////////////////////

#include <fstream>
#include <iostream>

#include "TString.h"
#include "TChain.h"
#include "TH1F.h"

#include "MArgs.h"
#include "MLog.h"

using namespace std;

Bool_t readDatacards(TString& filename);
void optimizeCuts();


//-----------------------------------------------------------------------------
// declaration of variables read from datacards
//-----------------------------------------------------------------------------

TString  onFile;
TString  offFile;
TString  outputFile;
ULong_t  nmaxevents=999999999;
Float_t  lowdistcut=0;
Float_t  uppdistcut=10;
Float_t  lowsizecut=1800;
Float_t  uppsizecut=1800;
Int_t    nsizebins=1;
UInt_t   onRate=1;
Float_t  lowlgcut=0.01;
Float_t  upplgcut=0.3;
Float_t  lgstep=0.005;
Float_t  lowwdcut=0.01;
Float_t  uppwdcut=0.3;
Float_t  wdstep=0.005;

//-----------------------------------------------------------------------------
// constants
//-----------------------------------------------------------------------------

const TString defaultcard="optimizecuts.datacard";
const Float_t conver   = 189./0.6; // conversion factor degrees to mm
const Int_t nbins=18; //number of bins in alpha plots
const Float_t frontiere=20.; // [deg] value below (above) wich signal (background is computed (normalized)

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

static void Usage()
{
  gLog <<endl;
  gLog << "Usage is:" << endl;
  gLog << "   optimizeCuts [-h] [-?] <datacards>" << endl << endl;
  gLog << "     <datacards>: datacards file name (default " << defaultcard <<")" << endl;
  gLog << "     -?/-h: This help" << endl << endl;
}

//-----------------------------------------------------------------------------
int main(int argc, char **argv)
{
  // evaluate arguments
  MArgs arg(argc, argv);
  if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
    {
      Usage();
      return -1;
    }

  TString datacard  = arg.GetArgumentStr(0);
  if(!datacard.Length())
    datacard = defaultcard;

  if(!readDatacards(datacard))
    {
      cout << "Error reading datacards. Stoping" << endl;
      return -1;
    }
  optimizeCuts();
}

//-----------------------------------------------------------------------------
void optimizeCuts()
{
  // create the TChains for ON and OFF data
  TChain* ton= new TChain("Parameters");
  TChain* toff= new TChain("Parameters");
  if(!ton->Add(onFile))
    {
      cout << "On-data file not found or not valid format" << endl;
      return;
    }
  if(!toff->Add(offFile))
    {
      cout << "Off-data file not found or not valid format" << endl;
      return;
    }

  ofstream fout;
  fout.open(outputFile.Data());

  // define aliases
  ton->SetAlias("length","MHillas.fLength*0.6/189");
  ton->SetAlias("width","MHillas.fWidth*0.6/189");
  //ton->SetAlias("length","MHillas.fLength");
  //ton->SetAlias("width","MHillas.fWidth");
  ton->SetAlias("dist","MHillasSrc.fDist*0.6/189");
  ton->SetAlias("conc","MNewImagePar.fConc");
  ton->SetAlias("size","MHillas.fSize");
  ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
  ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
  ton->SetAlias("leak","MNewImagePar.fInnerLeakage1");
  
  toff->SetAlias("length","MHillas.fLength*0.6/189");
  toff->SetAlias("width","MHillas.fWidth*0.6/189");
  //toff->SetAlias("length","MHillas.fLength");
  //toff->SetAlias("width","MHillas.fWidth");
  toff->SetAlias("dist","MHillasSrc.fDist*0.6/189");
  toff->SetAlias("conc","MNewImagePar.fConc");
  toff->SetAlias("size","MHillas.fSize");
  toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
  toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
  toff->SetAlias("leak","MNewImagePar.fInnerLeakage1");


  // general cut
  char gencut[256];
  char evecut[256];
  sprintf(evecut,"event%%%d==0",onRate);
  sprintf(gencut,"dist>%f && dist<%f",lowdistcut,uppdistcut);
  cout << "General cut " << gencut << " (" << evecut << " for on-data)" << endl;

  Bool_t isDiff = kFALSE;
  if(nsizebins<0)
    {
      nsizebins*=-1;
      isDiff=kTRUE;
    }
  else
    nsizebins+=1;

  // loop on size, width and length cuts
  for(Int_t isb=0;isb<nsizebins;isb++) // loop on size
    {      
      cout << isb << endl;      
      Float_t minsize = isDiff ? 
	TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) :
	TMath::Power(10,TMath::Log10(lowsizecut)+isb*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/(nsizebins-1));
      Float_t maxsize = isDiff ? TMath::Power(10,TMath::Log10(lowsizecut)+(isb+1)*(TMath::Log10(uppsizecut)-TMath::Log10(lowsizecut))/nsizebins) : 9999999999.;
      
      Float_t length=lowlgcut;
      while(length<=upplgcut) // loop on length
	{
	  Float_t width=lowwdcut;
	  while(width<=uppwdcut) // loop on width
	    {
	      TH1F* onhisto  = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
	      TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90);    
	      
	      // define the width/length cut
	      char varcut[256];
	      sprintf(varcut,"size>%f && size<%f && length<%f && width<%f",minsize,maxsize,length,width);
	      cout << "Cutting " << varcut << endl;
	      
	      // define the on/off data cut
	      char offcut[1024];
	      sprintf(offcut,"%s && %s",gencut,varcut);
	      char oncut[1024];
	      sprintf(oncut,"%s && %s",evecut,offcut);
	      
	      // Fill the histos
	      ton->Draw("alpha>>OnHist",oncut);
	      toff->Draw("alpha>>OffHist",offcut);
	      
	      // Normalize histos
	      const Int_t inibin = (Int_t)(frontiere/90.*nbins+1);
	      Float_t level=0;
	      Float_t leveloff=0;
	      for(Int_t ibin = inibin; ibin<=nbins;ibin++)
		{
		  level+=onhisto->GetBinContent(ibin);
		  leveloff+=offhisto->GetBinContent(ibin);
		}
	      offhisto->Sumw2(); // needed to compute correct errors after normalization
	      const Float_t norm = leveloff ? level/leveloff: 1;
	      cout << "Normalizing by factor " << norm <<endl;
	      offhisto->Scale(norm);    
	      
	      // compute significance/excess
	      Float_t sig=0,bg=0,esig=0,ebg=0;
	      Float_t significance=0;
	      const Int_t signbins = inibin-1;
	      for(Int_t ibin = 1; ibin<=signbins;ibin++)
		{
		  sig  += onhisto->GetBinContent(ibin);
		  esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
		  bg   += offhisto->GetBinContent(ibin);
		  ebg  += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
		}
	      Float_t error= TMath::Sqrt(esig+ebg);
	      significance = error ? (sig-bg)/error : 0; 
	      
	      cout << "Excess: " << sig-bg << endl;
	      cout << "N bkg: "  << bg << endl; 
	      cout << "Significance: " << significance << endl;	  
	      
	      // save the result in file
	      fout << minsize << '\t' 
		   << maxsize << '\t' 
		   << width << '\t' 
		   << length << '\t'
		   << sig-bg << '\t'
		   << significance << '\t' << endl;
	      
	      delete onhisto;
	      delete offhisto;
	      
	      width+=wdstep;
	    }
	  length+=lgstep;
	}
    }
  fout.close();
}
//-----------------------------------------------------------------------
  
Bool_t readDatacards(TString& filename)
{
  ifstream ifun(filename.Data());
  if(!ifun)
    {
      cout << "File " << filename << " not found" << endl;
      return kFALSE;
    }

  TString word;
  
  while(ifun >> word)
    {
      // skip comments
      if(word[0]=='/' && word[1]=='/')
	{
	  while(ifun.get()!='\n'); // skip line
	  continue;
	}
      
      // number of events
      if(strcmp(word.Data(),"NEVENTS")==0)
	ifun >> nmaxevents;
      
      // number of events
      if(strcmp(word.Data(),"ONRATE")==0)
	ifun >> onRate;
      
      // input file name (on data)
      if(strcmp(word.Data(),"ONFILES")==0)
	{
	  if(onFile.Length())
	    cout << "readDataCards Warning: overriding on-data file name" << endl;
	  ifun >> onFile;
	}

      // input file name (off data)
      if(strcmp(word.Data(),"OFFFILES")==0)
	{
	  if(offFile.Length())
	    cout << "readDataCards Warning: overriding off-data file name" << endl;
	  ifun >> offFile;
	}

      // output file name
      if(strcmp(word.Data(),"OUTFILE")==0)
	{
	  if(outputFile.Length())
	    cout << "readDataCards Warning: overriding output file name" << endl;
	  ifun >> outputFile;
	}

      // size cuts
      if(strcmp(word.Data(),"SIZECUTS")==0)
	{
	  ifun >> lowsizecut;
	  ifun >> uppsizecut;
	  ifun >> nsizebins;
	}

      // dist cuts
      if(strcmp(word.Data(),"DISTCUTS")==0)
	{
	  ifun >> lowdistcut ; 
	  ifun >> uppdistcut ; 
	}

      // length cuts
      if(strcmp(word.Data(),"LENGTHCUTS")==0)
	{
	  ifun >> lowlgcut;
	  ifun >> upplgcut;
	  ifun >> lgstep;
	}

      // width cuts
      if(strcmp(word.Data(),"WIDTHCUTS")==0)
	{
	  ifun >> lowwdcut;
	  ifun >> uppwdcut;
	  ifun >> wdstep;
	}

    }
  
  // check compulsory values
  if(!onFile.Length())
    {
      cout << "No on-data file name specified" << endl;
      return kFALSE;
    }
  
  if(!offFile.Length())
    {
      cout << "No off-data file name specified" << endl;
      return kFALSE;
    }
    
  if(!outputFile.Length())
    {
      cout << "No output file name specified" << endl;
      return kFALSE;
    }

  // Dump read values
  cout << "************************************************" << endl;
  cout << "* Datacards read from file " << filename << endl;
  cout << "************************************************" << endl;
  cout << "Maximum number of input events: " << nmaxevents << endl;
  cout << "On-data acceptance rate: 1/" << onRate << endl;
  cout << "On-data input file name(s): " << onFile << endl;
  cout << "Off-data input file name(s): " << offFile << endl;
  cout << "Output file name: " << outputFile << endl;
  cout << "Dist cuts (deg): [" << lowdistcut<< ","<<uppdistcut<< "]" << endl;
  cout << "Scanning size range: [" << lowsizecut << "," << uppsizecut << "]"<< " in " << nsizebins << " steps" << endl;
  cout << "Scanning length range (deg): [" << lowlgcut<< ","<<upplgcut<< "], with step (deg): "<< lgstep << endl;
  cout << "Scanning width range (deg): [" << lowwdcut<< ","<<uppwdcut<< "], with step (deg): "<< wdstep << endl;
  cout << "***********" << endl << endl;

  return kTRUE;
}