////////////////////////////////////////////////////////////////////////////////////
//
//             _____ 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  sizecut=1800;
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("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)");
  
  toff->SetAlias("length","MHillas.fLength*0.6/189");
  toff->SetAlias("width","MHillas.fWidth*0.6/189");
  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)");


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

  // loop on widht and length cuts
  Float_t length=lowlgcut;
  while(length<=upplgcut)
    {
      Float_t width=lowwdcut;
      while(width<=uppwdcut)
	{
	  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,"length<%f && width<%f",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 << 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 cut
      if(strcmp(word.Data(),"SIZECUT")==0)
	ifun >> sizecut;

      // 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 << "Size cut: " << sizecut << endl;
  cout << "Dist cuts (deg): [" << lowdistcut<< ","<<uppdistcut<< "]" << 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;
}
