Double_t ChiSquareNDof(TH1D *h1, TH1D *h2)
{
    Double_t chiq = 0.;
    Double_t chi;
    Double_t error;
    Int_t nbinsnozero = 0;

    Int_t nbins = h1->GetNbinsX();
    if (nbins != h2->GetNbinsX() || nbins == 0)
	return -1;

    for (UInt_t bin=1; bin<=nbins; bin++)
    {
	error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) +
			   h2->GetBinError(bin)*h2->GetBinError(bin));
	if (error != 0)
	{
	    chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error;
	    chiq += chi*chi;
	    nbinsnozero++;
	}
    }

    return (nbinsnozero>0?chiq/nbinsnozero:0);
}

Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[])
{

  Int_t result = -1;

  Int_t lowerbin = 0;
  Int_t upperbin = numberSizeBins;
  Int_t bin;

  Int_t count = 0;

  if (size >= sizeBins[0])
    {
      while (upperbin - lowerbin > 1 && count++<=numberSizeBins)
	{
	  bin = (upperbin+lowerbin)/2;
	  if (size >= sizeBins[bin])
	    lowerbin = bin;
	  else
	    upperbin = bin;
	}
      result = count<=numberSizeBins?lowerbin:-1;
    }

  return result;

}


void lightcurve(TString f_on_name = "../HillasFiles/Mrk421/*_H.root", 
                TString f_off_name = "../HillasFiles/OffMrk421/*_H.root")
{

  gROOT->Reset();
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetTimeOffset(-3600);

  // Configuration
  const Bool_t debug = kFALSE;
  const Double_t timebin = 480.; //[sec]
  TString psname = "./20040420_Mrk421.480s.ps";

  //Constanst
  const Double_t kConvDegToRad = TMath::Pi()/180.;
  const Double_t kSec = 1e6;  //[sec/microsec]
  const Double_t kDay = 24.*60.*60.;  //[Day/sec]

  UInt_t numberTimeBins = 1;
  TArrayI numberOnEvents(1);
  TArrayD numberBkgEventsToNormOn(1);
  TArrayD timeOn(1);

  TArrayD meanTimeBinOn(1);
  TArrayD widthTimeBinOn(1);

  TArrayD zenithMinimumOn(1);
  TArrayD zenithMaximumOn(1);
  
  TObjArray coszenithHistoOn;
  TObjArray alphaHistoOn;
  TObjArray srcposHistoOn;
  

  
  const Int_t numberSizeBins = 4;
  Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]

  //cuts

  Double_t widthmin[numberSizeBins]  = { 0.06, 0.06, 0.06, 0.06 }; 
  Double_t widthmax[numberSizeBins]  = { 0.10, 0.10, 0.12, 0.12 };
  Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 }; 
  Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 };
  Double_t distmin[numberSizeBins]   = { 0.30, 0.30, 0.30, 0.30,}; 
  Double_t distmax[numberSizeBins]   = { 1.20, 1.20, 1.20, 1.20,};

  //alpha plot integration for gammas
  Double_t sigexccmin = 0.;
  Double_t sigexccmax = 20.;
  Double_t bkgnormmin = 40.;
  Double_t bkgnormmax = 80.;
  
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit();
  
  //
  // Make a loop only for the ON data:
  //
  
  MParList plist_on;
  MTaskList tlist_on;
  plist_on.AddToList(&tlist_on);
  
  // ON containers
  MGeomCamMagic geomcam;
  MRawRunHeader runheader_on;
  MRawEvtHeader evtheader_on;
  MTime              time_on;
  MHillas          hillas_on;
  MHillasSrc    hillassrc_on;
  MSrcPosCam       srcpos_on;
  MReportDrive      drive_on;

  plist_on.AddToList(&geomcam);
  plist_on.AddToList(&runheader_on);
  plist_on.AddToList(&evtheader_on);
  plist_on.AddToList(&time_on);
  plist_on.AddToList(&hillas_on);
  plist_on.AddToList(&hillassrc_on);
  plist_on.AddToList(&srcpos_on);
  plist_on.AddToList(&drive_on);

  Int_t nbins_Size = 60;
  Double_t min_Size = log10(sizeBins[0])*0.8;
  Double_t max_Size = log10(1000000)*1.2;
  TH1F *hSize_on = new TH1F ("hSizeOn","",nbins_Size,min_Size,max_Size);

  Int_t nbins_Width = 20;
  Double_t min_Width = 0.;
  Double_t max_Width = widthmax[numberSizeBins-1]*1.2;
  TH1F *hWidth_on = new TH1F ("hWidthOn","",nbins_Width,min_Width,max_Width);

  Int_t nbins_Length = 20;
  Double_t min_Length = 0.;
  Double_t max_Length =  lengthmax[numberSizeBins-1]*1.2;
  TH1F *hLength_on = new TH1F ("hLengthOn","",nbins_Length,min_Length,max_Length);

  Int_t nbins_Dist = 20;
  Double_t min_Dist = 0.;
  Double_t max_Dist = distmax[numberSizeBins-1]*1.2;
  TH1F *hDist_on = new TH1F ("hDistOn","",nbins_Dist,min_Dist,max_Dist);

  Int_t nbins_abs = 18;
  Double_t minalpha_abs = 0.;
  Double_t maxalpha_abs =90.;
  TH1F *hAlpha_on_abs = new TH1F ("hAbsAlphaOn","",nbins_abs,minalpha_abs,maxalpha_abs);

  Int_t nbins = nbins_abs*2;
  Double_t minalpha = -90.;
  Double_t maxalpha =  90.;
  TH1F *hAlpha_on = new TH1F ("hAlphaOn","",nbins,minalpha,maxalpha);

  Int_t nbins_srcpos = 400;
  Double_t minsrcpos = -600.;
  Double_t maxsrcpos =  600.;  //[mm]  //!!! position precition 3mm ~ 0.01 deg
  TH2F *hSrcPos_on = new TH2F ("hSrcPosOn","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);

  //
  //tasks
  //
  
  MReadTree read_on("Parameters", f_on_name);
  read_on.DisableAutoScheme();
 
  MSrcPlace srcplace;
  MHillasSrcCalc csrc_on;
  
//   // prints
//   MPrint pevent("MRawEvtHeader");
//   MPrint phillas("MHillas");
//   MPrint phillassrc("MHillasSrc");
//   MPrint psrcpos("MSrcPosCam");
  
  //tasklist
  tlist_on.AddToList(&read_on);
  tlist_on.AddToList(&srcplace);
  tlist_on.AddToList(&csrc_on);
  
  // Create and setup the eventloop
  MEvtLoop loop_on;
  loop_on.SetParList(&plist_on);
  //loop_on.SetDisplay(display);
  
//   MProgressBar bar;
//   loop_on.SetProgressBar(&bar);
  
//   if (!loop_on.Eventloop())
//     return;

  if (!loop_on.PreProcess())
    return;
  
  Double_t mjdFirstEventofBin=0;
  Double_t mjdFirstEvent=0;

  MTime mtimeFirstEventofBin;
  MTime mtimeFirstEvent;

  Double_t mjdLastEvent=0;
  MTime  mtimeLastEvent=0;
  UInt_t   evtLastEvent;
  UInt_t   runLastEvent=0;
  MTime    mtimeLastEvent;

  Bool_t flag = kFALSE;

  zenithMinimumOn[numberTimeBins-1] = 100.;
  zenithMaximumOn[numberTimeBins-1] = 0.;
  timeOn[numberTimeBins-1] = 0;

  //create histos needed in the time bins

  TString alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1);
  TH1F *hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);

  TString srcposTitle =  Form("%s%02i","hSrcPosOn",numberTimeBins-1);
  TH2F *hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);

  Int_t nbins_coszenith = 200;
  Double_t mincoszenith = 0.;  
  Double_t maxcoszenith = 1.;  
  TString coszenithTitle =  Form("%s%02i","hCosZenithOn",numberTimeBins-1);
  TH1F *hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);

  while(tlist_on.Process())
    {
      // Compute live time 
      UInt_t   run = runheader_on.GetRunNumber();
      UInt_t   evt = evtheader_on.GetDAQEvtNumber();
      Double_t mjd = time_on.GetMjd();

      if (mjd == 0)
        {
	  
	  if (debug)
	    {
	      cout << "MTime::GetTime() == 0" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();
	    }
	  
	  flag = kTRUE;
       }
      else if (mjd-mjdLastEvent < 0 && flag)
        {
	  
	  if (debug)
	    {
	      cout << "mjd-mjdLastEvent < 0" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();	
	    }
	  
	  flag = kTRUE;
       }
      else if (mjd-mjdLastEvent == 0 && flag)
        {
	  
	  if (debug)
	    {
	      cout << "mjd-mjdLastEvent == 0" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();	
	    }

	  flag = kTRUE;
       }
      else if (evt-evtLastEvent <= 0)
        {
	  
	  if (debug)
	    {
	      cout << "evt-evtLastEvent <= 0" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();	
	    }

	  flag = kTRUE;
        }

      else if ((Int_t)mjd == mjd)
        {
	  
	  if (debug)
	    {
	      cout << "(Int_t)mjd == mjd" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();	
	    }

	  flag = kTRUE;
        }
      else
	{

	  if (flag && debug)
	    {
	    
	      cout << "flag = TRUE" << endl;
	      cout.precision(15);
	      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
	      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
	      mtimeLastEvent.Print();
	      time_on.Print();	
	      
	      flag = kFALSE;

	    }

	  if ( run !=  runLastEvent )
	    {
	      if ( mjdLastEvent != 0 )
		{
		  timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
		  
		  cout.precision(10);
		  cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
		  if (flag && debug)
		    {
		      cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl;
		      cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd <<  endl;
		      mtimeLastEvent.Print();
		      time_on.Print();	
		    }
		  
		}
	      
	      mjdFirstEvent = mjd;
	    }
	  
	  if (mjdFirstEventofBin == 0)
	    {
	      mjdFirstEvent = mjd;
	      mjdFirstEventofBin = mjd;
	    }
	  
	  evtLastEvent = evt;
	  runLastEvent = run;
	  mjdLastEvent = mjd;
	  mtimeLastEvent = time_on;

      
      
	  Double_t size = hillas_on.GetSize();
	  Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg();
	  Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg();
	  Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg();
	  Double_t alpha = hillassrc_on.GetAlpha();
	  Double_t srcposx = srcpos_on.GetX();
	  Double_t srcposy = srcpos_on.GetY();
	  Double_t zenith = drive_on.GetNominalZd();
	  

	  Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
	  
	  if (sizebin >= 0)
	    {
	      if (width > widthmin[sizebin] && width < widthmax[sizebin])
		{
		  if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
		    {
		      if (dist > distmin[sizebin] && dist < distmax[sizebin])
			{      

			  //General histos
			  hSize_on->Fill(log10(size));
			  hWidth_on->Fill(width);
			  hLength_on->Fill(length);
			  hDist_on->Fill(dist);
			  hAlpha_on_abs->Fill(TMath::Abs(alpha));
			  hAlpha_on->Fill(alpha);
			  hSrcPos_on->Fill(srcposx,srcposy);
			  
			  // Time bin histos
			  hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
			  hSrcPos_on_timebin->Fill(srcposx,srcposy);
			  hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad));
			  
			  if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1])
			    zenithMinimumOn[numberTimeBins-1] = zenith;
			  if (zenith>zenithMaximumOn[numberTimeBins-1])
			    zenithMaximumOn[numberTimeBins-1] = zenith;
			  
			}
		    }
		}
	    }
	  
  
      // Check if you are overload the maxim time bin
        if ((mjd-mjdFirstEventofBin)*kDay >= timebin)
	  {
	      //Compute the time on
	      timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
	      widthTimeBinOn[numberTimeBins-1] = (mjdLastEvent-mjdFirstEventofBin)/2;
	      meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
	      
	      cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
	      cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " widthTimeBinOn " <<  widthTimeBinOn[numberTimeBins-1] << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << ' ' <<  mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1] << endl;
	      
	      //Compute the number of on events
	      numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
	      numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);	      
	      //Initialize variables
	      mjdFirstEvent = mjd;
	      mjdFirstEventofBin = mjd;
	      
	      alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
	      srcposHistoOn.AddLast(hSrcPos_on_timebin);
	      coszenithHistoOn.AddLast(hCosZenith_on_timebin);
	      
// 	      cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;

	      //Increase the number of time bins and all needed arrays
	      numberTimeBins++;
	      
	      timeOn.Set(numberTimeBins);
	      numberOnEvents.Set(numberTimeBins);
	      numberBkgEventsToNormOn.Set(numberTimeBins);
	      widthTimeBinOn.Set(numberTimeBins);
	      meanTimeBinOn.Set(numberTimeBins);
	      zenithMinimumOn.Set(numberTimeBins);
	      zenithMaximumOn.Set(numberTimeBins);
	      
	      timeOn[numberTimeBins-1] = 0.;
	      zenithMinimumOn[numberTimeBins-1] = 100.;
	      zenithMaximumOn[numberTimeBins-1] = 0.;
	      
	      alphaTitle =  Form("%s%02i","hAlphaOn",numberTimeBins-1);         
	      hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);
	      
	      srcposTitle =  Form("%s%02i","hSrcPosOn",numberTimeBins-1);
	      hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);

	      coszenithTitle =  Form("%s%02i","hCosZenithOn",numberTimeBins-1);
	      hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith);

	      flag = kTRUE;
	  }
	}	
    }

  
  //Compute the time on for last time bin
  timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay;
  widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2;
  meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1];
  //Compute the number of on events for the last time bin
  numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
  numberBkgEventsToNormOn[numberTimeBins-1] =  (Double_t)  hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);	      

  cout.precision(10);
  cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl;
  cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << " widthTimeBinOn " <<  widthTimeBinOn[numberTimeBins-1] << endl;

  alphaHistoOn.AddLast(hAlpha_on_abs_timebin);
  srcposHistoOn.AddLast(hSrcPos_on_timebin);
  coszenithHistoOn.AddLast(hCosZenith_on_timebin);

  //  cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin <<  " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;

  loop_on.PostProcess();
  
  tlist_on.PrintStatistics();
  
  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    cout << bin << " timeOn " << timeOn[bin] << " min-max zenithOn " << zenithMinimumOn[bin] << "-" <<  zenithMaximumOn[bin] << " min-max cos(zenithOn) " << TMath::Cos(zenithMinimumOn[bin]*kConvDegToRad) << "-" <<  TMath::Cos(zenithMaximumOn[bin]*kConvDegToRad) << " numberOnEvents " << numberOnEvents[bin] << endl;


  //-----------------------OFF------------------------

  TObjArray alphaHistoOff;
  TObjArray srcposHistoOff;

  TArrayD timeOff(numberTimeBins);

  TArrayI numberOffEvents(numberTimeBins);

  TArrayD numberBkgEventsToNormOff(numberTimeBins);
  TArrayD onOffNormFactor(numberTimeBins);

  TArrayD numberExcessEvents(numberTimeBins);
  TArrayD errorNumberExcessEvents(numberTimeBins);

  TArrayD numberExcessEventsPerSec(numberTimeBins);
  TArrayD errorNumberExcessEventsPerSec(numberTimeBins);

  TArrayD numberExcessEventsPerMin(numberTimeBins);
  TArrayD errorNumberExcessEventsPerMin(numberTimeBins);

  timeOff.Set(numberTimeBins);
  
  TH1F* hAlpha_off_abs_timebin;
  TH2F* hSrcPos_off_timebin;
  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    {
      alphaTitle =  Form("%s%02i","hAlphaOff",bin);
      hAlpha_off_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs);;
      alphaHistoOff.AddLast(hAlpha_off_abs_timebin);
      //      cout << "hAlpha_off_abs_timebin " << hAlpha_off_abs_timebin <<  " alphaHistoOff [" << bin << "] " << alphaHistoOff[bin] << endl;

      srcposTitle =  Form("%s%02i","hSrcPosOff",bin);
      hSrcPos_off_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);
      srcposHistoOff.AddLast(hSrcPos_off_timebin);
      //      cout << "hSrcPos_off_timebin " << hSrcPos_off_timebin <<  " srcposHistoOff [" << bin << "] " << srcposHistoOff[bin] << endl;
    }

  // 
  // Make a loop only for the OFF data:
  //
  
  MParList plist_off;
  MTaskList tlist_off;
  plist_off.AddToList(&tlist_off);
  
  MRawRunHeader runheader_off;
  MRawEvtHeader evtheader_off;
  MTime              time_off;
  MHillas          hillas_off;
  MHillasSrc    hillassrc_off;
  MSrcPosCam       srcpos_off;
  MReportDrive      drive_off;
  
  plist_off.AddToList(&geomcam);
  plist_off.AddToList(&runheader_off);
  plist_off.AddToList(&evtheader_off);
  plist_off.AddToList(&time_off);
  plist_off.AddToList(&hillas_off);
  plist_off.AddToList(&hillassrc_off);
  plist_off.AddToList(&srcpos_off);
  plist_off.AddToList(&drive_off);

  TH1F *hSize_off      = new TH1F ("hSizeOff","",nbins_Size,min_Size,max_Size);
  TH1F *hWidth_off     = new TH1F ("hWidthOff","",nbins_Width,min_Width,max_Width);
  TH1F *hLength_off    = new TH1F ("hLengthOff","",nbins_Length,min_Length,max_Length);
  TH1F *hDist_off      = new TH1F ("hDistOff","",nbins_Dist,min_Dist,max_Dist);
  TH1F *hAlpha_off_abs = new TH1F ("hAbsAlphaOff","",nbins_abs,minalpha_abs,maxalpha_abs);
  TH1F *hAlpha_off     = new TH1F ("hAlphaOff","",nbins,minalpha,maxalpha);
  TH2F *hSrcPos_off    = new TH2F ("hSrcPosOff","",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos);

  //tasks
  MReadTree read_off("Parameters", f_off_name);
  read_off.DisableAutoScheme();

  // taks to set the src position in the off data for the time bins
  // --------------------------------------------------------------
  MSrcPlace srcplace_timebin;
  srcplace_timebin.SetCreateHisto(kFALSE);
  srcplace_timebin.SetMode(MSrcPlace::kOff);
  // --------------------------------------------------------------

  //  srcplace.SetMode(MSrcPlace::kOff);
  
  MHillasSrcCalc csrc_off;
  
  //tasklist
  tlist_off.AddToList(&read_off);
  tlist_off.AddToList(&srcplace_timebin); // This is just to run the preprocess correctely
  tlist_off.AddToList(&csrc_off); // This is just to run the preprocess correctely
  
  // Create and setup the eventloop
  MEvtLoop loop_off;
  loop_off.SetParList(&plist_off);
  //loop_off.SetDisplay(display);
  
//   MProgressBar bar_off;
//   loop_off.SetProgressBar(&bar_off);
  
//   if (!loop_off.Eventloop(numEntries))
//     return;
  
  if (!loop_off.PreProcess())
    return;

  while(tlist_off.Process())
    {

      //First read the variables source(i.e. time bin) independent
      Double_t zenith = drive_off.GetNominalZd();

      Double_t size = hillas_off.GetSize();
      Double_t width = hillas_off.GetWidth()*geomcam.GetConvMm2Deg();
      Double_t length = hillas_off.GetLength()*geomcam.GetConvMm2Deg();
      
      Int_t sizebin = GetBin(size,numberSizeBins,sizeBins);
      
      if (sizebin >= 0)
	{
	  if (width > widthmin[sizebin] && width < widthmax[sizebin])
	    {
	      if (length > lengthmin[sizebin] && length < lengthmax[sizebin])
		{
		  //General histos

		  if (dist > distmin[sizebin] && dist < distmax[sizebin])
		    {
		      hSize_off->Fill(log10(size));
		      hWidth_off->Fill(width);
		      hLength_off->Fill(length);

		      srcplace_timebin.SetPositionHisto(hSrcPos_on);
		      srcplace_timebin.CallProcess();
		      csrc_off.CallProcess();
		  
		      Double_t dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
		      Double_t alpha = hillassrc_off.GetAlpha();
		      Double_t srcposx = srcpos_off.GetX();
		      Double_t srcposy = srcpos_off.GetY();

		      hDist_off->Fill(dist);
		      hAlpha_off_abs->Fill(TMath::Abs(alpha));
		      hAlpha_off->Fill(alpha);
		      hSrcPos_off->Fill(srcposx,srcposy);
		    }

		  // Time bin histos
		  for (UInt_t bin=0; bin<numberTimeBins; bin++)
		    {
		      srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
		      srcplace_timebin.CallProcess();
		      csrc_off.CallProcess();
		      
		      Double_t dist = hillassrc_off.GetDist()*geomcam.GetConvMm2Deg();
		      Double_t alpha = hillassrc_off.GetAlpha();
		      Double_t srcposx = srcpos_off.GetX();
		      Double_t srcposy = srcpos_off.GetY();
		      
		      if (dist > distmin[sizebin] && dist < distmax[sizebin])
			{
			  ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
			  ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
			}
		    }
		}
	    }
	}
      
      //       if (zenith >  zenithMaximumOn[numberTimeBins-1])
      // 	{
      // 	  cout << "Exit loop because we arrive to zenith = " << zenith << endl;
      // 	  break;
      // 	}
    }
  
  loop_off.PostProcess();

  tlist_off.PrintStatistics();

  TArrayD meanTimeBinOnInSec(numberTimeBins);
  TArrayD widthTimeBinOnInSec(numberTimeBins);
  
  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    {
      
      meanTimeBinOnInSec[bin]  = (meanTimeBinOn[bin]-(Int_t)meanTimeBinOn[bin]);
      if (meanTimeBinOnInSec[bin] > 0.5)
	meanTimeBinOnInSec[bin] = meanTimeBinOnInSec[bin] - 1 ;
      meanTimeBinOnInSec[bin] *= kDay;

      widthTimeBinOnInSec[bin] = widthTimeBinOn[bin]*kDay;

      numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
      numberBkgEventsToNormOff[bin] =  (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
      if (numberOffEvents[bin] != 0 && numberBkgEventsToNormOff[bin] != 0)
	{
	  onOffNormFactor[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin];
	  numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin];
	  errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]);
	  numberExcessEventsPerSec[bin] = numberExcessEvents[bin]/timeOn[bin];
	  errorNumberExcessEventsPerSec[bin] = errorNumberExcessEvents[bin]/timeOn[bin];
	  numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin];
	  errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin];
	}
    }
  
  for (UInt_t bin=0; bin<numberTimeBins; bin++)
    {
      cout.precision(5);
      cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl;
      cout << " numberOnEvents\t " << numberOnEvents[bin] << endl;
      cout << " numberOffEvents\t " << numberOffEvents[bin] << endl;
      cout << " numberBkgEventsToNormOn\t " << numberBkgEventsToNormOn[bin] << endl;
      cout << " numberBkgEventsToNormOff\t " << numberBkgEventsToNormOff[bin] << endl;
      cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl;
      cout << " numberExcessEvents\t " <<  numberExcessEvents[bin] <<  " +- " << errorNumberExcessEvents[bin] << endl;
      cout << " Excess/Sec\t " << numberExcessEventsPerSec[bin] << " +- " << errorNumberExcessEventsPerSec[bin] << endl;
    }

  

  TString openpsname = psname + "(";
  TString closepsname = psname + ")";

  TCanvas *c0 = new TCanvas;
  c0->cd(1);
  TGraphErrors* lightcurvegraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),numberExcessEventsPerMin.GetArray(),widthTimeBinOnInSec.GetArray(),errorNumberExcessEventsPerMin.GetArray());
  lightcurvegraph->SetTitle("LightCurve");
  lightcurvegraph->SetMinimum(0.);
  lightcurvegraph->SetMarkerStyle(21);
  lightcurvegraph->SetMarkerSize(0.03);
  lightcurvegraph->Draw("AP");
  lightcurvegraph->GetYaxis()->SetTitle("Excess/min");
  lightcurvegraph->GetXaxis()->SetTitle("UTC Time");
  lightcurvegraph->GetXaxis()->SetTimeDisplay(1);
  c0->Print(openpsname);
      
  TCanvas** c = new TCanvas*[numberTimeBins];
  
  for (UInt_t bin=0; bin<numberTimeBins-1; bin++)
    {
      c[bin] = new TCanvas;
      c[bin]->cd(1);
      
      ((TH1F*)alphaHistoOn[bin])->Sumw2();
      ((TH1F*)alphaHistoOff[bin])->SetStats(0);
      ((TH1F*)alphaHistoOff[bin])->Sumw2();
      ((TH1F*)alphaHistoOff[bin])->Scale(onOffNormFactor[bin]); 
      ((TH1F*)alphaHistoOn[bin])->SetStats(0); //-> Do NOT show the legend with statistics
      ((TH1F*)alphaHistoOn[bin])->SetLineColor(kBlack);
      ((TH1F*)alphaHistoOn[bin])->SetMarkerStyle(21);
      ((TH1F*)alphaHistoOn[bin])->SetMarkerColor(kBlack);
      ((TH1F*)alphaHistoOn[bin])->SetMarkerSize(0.7);
      ((TH1F*)alphaHistoOff[bin])->SetFillColor(46);
      ((TH1F*)alphaHistoOff[bin])->SetLineColor(46);
      ((TH1F*)alphaHistoOff[bin])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
      ((TH1F*)alphaHistoOff[bin])->SetMinimum(0.);
      alphaTitle =  Form("%s%02i","hAlphaOnOff",bin);
      ((TH1F*)alphaHistoOn[bin])->SetTitle(alphaTitle);
      
      
      ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1P");
      ((TH1F*)alphaHistoOff[bin])->DrawCopy("HISTSAME");
      ((TH1F*)alphaHistoOff[bin])->DrawCopy("ESAME");
      ((TH1F*)alphaHistoOn[bin])->DrawCopy("E1PSAME");

      c[bin]->Print(psname);
    }

  c[numberTimeBins-1] = new TCanvas;
  c[numberTimeBins-1]->cd(1);
  
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->Sumw2();
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetStats(0);
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->Sumw2();
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->Scale(onOffNormFactor[numberTimeBins-1]); 
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetStats(0); //-> Do NOT show the legend with statistics
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetLineColor(kBlack);
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerStyle(21);
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerColor(kBlack);
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetMarkerSize(0.7);
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillColor(46);
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetLineColor(46);
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->SetMinimum(0.);
  alphaTitle =  Form("%s%02i","hAlphaOnOff",numberTimeBins-1);
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->SetTitle(alphaTitle);
  
  
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1P");
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("HISTSAME");
  ((TH1F*)alphaHistoOff[numberTimeBins-1])->DrawCopy("ESAME");
  ((TH1F*)alphaHistoOn[numberTimeBins-1])->DrawCopy("E1PSAME");
  
  c[numberTimeBins-1]->Print(psname);

//   TFile input("./prueba.root", "RECREATE");

//    for (UInt_t bin=0; bin<numberTimeBins; bin++)
//      {
//        ((TH1F*)alphaHistoOn[bin])->Write();
//        ((TH1F*)alphaHistoOff[bin])->Write();
//        ((TH2F*)srcposHistoOn[bin])->Write();
//        ((TH2F*)srcposHistoOff[bin])->Write();
//        ((TH1F*)coszenithHistoOn[bin])->Write();
//      }

//    input.Close();

  // ############################################################################
  // Calculate significance and excess: 
  // ############################################################################

  Double_t norm_on_abs  = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
  Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
  Double_t norm_off_abs  = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
  Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
  Double_t norm = norm_on_abs/norm_off_abs;

  char text_tit_alpha[256];
  sprintf(text_tit_alpha, " Alpha Plot On and Off ");
  hAlpha_off_abs->SetTitle(text_tit_alpha);
  hAlpha_on_abs->SetTitle(text_tit_alpha);

  Double_t excess  = exces_on_abs - exces_off_abs*norm;
  Double_t sign    = excess / sqrt( exces_on_abs + norm*norm*exces_off_abs );
  Double_t int_off = (Double_t) hAlpha_off_abs->Integral(1, 18);
  int hAlpha_on_entries  = (int) hAlpha_on_abs->GetEntries();
  int hAlpha_off_entries = (int) hAlpha_off_abs->GetEntries();
    
  cout << "---> Normalization F factor =\t" << norm <<endl;
  cout << "---> Excess =\t\t\t" << excess <<endl;
  cout << "---> Significancia =\t\t" << sign <<endl;    
  cout << "---> entries on   =\t\t" << hAlpha_on_entries  <<endl;
  cout << "---> entries off  =\t\t" << hAlpha_off_entries <<endl;
  cout << "---> integral off =\t\t" << int_off <<endl;

  Double_t shiftx;

  // ############################################################################
  // Draw SIZE
  // ############################################################################
  TCanvas *c1 = new TCanvas;

  gPad->cd();

  gPad->SetLogy();
  hSize_on->Sumw2();
  hSize_off->Sumw2();
  hSize_off->Scale(norm); 
  hSize_on->SetLineColor(kBlack);
  hSize_on->SetMarkerStyle(21);
  hSize_on->SetMarkerSize(0.7);
  hSize_on->SetMarkerColor(kBlack);
  hSize_off->SetFillColor(46);
  hSize_off->SetLineColor(46);
  hSize_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hSize_off->SetMinimum(0.1);
  hSize_on->SetMinimum(0.1);
  hSize_on->SetTitle("SIZE distribution");
  hSize_off->SetTitle("SIZE distribution");

  hSize_on->DrawCopy("E1P");

  // move stat box to make them all visible
  gPad->Update();
  TPaveStats* pavs_on_size = (TPaveStats*) hSize_on->GetListOfFunctions()->FindObject("stats");
  if(pavs_on_size){
    shiftx = pavs_on_size->GetX2NDC() - pavs_on_size->GetX1NDC();
    pavs_on_size->SetX1NDC(pavs_on_size->GetX1NDC() - shiftx);
    pavs_on_size->SetX2NDC(pavs_on_size->GetX2NDC() - shiftx);  
  }
  gPad->Modified();
  gPad->Update();

  hSize_off->DrawCopy("HISTSAME");
  hSize_off->DrawCopy("ESAME");

  gPad->Modified();
  gPad->Update();

  Double_t chisize = ChiSquareNDof((TH1D*)hSize_on,(TH1D*)hSize_off);

  Double_t x_label_pos  = log10(1000000)*0.7;
  Double_t y_label_pos  = log10((hSize_on->GetBinContent(hSize_on->GetMaximumBin()))/2.);
  Double_t textsize = 0.03;

  char text_size[256];
  sprintf(text_size,"ChiSquare/NDof = %4.2f",chisize);

  TLatex *tsize = new TLatex(x_label_pos, y_label_pos, text_size);
  tsize->SetTextSize(textsize);
//  tsize->Draw();

  gPad->Modified();
  gPad->Update();

  c1->Print(psname);

  // ############################################################################
  // DrawCopy DIST
  // ############################################################################
  TCanvas *c2 = new TCanvas;

  gPad->cd();

  hDist_on->Sumw2();
  hDist_off->Sumw2();
  hDist_off->Scale(norm); 
  hDist_on->SetLineColor(kBlack);
  hDist_on->SetMarkerStyle(21);
  hDist_on->SetMarkerSize(0.7);
  hDist_on->SetMarkerColor(kBlack);
  hDist_off->SetFillColor(46);
  hDist_off->SetLineColor(46);
  hDist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hDist_off->SetMinimum(0.);
  hDist_on->SetTitle("DIST distribution");
  hDist_off->SetTitle("DIST distribution");

  hDist_on->DrawCopy("E1P");

  // move stat box to make them all visible
  gPad->Update();
  TPaveStats* pavs_on_dist = (TPaveStats*) hDist_on->GetListOfFunctions()->FindObject("stats");
  if(pavs_on_dist){
    shiftx = pavs_on_dist->GetX2NDC() - pavs_on_dist->GetX1NDC();
    pavs_on_dist->SetX1NDC(pavs_on_dist->GetX1NDC() - shiftx);
    pavs_on_dist->SetX2NDC(pavs_on_dist->GetX2NDC() - shiftx);  
  }
  gPad->Modified();
  gPad->Update();

  hDist_off->DrawCopy("HISTSAME");
  hDist_off->DrawCopy("ESAME");
  hDist_on->DrawCopy("E1PSAME");

  Double_t chidist = ChiSquareNDof((TH1D*)hDist_on,(TH1D*)hDist_off);

  x_label_pos  = distmax[numberSizeBins-1]*0.7;
  y_label_pos  = hDist_on->GetBinContent(hDist_on->GetMaximumBin())/2.;

  char text_dist[256];
  sprintf(text_size,"ChiSquare/NDof = %4.2f",chidist);

  TLatex *tdist = new TLatex(x_label_pos, y_label_pos, text_dist);
  tdist->SetTextSize(textsize);
//  tdist->Draw();

  gPad->Modified();
  gPad->Update();

  c2->Print(psname);

   // ############################################################################
  // DrawCopy WIDTH
  // ############################################################################
  TCanvas *c3 = new TCanvas;

  gPad->cd();

  hWidth_off->Sumw2();
  hWidth_off->Scale(norm); 
  hWidth_on->SetLineColor(kBlack);
  hWidth_on->SetMarkerStyle(21);
  hWidth_on->SetMarkerSize(0.7);
  hWidth_on->SetMarkerColor(kBlack);
  hWidth_off->SetFillColor(46);
  hWidth_off->SetLineColor(46);
  hWidth_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hWidth_off->SetMinimum(0.);
  hWidth_on->SetTitle("WIDTH distribution");
  hWidth_off->SetTitle("WIDTH distribution");

  hWidth_on->DrawCopy("E1P");

  // move stat box to make them all visible
  gPad->Update();
  TPaveStats* pavs_on_width = (TPaveStats*) hWidth_on->GetListOfFunctions()->FindObject("stats");
  if(pavs_on_width){
    shiftx = pavs_on_width->GetX2NDC() - pavs_on_width->GetX1NDC();
    pavs_on_width->SetX1NDC(pavs_on_width->GetX1NDC() - shiftx);
    pavs_on_width->SetX2NDC(pavs_on_width->GetX2NDC() - shiftx);  
  }
  gPad->Modified();
  gPad->Update();

  hWidth_off->DrawCopy("HISTSAME");
  hWidth_off->DrawCopy("ESAME");
  hWidth_on->DrawCopy("E1PSAME");

  Double_t chiwidth = ChiSquareNDof((TH1D*)hWidth_on,(TH1D*)hWidth_off);

  x_label_pos  = widthmax[numberSizeBins-1]*0.7;
  y_label_pos  = hWidth_on->GetBinContent(hWidth_on->GetMaximumBin())/2.;

  char text_width[256];
  sprintf(text_size,"ChiSquare/NDof = %4.2f",chiwidth);

  TLatex *twidth = new TLatex(x_label_pos, y_label_pos, text_width);
  twidth->SetTextSize(textsize);
//  twidth->Draw();

  gPad->Modified();
  gPad->Update();
 
  c3->Print(psname);

  // ############################################################################
  // DrawCopy LENGTH
  // ############################################################################
  TCanvas *c4 = new TCanvas;
 
  gPad->cd();

  hLength_on->Sumw2();
  hLength_off->Sumw2();
  hLength_off->Scale(norm); 
  hLength_on->SetLineColor(kBlack);
  hLength_on->SetMarkerStyle(21);
  hLength_on->SetMarkerSize(0.7);
  hLength_on->SetMarkerColor(kBlack);
  hLength_off->SetFillColor(46);
  hLength_off->SetLineColor(46);
  hLength_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hLength_off->SetMinimum(0.);
  hLength_on->SetTitle("LENGTH distribution");
  hLength_off->SetTitle("LENGTH distribution");

  hLength_on->DrawCopy("E1P");

  // move stat box to make them all visible
  gPad->Update();
  TPaveStats* pavs_on_length = (TPaveStats*) hLength_on->GetListOfFunctions()->FindObject("stats");
  if(pavs_on_length){
    shiftx = pavs_on_length->GetX2NDC() - pavs_on_length->GetX1NDC();
    pavs_on_length->SetX1NDC(pavs_on_length->GetX1NDC() - shiftx);
    pavs_on_length->SetX2NDC(pavs_on_length->GetX2NDC() - shiftx);  
  }
  gPad->Modified();
  gPad->Update();

  hLength_off->DrawCopy("HISTSAME");
  hLength_off->DrawCopy("ESAME");
  hLength_on->DrawCopy("E1PSAME");

  Double_t chilength = ChiSquareNDof((TH1D*)hLength_on,(TH1D*)hLength_off);

  x_label_pos  = lengthmax[numberSizeBins-1]*0.7;
  y_label_pos  = hLength_on->GetBinContent(hLength_on->GetMaximumBin())/2.;

  char text_length[256];
  sprintf(text_size,"ChiSquare/NDof = %4.2f",chilength);

  TLatex *tlength = new TLatex(x_label_pos, y_label_pos, text_length);
  tlength->SetTextSize(textsize);
//  tlength->Draw();

  gPad->Modified();
  gPad->Update();

  c4->Print(psname);

 // ############################################################################
  // DrawCopy normalized ALPHA plot
  // ############################################################################
  TCanvas *c5 = new TCanvas;
  
  gPad->cd();

  hAlpha_on_abs->Sumw2();
  hAlpha_off_abs->SetStats(0);
  hAlpha_off_abs->Sumw2();
  hAlpha_off_abs->Scale(norm); 
  hAlpha_on_abs->SetStats(0); //-> Do NOT show the legend with statistics
  hAlpha_on_abs->SetLineColor(kBlack);
  hAlpha_on_abs->SetMarkerStyle(21);
  //hAlpha_on_abs->SetMarkerSize();
  hAlpha_on_abs->SetMarkerColor(kBlack);
  hAlpha_on_abs->SetMarkerSize(0.7);
  hAlpha_off_abs->SetFillColor(46);
  hAlpha_off_abs->SetLineColor(46);
  hAlpha_off_abs->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hAlpha_off_abs->SetMinimum(0.);
  hAlpha_on_abs->SetTitle("Alpha plot");
  hAlpha_off_abs->SetTitle("Alpha plot");

  
  hAlpha_on_abs->DrawCopy("E1P");
  hAlpha_off_abs->DrawCopy("HISTSAME");
  hAlpha_off_abs->DrawCopy("ESAME");
  hAlpha_on_abs->DrawCopy("E1PSAME");


   //draw the LEGEND with excess and significance values in the alpha plot:
  char text_Fnorm[256], text_excess[256], text_sign[256];
  char text_entries_on[256], text_entries_off[256], text_integral_off[256];
  int hAlpha_on_entries  = (int) hAlpha_on_abs->GetEntries();
  int hAlpha_off_entries = (int) hAlpha_off_abs->GetEntries();
  sprintf(text_Fnorm,       " F norm =       %.3f", norm);
  sprintf(text_excess,      " Excess =       %.3f", excess);
  sprintf(text_sign,        " Significance = %.3f", sign);
  sprintf(text_entries_on,  " Entries ON   = %d",  hAlpha_on_entries);
  sprintf(text_entries_off, " Entries OFF  = %d",  hAlpha_off_entries);
  sprintf(text_integral_off," Integral OFF = %d",  int_off);
  
  x_label_pos  = 90.*0.7;
  y_label_pos  = (hAlpha_on_abs->GetBinContent(hAlpha_on_abs->GetMaximumBin())); //2.;
  Double_t y_label_step = y_label_pos / 8.;

  TLatex *t0 = new TLatex(x_label_pos, y_label_pos - y_label_step*0, text_Fnorm);
  t0->SetTextSize(textsize);
  t0->Draw();
  TLatex *t1 = new TLatex(x_label_pos, y_label_pos - y_label_step*1, text_excess);
  t1->SetTextSize(textsize);
  t1->Draw();
  TLatex *t2 = new TLatex(x_label_pos, y_label_pos - y_label_step*2, text_sign);
  t2->SetTextSize(textsize);
  t2->Draw();
  TLatex *t3 = new TLatex(x_label_pos, y_label_pos - y_label_step*3, text_entries_on);
  t3->SetTextSize(textsize);
  t3->Draw();
  TLatex *t4 = new TLatex(x_label_pos, y_label_pos - y_label_step*4, text_entries_off);
  t4->SetTextSize(textsize);
  t4->Draw();
  TLatex *t5 = new TLatex(x_label_pos, y_label_pos - y_label_step*5, text_integral_off);
  t5->SetTextSize(textsize);
  t5->Draw();
  

  Double_t chialpha = ChiSquareNDof((TH1D*)hAlpha_on_abs,(TH1D*)hAlpha_off_abs);

  y_label_pos  = (hAlpha_on_abs->GetBinContent(hAlpha_on_abs->GetMaximumBin()))/2.;

  char text_alpha[256];
  sprintf(text_size,"ChiSquare/NDof = %4.2f",chialpha);

  TLatex *talpha = new TLatex(x_label_pos, y_label_pos, text_alpha);
  talpha->SetTextSize(textsize);
//  talpha->Draw();

  gPad->Modified();
  gPad->Update();

  c5->Print(psname);

  // ############################################################################
  // DrawCopy normalized alpha histos for alpha form -90 to 90 deg.
  // ############################################################################
  TCanvas *c6 = new TCanvas;

  gPad->cd();

  hAlpha_on->Sumw2();
  hAlpha_off->SetStats(0);
  hAlpha_off->Sumw2();
  hAlpha_off->Scale(norm); 
  hAlpha_off->SetFillColor(46);
  hAlpha_off->SetLineColor(46);
  hAlpha_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
  hAlpha_off->SetMinimum(0.); 
  hAlpha_on->SetStats(0); //-> Do NOT show the legend with statistics
  hAlpha_on->SetLineColor(kBlack);
  hAlpha_on->SetMarkerStyle(21);
  hAlpha_on->SetMarkerSize(0.7);
  hAlpha_on->SetMarkerColor(kBlack);
  hAlpha_on->SetTitle("Alpha plot form -90 to 90 deg");
  hAlpha_off->SetTitle("Alpha plot form -90 to 90 deg");

  hAlpha_on->DrawCopy("E1P");
  hAlpha_off->DrawCopy("HISTSAME");
  hAlpha_off->DrawCopy("ESAME");
  hAlpha_on->DrawCopy("E1PSAME");

  Double_t chialpha90 = ChiSquareNDof((TH1D*)hAlpha_on,(TH1D*)hAlpha_off);

  x_label_pos  = 90.*0.5;
  y_label_pos  = hAlpha_on->GetBinContent(hAlpha_on->GetMaximumBin())/2.;

  char text_alpha90[256];
  sprintf(text_alpha90,"ChiSquare/NDof = %4.2f",chialpha90);

  TLatex *talpha90 = new TLatex(x_label_pos, y_label_pos, text_alpha90);
  talpha90->SetTextSize(textsize);
//  talpha90->Draw();

  gPad->Update();
  gPad->Modified();

  c6->Print(psname);

  cout << "---> ChiSquare/NDof [Size] =\t\t" << chisize << endl;
  cout << "---> ChiSquare/NDof [Dist] =\t\t" << chidist << endl;
  cout << "---> ChiSquare/NDof [Width] =\t\t" << chiwidth << endl;
  cout << "---> ChiSquare/NDof [Length] =\t\t" << chilength << endl;
  cout << "---> ChiSquare/NDof [Abs(Alpha)] =\t" << chialpha << endl;
  cout << "---> ChiSquare/NDof [Alpha] =\t\t" << chialpha90 << endl;


  TCanvas *c7 = new TCanvas;
  hSrcPos_on->DrawCopy("BOX");
  c7->Print(psname);

  TCanvas *c8 = new TCanvas;
  hSrcPos_off->DrawCopy("BOX");
  c8->Print(closepsname);

  cout << "Done!!" <<endl;
  
}



