Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C	(revision 4375)
+++ trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C	(revision 4375)
@@ -0,0 +1,1431 @@
+Bool_t debug = kFALSE;
+
+
+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 chiq/nbinsnozero;
+}
+
+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")
+{
+
+  //Constanst
+  const Double_t ConvDegToRad = TMath::ACos(-1)/180.;
+
+  // Time binning
+  const Double_t timebin = 500.; //[sec]
+  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;
+  
+
+  //none cuts
+//   Double_t sizemin   = 2000.; //[Photons]
+//   Double_t sizemax   = 10000000000.; //[Photons]
+//   Double_t widthmin  = 0.; 
+//   Double_t widthmax  = 2.0;
+//   Double_t lengthmin = 0.; 
+//   Double_t lengthmax = 2.0;
+//   Double_t distmin   = 0.; 
+//   Double_t distmax   = 4.0;
+//   Double_t alphamin   = 0.; 
+//   Double_t alphamax   = 90.;
+
+  
+  const Int_t numberSizeBins = 4;
+  Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons]
+
+  //cuts
+  Double_t sizemin = 2000.; //[Photons]
+  Double_t sizemax = 10E9.; //[Photons]
+
+  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 = 30.;
+  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;
+  MSrcPosCam source_on;
+  MHillas hillas;
+  MHillasSrc hillasscr;
+  
+  plist_on.AddToList(&geomcam);
+  plist_on.AddToList(&source_on);
+  plist_on.AddToList(&hillas);
+  plist_on.AddToList(&hillasscr);
+  
+  //create some 1-dim histo to test only for the ON distribution of dist, width , length, size...
+  MH3 hDist_on("MHillasSrc.fDist/MGeomCam.fConvMm2Deg");
+  hDist_on.SetName("Dist_on");
+  plist_on.AddToList(&hDist_on);
+  MBinning binsDist_on("BinningDist_on");
+  Int_t nbins_Dist = 20;
+  Double_t min_Dist = 0.;
+  Double_t max_Dist = distmax[0]*1.2;
+  binsDist_on.SetEdges(nbins_Dist, min_Dist, max_Dist);
+  plist_on.AddToList(&binsDist_on);
+  
+  MH3 hWidth_on("MHillas.fWidth/MGeomCam.fConvMm2Deg");
+  hWidth_on.SetName("Width_on");
+  plist_on.AddToList(&hWidth_on);
+  MBinning binsWidth_on("BinningWidth_on");
+  Int_t nbins_Width = 20;
+  Double_t min_Width = 0.;
+  Double_t max_Width = widthmax[0]*1.2;
+  binsWidth_on.SetEdges(nbins_Width, min_Width, max_Width);
+  plist_on.AddToList(&binsWidth_on);
+  
+  MH3 hLength_on("MHillas.fLength/MGeomCam.fConvMm2Deg");
+  hLength_on.SetName("Length_on");
+  plist_on.AddToList(&hLength_on);
+  MBinning binsLength_on("BinningLength_on");
+  Int_t nbins_Length = 20;
+  Double_t min_Length = 0.;
+  Double_t max_Length =  lengthmax[0]*1.2;
+  binsLength_on.SetEdges(nbins_Length, min_Length, max_Length);
+  plist_on.AddToList(&binsLength_on);
+  
+  MH3 hSize_on("log10(MHillas.fSize)");
+  hSize_on.SetName("Size_on");
+  plist_on.AddToList(&hSize_on);
+  MBinning binsSize_on("BinningSize_on");
+  Int_t nbins_Size = 60;
+  Double_t min_Size = log10(sizemin)*0.8;
+  Double_t max_Size = log10(1000000)*1.2;
+  binsSize_on.SetEdges(nbins_Size, min_Size, max_Size);
+  plist_on.AddToList(&binsSize_on);
+    
+  //create a histo to fill the alpha values: one alpha plot form 0 to +90 deg in abs value
+  MH3 hAlpha_on_abs("abs(MHillasSrc.fAlpha)");
+  hAlpha_on_abs.SetName("Alpha_on_abs");
+  plist_on.AddToList(&hAlpha_on_abs);
+  MBinning binsAlpha_on_abs("BinningAlpha_on_abs");
+  Int_t nbins_abs = 9;
+  Double_t minalpha_abs = 0.;
+  Double_t maxalpha_abs =90.;
+  binsAlpha_on_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs);
+  plist_on.AddToList(&binsAlpha_on_abs);
+  
+  //create a histo to fill the alpha values: one alpha plot form -90 to +90 deg.
+  MH3 hAlpha_on("MHillasSrc.fAlpha");
+  hAlpha_on.SetName("Alpha_on");
+  plist_on.AddToList(&hAlpha_on);
+  MBinning binsAlpha_on("BinningAlpha_on");
+  Int_t nbins = nbins_abs*2;
+  Double_t minalpha = -90.;
+  Double_t maxalpha =  90.;
+  binsAlpha_on.SetEdges(nbins, minalpha, maxalpha);
+  plist_on.AddToList(&binsAlpha_on);
+    
+
+  //create a histo to fill the source position 
+  MH3 hSrcPos_on("MSrcPosCam.fX","MSrcPosCam.fY");
+  hSrcPos_on.SetName("SrcPos_on");
+  plist_on.AddToList(&hSrcPos_on);
+  MBinning binsSrcPos_onX("BinningSrcPos_onX");
+  MBinning binsSrcPos_onY("BinningSrcPos_onY");
+  Int_t nbins_srcpos = 400;
+  Double_t minsrcpos = -600.;
+  Double_t maxsrcpos =  600.;
+  binsSrcPos_onX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
+  binsSrcPos_onY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
+  plist_on.AddToList(&binsSrcPos_onX);
+  plist_on.AddToList(&binsSrcPos_onY);
+  
+  MH3 hDAQEvtNumber_on("MRawEvtHeader.fDAQEvtNumber");
+  hDAQEvtNumber_on.SetName("DAQEvtNumber_on");
+  plist_on.AddToList(&hDAQEvtNumber_on);
+  MBinning binsDAQEvtNumber_onX("BinningDAQEvtNumber_onX");
+  Int_t nbins_evtnum = 1000;
+  Double_t minevtnum =  0.;
+  Double_t maxevtnum =  1000.;
+  binsDAQEvtNumber_onX.SetEdges(nbins_evtnum,minevtnum,maxevtnum);
+  plist_on.AddToList(&binsDAQEvtNumber_onX);
+  
+
+  //
+  //tasks
+  //
+  
+  MReadTree read_on("Parameters", f_on_name);
+  read_on.DisableAutoScheme();
+ 
+  MSrcPlace srcplace;
+
+  //cuts
+  TString sizestr = "(MHillas.fSize < ";
+  sizestr += sizemin;
+  sizestr += ") || (";
+  sizestr += "MHillas.fSize > ";
+  sizestr += sizemax;
+  sizestr += ")";
+  MF sizefilter(sizestr);
+
+  TString widthstr = "({MHillas.fWidth/MGeomCam.fConvMm2Deg} < ";
+  widthstr += widthmin[0];
+  widthstr += ") || (";
+  widthstr += "{MHillas.fWidth/MGeomCam.fConvMm2Deg} > ";
+  widthstr += widthmax[0];
+  widthstr += ")";
+  MF widthfilter(widthstr);
+  
+  TString lengthstr = "({MHillas.fLength/MGeomCam.fConvMm2Deg} < ";
+  lengthstr += lengthmin[0];
+  lengthstr += ") || (";
+  lengthstr += "{MHillas.fLength/MGeomCam.fConvMm2Deg} > ";
+  lengthstr += lengthmax[0];
+  lengthstr += ")";
+  MF lengthfilter(lengthstr);
+  
+  TString diststr = "({MHillasSrc.fDist/MGeomCam.fConvMm2Deg} < ";
+  diststr += distmin[0];
+  diststr += ") || (";
+  diststr += "{MHillasSrc.fDist/MGeomCam.fConvMm2Deg} > ";
+  diststr += distmax[0];
+  diststr += ")";
+  MF distfilter(diststr);
+  
+  MF evenfilter("{MRawEvtHeader.fDAQEvtNumber%3}<0.5");
+  MF oddfilter("{MRawEvtHeader.fDAQEvtNumber%3}>0.5");
+  
+  MContinue cont_size(&sizefilter);
+  MContinue cont_width(&widthfilter);
+  MContinue cont_length(&lengthfilter);
+  MContinue cont_dist(&distfilter);
+  MContinue cont_even(&evenfilter);
+  MContinue cont_odd(&oddfilter);
+  
+  MHillasSrcCalc csrc_on;
+  
+  // fill all histograms
+  MFillH falpha_on_abs(&hAlpha_on_abs);
+  MFillH falpha_on(&hAlpha_on);
+  MFillH fdist_on(&hDist_on);
+  MFillH fwidth_on(&hWidth_on);
+  MFillH flength_on(&hLength_on);
+  MFillH fsize_on(&hSize_on);
+  MFillH fsrcpos_on(&hSrcPos_on);
+  MFillH fevtnum_on(&hDAQEvtNumber_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);
+  tlist_on.AddToList(&fsrcpos_on);
+  //  tlist_on.AddToList(&cont_odd);
+  tlist_on.AddToList(&cont_size);
+  tlist_on.AddToList(&cont_width);
+  tlist_on.AddToList(&cont_length);
+  tlist_on.AddToList(&cont_dist);
+  tlist_on.AddToList(&falpha_on_abs);
+  tlist_on.AddToList(&falpha_on);
+  tlist_on.AddToList(&fdist_on);
+  tlist_on.AddToList(&fwidth_on);
+  tlist_on.AddToList(&flength_on);
+  tlist_on.AddToList(&fsize_on);
+  //  tlist_on.AddToList(&fevtnum_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;
+  
+  MRawRunHeader* runheader_on = plist_on.FindObject("MRawRunHeader");
+  MRawEvtHeader* evtheader_on = plist_on.FindObject("MRawEvtHeader");
+  MTime* time_on = plist_on.FindObject("MTime");
+  MHillas* hillas_on = plist_on.FindObject("MHillas");
+  MHillasSrc* hillassrc_on = plist_on.FindObject("MHillasSrc");
+  MReportDrive* drive_on = plist_on.FindObject("MReportDrive");
+  MSrcPosCam* srcpos_on = plist_on.FindObject("MSrcPosCam");
+  
+
+  Double_t mjdFirstEventofBin=0;
+  Double_t mjdFirstEvent=0;
+
+  Double_t mjdLastEvent=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);
+  Int_t nbins_srcpos = 200;
+  Double_t minsrcpos = -600.; //[mm]
+  Double_t maxsrcpos = 600.;  //[mm]  //!!! position precition 3mm ~ 0.01 deg
+  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])
+			{      
+			  hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha));
+			  hSrcPos_on_timebin->Fill(srcposx,srcposy);
+			  hCosZenith_on_timebin->Fill(TMath::Cos(zenith*ConvDegToRad));
+			  
+			  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] = (mjd-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]*ConvDegToRad) << "-" <<  TMath::Cos(zenithMaximumOn[bin]*ConvDegToRad) << " 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);
+  
+  MSrcPosCam source_off;
+  
+  plist_off.AddToList(&geomcam);
+  plist_off.AddToList(&source_off);
+  plist_off.AddToList(&hillas);
+  plist_off.AddToList(&hillasscr);
+  
+  //create some 1-dim histo to test only for the OFF distribution of dist, width , length, size...
+  MH3 hDist_off("MHillasSrc.fDist/MGeomCam.fConvMm2Deg");
+  hDist_off.SetName("Dist_off");
+  plist_off.AddToList(&hDist_off);
+  MBinning binsDist_off("BinningDist_off");
+  binsDist_off.SetEdges(nbins_Dist, min_Dist, max_Dist);
+  plist_off.AddToList(&binsDist_off);
+  
+  MH3 hWidth_off("MHillas.fWidth/MGeomCam.fConvMm2Deg");
+  hWidth_off.SetName("Width_off");
+  plist_off.AddToList(&hWidth_off);
+  MBinning binsWidth_off("BinningWidth_off");
+  binsWidth_off.SetEdges(nbins_Width, min_Width, max_Width);
+  plist_off.AddToList(&binsWidth_off);
+  
+  MH3 hLength_off("MHillas.fLength/MGeomCam.fConvMm2Deg");
+  hLength_off.SetName("Length_off");
+  plist_off.AddToList(&hLength_off);
+  MBinning binsLength_off("BinningLength_off");
+  binsLength_off.SetEdges(nbins_Length, min_Length, max_Length);
+  plist_off.AddToList(&binsLength_off);
+  
+  MH3 hSize_off("log10(MHillas.fSize)");
+  hSize_off.SetName("Size_off");
+  plist_off.AddToList(&hSize_off);
+  MBinning binsSize_off("BinningSize_off");
+  binsSize_off.SetEdges(nbins_Size, min_Size, max_Size);
+  plist_off.AddToList(&binsSize_off);
+  
+  //create a histo to fill the alpha values: from 0 to 90 deg -> abs value
+  MH3 hAlpha_off_abs("abs(MHillasSrc.fAlpha)");
+  hAlpha_off_abs.SetName("Alpha_off_abs");
+  plist_off.AddToList(&hAlpha_off_abs);
+  MBinning binsAlpha_off_abs("BinningAlpha_off_abs");
+  binsAlpha_off_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs);
+  plist_off.AddToList(&binsAlpha_off_abs);
+  
+  //create a histo to fill the alpha values: from -90 to 90 deg
+  MH3 hAlpha_off("MHillasSrc.fAlpha");
+  hAlpha_off.SetName("Alpha_off");
+  plist_off.AddToList(&hAlpha_off);
+  MBinning binsAlpha_off("BinningAlpha_off");
+  binsAlpha_off.SetEdges(nbins, minalpha, maxalpha);
+  plist_off.AddToList(&binsAlpha_off);
+  
+  
+  MH3 hSrcPos_off("MSrcPosCam.fX","MSrcPosCam.fY");
+  hSrcPos_off.SetName("SrcPos_off");
+  plist_off.AddToList(&hSrcPos_off);
+  MBinning binsSrcPos_offX("BinningSrcPos_offX");
+  MBinning binsSrcPos_offY("BinningSrcPos_offY");
+  binsSrcPos_offX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
+  binsSrcPos_offY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos);
+  plist_off.AddToList(&binsSrcPos_offX);
+  plist_off.AddToList(&binsSrcPos_offY);
+  
+  MH3 hDAQEvtNumber_off("MRawEvtHeader.fDAQEvtNumber");
+  hDAQEvtNumber_off.SetName("DAQEvtNumber_off");
+  plist_off.AddToList(&hDAQEvtNumber_off);
+  MBinning binsDAQEvtNumber_offX("BinningDAQEvtNumber_offX");
+  Int_t nbins_evtnum = 100;
+  Double_t minevtnum =  0.;
+  Double_t maxevtnum =  100.;
+  binsDAQEvtNumber_offX.SetEdges(nbins_evtnum,minevtnum,maxevtnum);
+  plist_off.AddToList(&binsDAQEvtNumber_offX);
+  
+  //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;
+  
+  // fill all histograms
+  MFillH falpha_off_abs(&hAlpha_off_abs);
+  MFillH falpha_off(&hAlpha_off);
+  MFillH fdist_off(&hDist_off);
+  MFillH fwidth_off(&hWidth_off);
+  MFillH flength_off(&hLength_off);
+  MFillH fsize_off(&hSize_off);
+  MFillH fsrcpos_off(&hSrcPos_off);
+  MFillH fevtnum_off(&hDAQEvtNumber_off);
+  
+  //tasklist
+  tlist_off.AddToList(&read_off);
+  tlist_off.AddToList(&srcplace);
+  tlist_off.AddToList(&csrc_off);
+  tlist_off.AddToList(&fsrcpos_off);
+//   tlist_off.AddToList(&cont_even);
+//   tlist_off.AddToList(&cont_size);
+//   tlist_off.AddToList(&cont_width);
+//   tlist_off.AddToList(&cont_length);
+//   tlist_off.AddToList(&cont_dist);
+//   tlist_off.AddToList(&falpha_off_abs);
+//   tlist_off.AddToList(&falpha_off);
+//   tlist_off.AddToList(&fdist_off);
+//   tlist_off.AddToList(&fwidth_off);
+//   tlist_off.AddToList(&flength_off);
+//   tlist_off.AddToList(&fsize_off);
+//   tlist_off.AddToList(&fevtnum_off);
+  tlist_off.AddToList(&srcplace_timebin); // 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;
+
+
+  MHillas* hillas_off = plist_off.FindObject("MHillas");
+  MHillasSrc* hillassrc_off = plist_off.FindObject("MHillasSrc");
+  MReportDrive* drive_off = plist_off.FindObject("MReportDrive");
+  MSrcPosCam* srcpos_off = plist_off.FindObject("MSrcPosCam");
+
+  TArrayI sizecut(numberTimeBins);
+  TArrayI widthcut(numberTimeBins);
+  TArrayI lengthcut(numberTimeBins);
+  TArrayI distcut(numberTimeBins);
+  TArrayI alphacut(numberTimeBins);
+  
+  sizecut.Reset(0);
+  widthcut.Reset(0);
+  lengthcut.Reset(0);
+  distcut.Reset(0);
+  alphacut.Reset(0);
+
+  while(tlist_off.Process())
+    {
+
+      //Now fill just the alpha plot in the time/zenith bins defined in on
+      Double_t zenith = drive_off->GetNominalZd();
+      
+      for (UInt_t bin=0; bin<numberTimeBins; bin++)
+	{
+ 	  srcplace_timebin.SetPositionHisto((TH2F*)srcposHistoOn[bin]);
+ 	  srcplace_timebin.CallProcess();
+ 	  csrc_off.CallProcess();
+	  
+	  Double_t size = hillas_off->GetSize();
+	  Double_t width = hillas_off->GetWidth()*geomcam.GetConvMm2Deg();
+	  Double_t length = hillas_off->GetLength()*geomcam.GetConvMm2Deg();
+	  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();
+	  
+	  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])
+			{
+			  ((TH1F*)alphaHistoOff[bin])->Fill(TMath::Abs(alpha));
+			  ((TH2F*)srcposHistoOff[bin])->Fill(srcposx,srcposy);
+			}
+		      else
+			distcut[bin]++;	    
+		    }
+		  else	
+		    lengthcut[bin]++;
+		}
+	      else
+		widthcut[bin]++;
+	    }
+	  else
+	    sizecut[bin]++;
+	}
+      
+//       if (zenith >  zenithMaximumOn[numberTimeBins-1])
+// 	{
+// 	  cout << "Exit loop because we arrive to zenith = " << zenith << endl;
+// 	  break;
+// 	}
+    }
+  
+  loop_off.PostProcess();
+
+  tlist_off.PrintStatistics();
+
+//   for (UInt_t bin=0; bin<numberTimeBins; bin++)
+//     {
+//       cout <<  bin << ' '  << ((TH1F*)alphaHistoOff[bin])->GetEntries() << ' ' << ((TH2F*)srcposHistoOff[bin])->GetEntries() << endl;
+//       cout <<  "sizecut " << sizecut[bin] << endl;
+//       cout <<  "width   " << widthcut[bin] << endl;
+//       cout <<  "length  " << lengthcut[bin] << endl;
+//       cout <<  "dist    " << distcut[bin] << endl;
+//       cout <<  "alpha   " << alphacut[bin] << endl;
+//     }
+  
+  for (UInt_t bin=0; bin<numberTimeBins; bin++)
+    {
+      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 psname = "./prueba.ps";
+  TString openpsname = psname + "(";
+  TString closepsname = psname + ")";
+
+  TCanvas *c1 = new TCanvas;
+  c1->cd(1);
+  TGraphErrors* lightcurvegraph = new TGraphErrors(numberTimeBins,meanTimeBinOn.GetArray(),numberExcessEventsPerMin.GetArray(),widthTimeBinOn.GetArray(),errorNumberExcessEventsPerMin.GetArray());
+  lightcurvegraph->SetTitle("LightCurve");
+  lightcurvegraph->SetMinimum(0.);
+  lightcurvegraph->SetMarkerStyle(21);
+  lightcurvegraph->SetMarkerSize(0.03);
+  lightcurvegraph->Draw("AP");
+  c1->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(closepsname);
+
+  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();
+
+  if(kFALSE)
+    {
+
+  // ############################################################################
+  //  look for the histograms
+  // ############################################################################
+  
+  TH1F *hist_size_on = (TH1F*)hSize_on.GetHist();
+  TH1F *hist_size_off = (TH1F*)hSize_off.GetHist();
+  
+  TH1F *hist_dist_on = (TH1F*)hDist_on.GetHist();
+  TH1F *hist_dist_off = (TH1F*)hDist_off.GetHist();
+
+  TH1F *hist_width_on = (TH1F*)hWidth_on.GetHist();
+  TH1F *hist_width_off = (TH1F*)hWidth_off.GetHist();
+
+  TH1F *hist_length_on = (TH1F*)hLength_on.GetHist();
+  TH1F *hist_length_off = (TH1F*)hLength_off.GetHist();
+
+  TH1F *hist_on_abs = (TH1F*)hAlpha_on_abs.GetHist();
+  TH1F *hist_off_abs = (TH1F*)hAlpha_off_abs.GetHist();
+
+  TH1F *hist_on = (TH1F*)hAlpha_on.GetHist();
+  TH1F *hist_off = (TH1F*)hAlpha_off.GetHist();
+
+
+  // ############################################################################
+  // Calculate significance and excess: 
+  // ############################################################################
+
+  Double_t norm_on_abs  = (Double_t) hist_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
+  Double_t exces_on_abs = (Double_t) hist_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1);
+  Double_t norm_off_abs  = (Double_t) hist_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1);
+  Double_t exces_off_abs = (Double_t) hist_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 ");
+  hist_off_abs->SetTitle(text_tit_alpha);
+  hist_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) hist_off_abs->Integral(1, 18);
+  int hist_on_entries  = (int) hist_on_abs->GetEntries();
+  int hist_off_entries = (int) hist_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" << hist_on_entries  <<endl;
+  cout << "---> entries off  =\t\t" << hist_off_entries <<endl;
+  cout << "---> integral off =\t\t" << int_off <<endl;
+
+  Double_t shiftx;
+
+  //
+  //Create the display -> from now on, all histos are plotted
+  MStatusDisplay *display = new MStatusDisplay;
+  display->SetUpdateTime(3000);
+  display->Resize(850,700);
+  
+  // ############################################################################
+  // Draw SIZE
+  // ############################################################################
+  display->AddTab("SIZE");
+
+  gPad->cd();
+
+  gPad->SetLogy();
+  hist_size_on->Sumw2();
+  hist_size_off->Sumw2();
+  hist_size_off->Scale(norm); 
+  hist_size_on->SetLineColor(kBlack);
+  hist_size_on->SetMarkerStyle(21);
+  hist_size_on->SetMarkerSize(0.7);
+  hist_size_on->SetMarkerColor(kBlack);
+  hist_size_off->SetFillColor(46);
+  hist_size_off->SetLineColor(46);
+  hist_size_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_size_off->SetMinimum(0.1);
+  hist_size_on->SetMinimum(0.1);
+  hist_size_on->SetTitle("SIZE distribution");
+  hist_size_off->SetTitle("SIZE distribution");
+
+  hist_size_on->DrawCopy("E1P");
+
+  // move stat box to make them all visible
+  gPad->Update();
+  TPaveStats* pavs_on_size = (TPaveStats*) hist_size_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();
+
+  hist_size_off->DrawCopy("HISTSAME");
+  hist_size_off->DrawCopy("ESAME");
+
+  gPad->Modified();
+  gPad->Update();
+
+  Double_t chisize = ChiSquareNDof((TH1D*)hist_size_on,(TH1D*)hist_size_off);
+
+  Double_t x_label_pos  = log10(1000000)*0.7;
+  Double_t y_label_pos  = log10((hist_size_on->GetBinContent(hist_size_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();
+
+  // ############################################################################
+  // DrawCopy DIST
+  // ############################################################################
+  display->AddTab("DIST");
+
+  gPad->cd();
+
+  hist_dist_on->Sumw2();
+  hist_dist_off->Sumw2();
+  hist_dist_off->Scale(norm); 
+  hist_dist_on->SetLineColor(kBlack);
+  hist_dist_on->SetMarkerStyle(21);
+  hist_dist_on->SetMarkerSize(0.7);
+  hist_dist_on->SetMarkerColor(kBlack);
+  hist_dist_off->SetFillColor(46);
+  hist_dist_off->SetLineColor(46);
+  hist_dist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_dist_off->SetMinimum(0.);
+  hist_dist_on->SetTitle("DIST distribution");
+  hist_dist_off->SetTitle("DIST distribution");
+
+  hist_dist_on->DrawCopy("E1P");
+
+  // move stat box to make them all visible
+  gPad->Update();
+  TPaveStats* pavs_on_dist = (TPaveStats*) hist_dist_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();
+
+  hist_dist_off->DrawCopy("HISTSAME");
+  hist_dist_off->DrawCopy("ESAME");
+  hist_dist_on->DrawCopy("E1PSAME");
+
+  Double_t chidist = ChiSquareNDof((TH1D*)hist_dist_on,(TH1D*)hist_dist_off);
+
+  x_label_pos  = distmax[0]*0.7;
+  y_label_pos  = hist_dist_on->GetBinContent(hist_dist_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();
+
+   // ############################################################################
+  // DrawCopy WIDTH
+  // ############################################################################
+  display->AddTab("WIDTH");
+
+  gPad->cd();
+
+  hist_width_off->Sumw2();
+  hist_width_off->Scale(norm); 
+  hist_width_on->SetLineColor(kBlack);
+  hist_width_on->SetMarkerStyle(21);
+  hist_width_on->SetMarkerSize(0.7);
+  hist_width_on->SetMarkerColor(kBlack);
+  hist_width_off->SetFillColor(46);
+  hist_width_off->SetLineColor(46);
+  hist_width_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_width_off->SetMinimum(0.);
+  hist_width_on->SetTitle("WIDTH distribution");
+  hist_width_off->SetTitle("WIDTH distribution");
+
+  hist_width_on->DrawCopy("E1P");
+
+  // move stat box to make them all visible
+  gPad->Update();
+  TPaveStats* pavs_on_width = (TPaveStats*) hist_width_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();
+
+  hist_width_off->DrawCopy("HISTSAME");
+  hist_width_off->DrawCopy("ESAME");
+  hist_width_on->DrawCopy("E1PSAME");
+
+  Double_t chiwidth = ChiSquareNDof((TH1D*)hist_width_on,(TH1D*)hist_width_off);
+
+  x_label_pos  = widthmax[0]*0.7;
+  y_label_pos  = hist_width_on->GetBinContent(hist_width_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();
+ 
+  // ############################################################################
+  // DrawCopy LENGTH
+  // ############################################################################
+  display->AddTab("LENGTH");
+ 
+  gPad->cd();
+
+  hist_length_on->Sumw2();
+  hist_length_off->Sumw2();
+  hist_length_off->Scale(norm); 
+  hist_length_on->SetLineColor(kBlack);
+  hist_length_on->SetMarkerStyle(21);
+  hist_length_on->SetMarkerSize(0.7);
+  hist_length_on->SetMarkerColor(kBlack);
+  hist_length_off->SetFillColor(46);
+  hist_length_off->SetLineColor(46);
+  hist_length_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_length_off->SetMinimum(0.);
+  hist_length_on->SetTitle("LENGTH distribution");
+  hist_length_off->SetTitle("LENGTH distribution");
+
+  hist_length_on->DrawCopy("E1P");
+
+  // move stat box to make them all visible
+  gPad->Update();
+  TPaveStats* pavs_on_length = (TPaveStats*) hist_length_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();
+
+  hist_length_off->DrawCopy("HISTSAME");
+  hist_length_off->DrawCopy("ESAME");
+  hist_length_on->DrawCopy("E1PSAME");
+
+  Double_t chilength = ChiSquareNDof((TH1D*)hist_length_on,(TH1D*)hist_length_off);
+
+  x_label_pos  = lengthmax[0]*0.7;
+  y_label_pos  = hist_length_on->GetBinContent(hist_length_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();
+
+  // ############################################################################
+  // DrawCopy normalized ALPHA plot
+  // ############################################################################
+  display->AddTab("ALPHA");
+  
+  gPad->cd();
+
+  hist_on_abs->Sumw2();
+  hist_off_abs->SetStats(0);
+  hist_off_abs->Sumw2();
+  hist_off_abs->Scale(norm); 
+  hist_on_abs->SetStats(0); //-> Do NOT show the legend with statistics
+  hist_on_abs->SetLineColor(kBlack);
+  hist_on_abs->SetMarkerStyle(21);
+  //hist_on_abs->SetMarkerSize();
+  hist_on_abs->SetMarkerColor(kBlack);
+  hist_on_abs->SetMarkerSize(0.7);
+  hist_off_abs->SetFillColor(46);
+  hist_off_abs->SetLineColor(46);
+  hist_off_abs->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_off_abs->SetMinimum(0.);
+  hist_on_abs->SetTitle("Alpha plot");
+  hist_off_abs->SetTitle("Alpha plot");
+
+  
+  hist_on_abs->DrawCopy("E1P");
+  hist_off_abs->DrawCopy("HISTSAME");
+  hist_off_abs->DrawCopy("ESAME");
+  hist_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 hist_on_entries  = (int) hist_on_abs->GetEntries();
+  int hist_off_entries = (int) hist_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",  hist_on_entries);
+  sprintf(text_entries_off, " Entries OFF  = %d",  hist_off_entries);
+  sprintf(text_integral_off," Integral OFF = %d",  int_off);
+  
+  x_label_pos  = alphamax[0]*0.7;
+  y_label_pos  = (hist_on_abs->GetBinContent(hist_on_abs->GetMaximumBin()))/1.6; //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*)hist_on_abs,(TH1D*)hist_off_abs);
+
+  y_label_pos  = (hist_on_abs->GetBinContent(hist_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();
+
+  // ############################################################################
+  // DrawCopy normalized alpha histos for alpha form -90 to 90 deg.
+  // ############################################################################
+  display->AddTab("ALPHA +-90");
+
+  gPad->cd();
+
+  hist_on->Sumw2();
+  hist_off->SetStats(0);
+  hist_off->Sumw2();
+  hist_off->Scale(norm); 
+  hist_off->SetFillColor(46);
+  hist_off->SetLineColor(46);
+  hist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines
+  hist_off->SetMinimum(0.); 
+  hist_on->SetStats(0); //-> Do NOT show the legend with statistics
+  hist_on->SetLineColor(kBlack);
+  hist_on->SetMarkerStyle(21);
+  hist_on->SetMarkerSize(0.7);
+  hist_on->SetMarkerColor(kBlack);
+  hist_on->SetTitle("Alpha plot form -90 to 90 deg");
+  hist_off->SetTitle("Alpha plot form -90 to 90 deg");
+
+  hist_on->DrawCopy("E1P");
+  hist_off->DrawCopy("HISTSAME");
+  hist_off->DrawCopy("ESAME");
+  hist_on->DrawCopy("E1PSAME");
+
+  Double_t chialpha90 = ChiSquareNDof((TH1D*)hist_on,(TH1D*)hist_off);
+
+  x_label_pos  = alphamax[0]*0.5;
+  y_label_pos  = hist_on->GetBinContent(hist_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();
+
+
+  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;
+
+
+  display->AddTab("SRCPOS ON");
+  TH2F *hist_srcpos_on = (TH2F*)hSrcPos_on.GetHist();
+  hist_srcpos_on->DrawCopy("BOX");
+
+  display->AddTab("SRCPOS OFF");
+  TH2F *hist_srcpos_off = (TH2F*)hSrcPos_off.GetHist();
+  hist_srcpos_off->DrawCopy("BOX");
+
+//   display->AddTab("EVTNUM ON");
+//   TH1F *hist_evtnum_on = (TH1F*)hDAQEvtNumber_on.GetHist();
+//   hist_evtnum_on->DrawCopy();
+
+//   display->AddTab("EVTNUM OFF");
+//   TH1F *hist_evtnum_off = (TH1F*)hDAQEvtNumber_off.GetHist();
+//   hist_evtnum_off->DrawCopy();
+
+    }  
+
+  cout << "Done!!" <<endl;
+  
+}
+
+
+