////////////////////////////////////////////////////////////////////////////////////
//
//             _____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 "TH1F.h"
#include "TFile.h"

#include "MHillasSrcCalc.h"
#include "MSrcPosCam.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;
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

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

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];
  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);

    }
  }
  

  
  /******************************************************************
                           FILL ON-DATA HISTOS
  *******************************************************************/
  
  // source dependent hillas containers/tasks
  MHillasSrcCalc* csrc1[binning][binning];
  MSrcPosCam*     source[binning][binning];
  MHillasSrc*     hillasSrc[binning][binning];
  MSrcRotate*     srcrotate[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);
    
  Int_t k,l;  
  char sourceName[100];  
  char hillasSrcName[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)
	    {
	      k = -i;
	      l = -j;
	      sprintf(sourceName,    "MSrcPosCam_Min%dMin%d", k, l);
	      sprintf(hillasSrcName, "MHillasSrc_Min%dMin%d", k, l);
	    }
	  
	  if (i<0 && j>=0)
	    {
	      k = -i;
	      sprintf(sourceName,    "MSrcPosCam_Min%d%d", k, j);
	      sprintf(hillasSrcName, "MHillasSrc_Min%d%d", k, j);
	    }
	  
	  if (i>=0 && j<0)
	    {
	      l = -j;
	      sprintf(sourceName,    "MSrcPosCam_%dMin%d", i, l);
	      sprintf(hillasSrcName, "MHillasSrc_%dMin%d", i, l);
	    }
	  
	  if (i>=0 && j>= 0)
	    {
	      sprintf(sourceName,    "MSrcPosCam_%d%d", i, j);
	      sprintf(hillasSrcName, "MHillasSrc_%d%d", i, j);
	    }

	  // include containers in parameter list
	  source[i+ival][j+ival] = new MSrcPosCam(sourceName);
	  source[i+ival][j+ival]->SetXY(xpos, ypos);	  
	  plist.AddToList(source[i+ival][j+ival]);
	  source[i+ival][j+ival]->Print();
	  	  
	  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
	  if(kRotate)
	    {
	      srcrotate[i+ival][j+ival] = new MSrcRotate(sourceName);
	      srcrotate[i+ival][j+ival]->SetRAandDEC(fRA,fDEC);	  
	      tlist.AddToList(srcrotate[i+ival][j+ival]);  
	    }

	  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(abs(alpha));
	      }	
	  }
      if(++nread>nmaxevents) break;
    }
      
  // end 
  evtloop.PostProcess();    
  tlist.PrintStatistics();


  /******************************************************************
                        FILL OFF-DATA HISTOS
  *******************************************************************/
  
  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(abs(alpha));	      
	      } 	
	  }
      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();
    }
  }

  f->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;

      // 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 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;
  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;
}
