////////////////////////////////////////////////////////////////////////////////////
//
//             _____False Source Method macro_____
//
//  Take as input root files with hillas parameters, perform derotation (optional)
//  and generates the alpha histograms for On and Off for different positions of the
//  source. These histos can be ploted with signal.C and signalPoint.C
//
//                 Ester Aliu   <aliu@ifae.es>
//                 Oscar Blanch <blanch@ifae.es>
//                 Javier Rico  <jrico@ifae.es>
////////////////////////////////////////////////////////////////////////////////////

#include <fstream>
#include <iostream>

#include "TString.h"
#include "TH2F.h"
#include "TFile.h"

#include "MHillasSrcCalc.h"
#include "MSrcPosCam.h"
#include "MSrcTranslate.h"
#include "MHillasSrc.h"
#include "MSrcRotate.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MHillas.h"
#include "MReadTree.h"
#include "MEvtLoop.h"
#include "MLog.h"
#include "MArgs.h"

using namespace std;

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

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

TString onFile;
TString offFile;
TString outputFile;
Bool_t  kRotate=1;
Bool_t  kSrcPolicy=kFALSE;
Double_t fRA= -1.;
Double_t fDEC= -1.;
ULong_t nmaxevents=999999999;

// cuts (default values to be overriden by datacards)
Float_t minsize   = 0;        // minimum size (# of photons)
Float_t maxsize   = 9999999;  // maximum size (# of photons)
Float_t mindist   = 0;        // minimum distance cut (deg)
Float_t maxdist   = 10;       // maximum distance cut (deg)
Float_t minwidth  = 0.;       // minimum width cut (deg)
Float_t maxwidth  = 10;       // maximum width cut (deg)
Float_t minlength = 0;        // minimum length cut (deg)
Float_t maxlength = 10;       // maximum length cut (deg)
Float_t maxXY     = maxdist;  // maximum X and Y distance from the camera center

// binning
Int_t binning = 21;           // number of bins in false source search (one coordinate) 
Int_t nbin    = 18;           // number of bins for alpha plots
Float_t field = 2.;           // width of the examined field (degrees)

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

const TString defaultcard="falsesource.datacard";

const Float_t alphamin = 0.;       // minimum value in alpha (degrees) 
const Float_t alphamax = 90.;      // maximum value in alpha (degrees)
const Float_t conver   = 189./0.6; // conversion factor degrees to mm
const Float_t srclimit = field*conver;
const Float_t srcstep  = 3.;
const Int_t   nsrcbin  = (Int_t) (srclimit/srcstep);

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

static void Usage()
{
  gLog <<endl;
  gLog << "Usage is:" << endl;
  gLog << "   falseSource [-h] [-?] <datacards>" << endl << endl;
  gLog << "     <datacards>: datacards file name (dafault falsesource.datacards)" << 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;
    }
  falseSource();
}

void falseSource()
{
  // CUTS

  // variable declaration  
  Float_t xpos;
  Float_t ypos;  
  Char_t name[20]; 
  Char_t title[50];   
  Float_t stepDegrees = field/(binning-1);    // degrees
  Float_t step        = stepDegrees*conver;   // mm
  Int_t   ival        =  (Int_t) ceil((Float_t)(binning-1)/2); // index bound
  
  // Hillas variables
  Float_t alpha;
  Float_t size;
  Float_t dist;  
  Float_t length;
  Float_t width;
  Float_t meanX;
  Float_t meanY;
  
  // create the histograms (empty)
  TH1F* hOnAlpha[binning][binning];
  TH1F* hOffAlpha[binning][binning];
  TH2F* hOnSrcPos[binning];
  TH2F* hOffSrcPos[binning];

  for(Int_t i = -ival; i <= ival ; i++){  
    for(Int_t j = -ival; j <= ival ; j++){
      
      sprintf(name,"hOnAlpha[%d][%d]", i, j);
      sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j);
      hOnAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax);
      
      sprintf(name,"hOffAlpha[%d][%d]", i, j);
      sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j);
      hOffAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax);

      if(j==0)
	{
	  sprintf(name,"hOnSrcPos[%d]", i);
	  sprintf(title,"Source position (On data) (%d)", i);
	  hOnSrcPos[i+ival] = new TH2F(name, title, nsrcbin, -srclimit, srclimit, nsrcbin, -srclimit, srclimit);
	  
	  sprintf(name,"hOffSrcPos[%d]", i);
	  sprintf(title,"Source position (Off data) (%d)", i);
	  hOffSrcPos[i+ival] = new TH2F(name, title, nsrcbin, -srclimit, srclimit, nsrcbin, -srclimit, srclimit);
	}
    }
  }
  

  
  /******************************************************************
                           FILL ON-DATA HISTOS
  *******************************************************************/
  
  // source dependent hillas containers/tasks
  MSrcPosCam*     source[binning][binning];
  MHillasSrc*     hillasSrc[binning][binning];
  MSrcTranslate*  srctranslate[binning][binning];
  MSrcRotate*     srcrotate[binning][binning];
  MHillasSrcCalc* csrc1[binning][binning];

  // normal containers/classes
  MParList  plist;  
  MTaskList tlist;
  MHillas   hillas;

  plist.AddToList(&tlist);
  plist.AddToList(&hillas);
    
  MReadTree read("Parameters", onFile);
  read.DisableAutoScheme();  
  tlist.AddToList(&read);
    
  char sourceName[100];  
  char hillasSrcName[100];
  char translateName[100];  
  char rotateName[100];

  // create the tasks/containers for the different source positions
  for(Int_t i=-ival;i<=ival;i++)
    {    
      xpos = step*i;      
      for(Int_t j=-ival;j<=ival;j++)
	{
	  ypos = step*j;
	  
	  // name the different containers
	  if (i<0 && j<0)
	    {
	      Int_t k = -i;
	      Int_t l = -j;
	      sprintf(sourceName,    "MSrcPosCam_Min%dMin%d", k, l);
	      sprintf(hillasSrcName, "MHillasSrc_Min%dMin%d", k, l);
	      sprintf(translateName, "MSrcTranslate_Min%dMin%d", k, l);
	      sprintf(rotateName,    "MSrcRotate_Min%dMin%d", k, l);
	    }
	  
	  if (i<0 && j>=0)
	    {
	      Int_t k = -i;
	      sprintf(sourceName,    "MSrcPosCam_Min%d%d", k, j);
	      sprintf(hillasSrcName, "MHillasSrc_Min%d%d", k, j);
	      sprintf(translateName, "MSrcTranslate_Min%d%d", k, j);
	      sprintf(rotateName,    "MSrcRotate_Min%d%d", k, j);
	    }
	  
	  if (i>=0 && j<0)
	    {
	      Int_t l = -j;
	      sprintf(sourceName,    "MSrcPosCam_%dMin%d", i, l);
	      sprintf(hillasSrcName, "MHillasSrc_%dMin%d", i, l);
	      sprintf(translateName, "MSrcTranslate_%dMin%d", i, l);
	      sprintf(rotateName,    "MSrcRotate_%dMin%d", i, l);
	    }
	  
	  if (i>=0 && j>= 0)
	    {
	      sprintf(sourceName,    "MSrcPosCam_%d%d", i, j);
	      sprintf(hillasSrcName, "MHillasSrc_%d%d", i, j);
	      sprintf(translateName, "MSrcTranslate_%d%d", i, j);
	      sprintf(rotateName,    "MSrcRotate_%d%d", i, j);
	    }

	  // include containers in parameter list
	  source[i+ival][j+ival] = new MSrcPosCam(sourceName);
	  plist.AddToList(source[i+ival][j+ival]);
	  	  
	  hillasSrc[i+ival][j+ival] = new MHillasSrc(hillasSrcName);	  
	  plist.AddToList(hillasSrc[i+ival][j+ival]);      
	  
	  // define the different tasks and include them in the task list
	  srctranslate[i+ival][j+ival] = new MSrcTranslate("MSrcPosCam",sourceName,"MDCA",translateName);	  
	  srctranslate[i+ival][j+ival]->SetTranslation(xpos,ypos);	  
	  srctranslate[i+ival][j+ival]->SetRelativeTranslation(kSrcPolicy);
	  srctranslate[i+ival][j+ival]->SetInternalHistoBinSize(5.);
	  srctranslate[i+ival][j+ival]->SetMode(MSrcPlace::kOn);
      
	  tlist.AddToList(srctranslate[i+ival][j+ival]);  

	  if(kRotate && !kSrcPolicy)
	    {
	      srcrotate[i+ival][j+ival] = new MSrcRotate(sourceName,sourceName,"MDCA",rotateName);
	      srcrotate[i+ival][j+ival]->SetRAandDEC(fRA,fDEC);	  

	      srctranslate[i+ival][j+ival]->SetCreateHisto(kFALSE);
	      srcrotate[i+ival][j+ival]->SetMode(MSrcPlace::kOn);	      
	      srcrotate[i+ival][j+ival]->SetInternalHistoBinSize(5.);

	      tlist.AddToList(srcrotate[i+ival][j+ival]);  
	    }
	  else
	    srcrotate[i+ival][j+ival] = NULL;


	  TString HilName = "MHillas";	  
	  csrc1[i+ival][j+ival] = new MHillasSrcCalc(sourceName, hillasSrcName);
	  csrc1[i+ival][j+ival]->SetInput(HilName);  
	  
	  tlist.AddToList(csrc1[i+ival][j+ival]);  
	  
	} // loop on j (y coordinate)
    } // loop on i (x coordinate)
  
  
  // Eventloop
  MEvtLoop evtloop;  
  evtloop.SetParList(&plist);
  //  MProgressBar bar;
  //  evtloop.SetProgressBar(&bar);  
  if (!evtloop.PreProcess())
    return;
  
  // select the events passing the cuts
  UInt_t nread=0;
  while (tlist.Process())
    {
      for(Int_t i = -ival; i <= ival ; i++)	
	for(Int_t j = -ival; j <= ival ; j++)
	  {
	    width  = hillas.GetWidth()/conver;
	    length = hillas.GetLength()/conver;
	    size   = hillas.GetSize();
	    dist   = hillasSrc[i+ival][j+ival]->GetDist()/conver;
	    meanX  = hillas.GetMeanX()/conver;
	    meanY  = hillas.GetMeanY()/conver;
	    
	    if (width<maxwidth && length<maxlength && dist>mindist && dist<maxdist && size>minsize && meanX<maxXY && meanY<maxXY)
	      {
		alpha = hillasSrc[i+ival][j+ival]->GetAlpha();
		hOnAlpha[i+ival][j+ival]->Fill(TMath::Abs(alpha));
	      }	

	    // fill source position histo
	    if(j==0)
	      hOnSrcPos[i+ival]->Fill(source[i+ival][j+ival]->GetX(),source[i+ival][j+ival]->GetY());
	  }
      if(++nread>nmaxevents) break;
    }
      
  // end 
  evtloop.PostProcess();    
  tlist.PrintStatistics();


  /******************************************************************
                        FILL OFF-DATA HISTOS
  *******************************************************************/
 for(Int_t i=-ival;i<=ival;i++)
   for(Int_t j=-ival;j<=ival;j++)
     {
       if(srcrotate[i+ival][j+ival])
	 srcrotate[i+ival][j+ival]->SetMode(MSrcPlace::kOff);
       srctranslate[i+ival][j+ival]->SetMode(MSrcPlace::kOff);
     }

  
  MReadTree read2("Parameters", offFile);
  read2.DisableAutoScheme();
  
  tlist.AddToListBefore(&read2, &read, "All");
  tlist.RemoveFromList(&read);
  
  if (!evtloop.PreProcess())
    return;

  nread=0;
  while (tlist.Process())
    {
      for(Int_t i = -ival; i <= ival ; i++)
	for(Int_t j = -ival; j <= ival ; j++)
	  {	    
	    width  = hillas.GetWidth()/conver;
	    length = hillas.GetLength()/conver;
	    size   = hillas.GetSize();   
	    dist   = hillasSrc[i+ival][j+ival]->GetDist()/conver;
	    meanX  = hillas.GetMeanX()/conver;
	    meanY  = hillas.GetMeanY()/conver;
	    
	    if (width<maxwidth && length<maxlength && dist>mindist && dist<maxdist && size>minsize && meanX<maxXY && meanY<maxXY)
	      {
		alpha = hillasSrc[i+ival][j+ival]->GetAlpha();
		hOffAlpha[i+ival][j+ival]->Fill(TMath::Abs(alpha));	      
	      } 	

	    // fill source position histo
	    if(j==0)
	      hOffSrcPos[i+ival]->Fill(source[i+ival][j+ival]->GetX(),source[i+ival][j+ival]->GetY());
	  }
      if(++nread>nmaxevents) break;
    }
        
  evtloop.PostProcess();    

  // Save results
  TFile *f = new TFile(outputFile, "RECREATE");
  
  for(Int_t i = -ival; i <= ival ; i++){  
    for(Int_t j = -ival; j <= ival ; j++){
      hOnAlpha[i+ival][j+ival]->Write();
      hOffAlpha[i+ival][j+ival]->Write();
      if(j==0)
	{
	  hOnSrcPos[i+ival]->Write();
	  hOffSrcPos[i+ival]->Write();
	}
    }
  }

  f->Close();

 for(Int_t i=-ival;i<=ival;i++)
   for(Int_t j=-ival;j<=ival;j++)
     {
       if(srcrotate[i+ival][j+ival])
	 delete srcrotate[i+ival][j+ival];
       delete srctranslate[i+ival][j+ival];
       delete source[i+ival][j+ival];
       delete hillasSrc[i+ival][j+ival];
       delete csrc1[i+ival][j+ival];
     }
}
//-----------------------------------------------------------------------
  
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;

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

      // off-data file name
      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;
	}

      // rotation flag
      if(strcmp(word.Data(),"ROTFLAG")==0)
	ifun >> kRotate;

      // source movement policy flag
      if(strcmp(word.Data(),"SRCABS")==0)
	ifun >> kSrcPolicy;	  


      // source celestial coordinates 
      if(strcmp(word.Data(),"SRCCOORDS")==0)
	{
	  ifun >> fRA;
	  ifun >> fDEC;	  
	}

      // field width
      if(strcmp(word.Data(),"FIELD")==0)
	ifun >> field;

      // binning
      if(strcmp(word.Data(),"BINNING")==0)
	ifun >> binning;

      // Number of bins in alpha plots
      if(strcmp(word.Data(),"NBIN")==0)
	ifun >> nbin;

      // size cut
      if(strcmp(word.Data(),"SIZECUT")==0)
	{
	  ifun >> minsize;
	  ifun >> maxsize;	  
	}

      // dist cut
      if(strcmp(word.Data(),"DISTCUT")==0)
	{
	  ifun >> mindist;
	  ifun >> maxdist;	  
	}

      // width cut
      if(strcmp(word.Data(),"WIDTHCUT")==0)
	{
	  ifun >> minwidth;
	  ifun >> maxwidth;	  
	}
      
      // length cut
      if(strcmp(word.Data(),"LENGTHCUT")==0)
	{
	  ifun >> minlength;
	  ifun >> maxlength;	  
	}

      // maxX and maxY upper cut
      if(strcmp(word.Data(),"CENTERCUT")==0)
	ifun >> maxXY;	  
    }

  // 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 file name(s): " << onFile << endl;
  cout << "Off-data file name(s): " << offFile << endl;
  cout << "Output file name: " << outputFile << endl;
  if(kSrcPolicy)
    cout << "Source position displaced with respect to the original existing one (no rotation applied)" << endl;
  else
    {
      cout << "De-rotation flag " << kRotate << endl;
      cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl;
    }
  cout << "Field width (degrees): " << field << endl;
  cout << "Field binning: " << binning << endl;
  cout << "Number of alpha plot bins: " << nbin << endl;
  cout << "Size cuts (# of photons): ("<<minsize<<","<<maxsize<<")" << endl;
  cout << "Dist cuts (degrees): ("<<mindist<<","<<maxdist<<")" << endl;
  cout << "Length cuts (degrees): ("<<minlength<<","<<maxlength<<")" << endl;
  cout << "Width cuts (degrees): ("<<minwidth<<","<<maxwidth<<")" << endl;
  cout << "maxX and maxY upper cut (degrees): " << maxXY << endl;
  cout << "***********" << endl << endl;

  return kTRUE;
}