/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Raquel de los Reyes, 03/2004 <mailto:reyes@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== *//////////////////////////////////////////////////////////////////////////////
//
// This macro made the check of the DAQ files (.raw files).
//
// WARNING!: This macro only works if the directory is alphabetically ordered
//
///////////////////////////////////////////////////////////////////////////
#include "MAGIC.h"

void DataAnalysis(const MReadMarsFile *pedfiles,const MReadMarsFile *calfiles,const MReadMarsFile *datafiles,TString directory="./")
{

  TVirtualPad *c;

  Bool_t batchmode = kTRUE;

  if(!batchmode)
    {
      MStatusDisplay *d = new MStatusDisplay;
      d->SetLogStream(&gLog, kTRUE);            // Disables output to stdout  
    }

  // 
  // Set up the source run-range
  //
  TRegexp run("_[0-9][0-9][0-9][0-9][0-9]_");
  TRegexp type("_[A-Z]_");
  TString ped = "", cal = "", data = "", source = "";
  if(pedfiles->GetEntries()!=0)
    source = pedfiles->GetFileName();
  else if(calfiles->GetEntries()!=0)
    source = calfiles->GetFileName();
  else
    source = datafiles->GetFileName();

  ped = pedfiles->GetFileName();
  if(calfiles->GetEntries()!=0)
    cal = calfiles->GetFileName();
  if(datafiles->GetEntries()!=0)
    data = datafiles->GetFileName();
  
  ped = ped(ped.Index(run)+1,5);
  cal = cal(cal.Index(run)+1,5);
  data = data(data.Index(run)+1,5);

  source = source(source.Index(type)+3,source.Length());
  source.Remove(source.Last('_'),source.Length());

//    cout << ped << endl;
//    cout << cal << endl;
//    cout << data << endl;
//    cout << source << endl;

  // **********************************************************************
  // ***************************** PEDESTALS ******************************
  // **********************************************************************
  if(pedfiles->GetEntries()==0)
    {
      cout << "Warning, no entries found in pedestal run(s)!!!"<< endl;
      break;
    }
  MParList plist;
  MTaskList tlist;
  plist.AddToList(&tlist);
  
//    pedfiles->Print("all");
  //
  // Now setup the tasks and tasklist for the pedestals:
  // ---------------------------------------------------
  //
  
  //tasks
  pedfiles->DisableAutoScheme();
  tlist.AddToList(pedfiles);
  
  // pedestal tasks and containers
  MGeomCamMagic  geomcam;
  plist.AddToList(&geomcam);
  MPedestalCam   pedcam;
  plist.AddToList(&pedcam);
  
  MGeomApply geomapl;
  tlist.AddToList(&geomapl);
  MPedCalcPedRun pedcalc;
  tlist.AddToList(&pedcalc);
  
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
  if(!batchmode)
    evtloop.SetDisplay(d);
  
  //
  // Execute first analysis
  //
  cout << "*************************" << endl;
  cout << "** COMPUTING PEDESTALS **" << endl;
  cout << "*************************" << endl;
  if (!evtloop.Eventloop())
    return;
  
  tlist.PrintStatistics();
  
  //
  // Display the pedestal checking plots
  //
  if ((d = evtloop.GetDisplay()))
    TCanvas &c1 = d.AddTab("PEDESTALS");
  else
    TCanvas *c1 = new TCanvas();
  MHCamera meanped(geomcam,"MPedestalCam;ped","Pedestals");
  meanped.SetCamContent(pedcam, 0);
  meanped.SetCamError(pedcam, 1);
  MHCamera rmsped(geomcam,"MPedestalCam;var","Sigma Pedestal");
  rmsped.SetCamContent(pedcam, 2);
  c1->Divide(1,2);
  gPad->SetBorderMode(0);
  c1->cd(1);
  meanped.DrawClone("nonewEPhist");
  c1->cd(2);
  gPad->SetBorderMode(0);
  c = gPad;
  c->Divide(2,1);
  c->cd(1);
  meanped.DrawClone("nonewpixelindex");
  c->cd(2);
  rmsped.DrawClone("nonewpixelindex");

  //
  // Save canvas/display into a postcript file
  //

  if ((d = evtloop.GetDisplay()))
    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
  else
    c1->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps(");

  if(calfiles->GetEntries()==0)
    c1->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps]");

  // **********************************************************************
  // ***************************** CALIBRATION ****************************
  // **********************************************************************
  //
  // Create a empty Parameter List and an empty Task List
  //
  if(calfiles->GetEntries()==0)
    {
      cout << "Warning, no entries found in calibration run(s)!!!"<< endl;
      break;
    }
  MParList  plist2;
  MTaskList tlist2;
  plist2.AddToList(&tlist2);
  
//    calfiles->Print("all");
  //
  // Now setup the new tasks and tasklist for the calibration
  // ---------------------------------------------------
  //
  
  calfiles->DisableAutoScheme();
  tlist2.AddToList(calfiles);  
  
  MExtractedSignalCam        sigcam;
  MArrivalTimeCam            timecam;
  MBadPixelsCam              badcam;
  MCalibrationChargeCam      calcam;
  MCalibrationChargePINDiode pindiode;
  MCalibrationChargeBlindPix blindpix;
  
  MHCalibrationRelTimeCam     histtime;
  MHCalibrationChargeCam      histcharge;
  MHCalibrationChargePINDiode histpin;
  MHCalibrationChargeBlindPix histblind;
  //
  // As long, as we don't have digital modules,
  // we have to set the color of the pulser LED by hand
  //
  blindpix.SetColor(kCT1);
  pindiode.SetColor(kCT1);
  //
  // Get the previously created MPedestalCam and MGeomCamMagic 
  // into the new Parameter List
  //
  plist2.AddToList(&pedcam);
  plist2.AddToList(&geomcam);

  plist2.AddToList(&sigcam);
  plist2.AddToList(&timecam);
  plist2.AddToList(&badcam);
  plist2.AddToList(&calcam);
  plist2.AddToList(&histtime);
  plist2.AddToList(&histcharge);
  plist2.AddToList(&histpin);
  plist2.AddToList(&histblind);
  
  //
  // We saw that the signal jumps between slices,
  // thus take the sliding window
  //
  MExtractSignal2        sigcalc2;
  MExtractPINDiode       pincalc;
  MExtractBlindPixel     blindcalc;
  MArrivalTimeCalc2      timecalc;
  MCalibrationChargeCalc calcalc;
  
  MFillH filltime( "MHCalibrationRelTimeCam"    , "MArrivalTimeCam");
  MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); 
  MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
  MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");
  
  // Skip the MFillH draw function
  filltime.SetBit(MFillH::kDoNotDisplay);
  fillpin.SetBit(MFillH::kDoNotDisplay);
  fillblind.SetBit(MFillH::kDoNotDisplay);
  fillcam.SetBit(MFillH::kDoNotDisplay);
  
  //
  // Skip the HiGain vs. LoGain calibration
  //
  calcalc.SkipHiLoGainCalibration();
  
  //
  // Apply a filter against cosmics
  // (was directly in MCalibrationCalc in earlier versions)
  //
  MFCosmics            cosmics;
  MContinue            cont(&cosmics);
  //
  // Get the previously created MGeomApply into the new Task List
  //  
  tlist2.AddToList(&geomapl);

  tlist2.AddToList(&sigcalc2);
  tlist2.AddToList(&pincalc);
  tlist2.AddToList(&blindcalc);
  //
  // In case, you want to skip the cosmics rejection,
  // uncomment the next line
  //
  tlist2.AddToList(&cont);
  //
  // In case, you want to skip the somewhat lengthy calculation
  // of the arrival times using a spline, uncomment the next two lines
  //
  tlist2.AddToList(&timecalc);
  tlist2.AddToList(&filltime);

  tlist2.AddToList(&fillpin);
  tlist2.AddToList(&fillblind);
  tlist2.AddToList(&fillcam);
  tlist2.AddToList(&calcalc);
  
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop2;
  evtloop2.SetParList(&plist2);
  if(!batchmode)
    evtloop2.SetDisplay(d);
  
  //
  // Execute second analysis
  //
  cout << "***********************************" << endl;
  cout << "** COMPUTING CALIBRATION FACTORS **" << endl;
  cout << "***********************************" << endl;
  if (!evtloop2.Eventloop())
    return;
  
  tlist2.PrintStatistics();
  
  //
  // print the most important results of all pixels to a file
  //
  //    MLog gauglog;
  //    gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
  //    calcam.SetLogStream(&gauglog);
  //    calcam.Print();
  //    calcam.SetLogStream(&gLog);

  TH1 *hist;

  //
  // Display the calibration checking plots
  //      
  // ---------------------- Conversion factor ---------------------------
  if ((d = evtloop2.GetDisplay()))
    TCanvas &c2 = d.AddTab("CONV FACTOR");
  else
    TCanvas *c2 = new TCanvas();
  c2->Divide(1,2);
  c2->cd(1);
  gPad->SetBorderMode(0);
  c = gPad;
  c->Divide(2,1);
  gPad->SetBorderMode(0);
  c->cd(1);
  gPad->SetBorderMode(0);
  MHCamera ConvFfactor(geomcam,"Cal;FFactorConv","Conversion Factor to photons (F-Factor method)");
  ConvFfactor.SetCamContent(calcam, 11);
  ConvFfactor.DrawClone("nonewpixelindex");
  c->cd(2);
  MHCamera ConvBlindPix(geomcam,"Cal;BlindPixConv","Conversion Factor to photons (Blind Pixel method)");
  ConvBlindPix.SetCamContent(calcam, 17);
  ConvBlindPix.DrawClone("nonewpixelindex");
  c2->cd(2);
  gPad->SetBorderMode(0);
  MHCamera ConvPINDiode(geomcam,"Cal;PINDiodeConv","Conversion Factor to photons (PINDiode method)");
  ConvPINDiode.SetCamContent(calcam, 23);
  ConvPINDiode.DrawClone("nonewpixelindex");
  // ----------------- Pixels with defects ------------------------------
  if ((d = evtloop2.GetDisplay()))
    TCanvas &c3 = d.AddTab("DEFECT PIXELS");
  else
    TCanvas *c3 = new TCanvas();
  c3->Divide(3,2);
  c3->cd(1);
  MHCamera PixExc(geomcam,"Cal;Excluded","Pixels previously excluded");
  PixExc.SetCamContent(calcam, 27);
  PixExc.DrawClone("nonewpixelindex");
  c3->cd(2);
  MHCamera PixNoFit(geomcam,"Cal;NotFitted","Pixels that could not be fitted");
  PixNoFit.SetCamContent(calcam, 28);
  PixNoFit.DrawClone("nonewpixelindex");
  c3->cd(3);
  MHCamera PixNoValid(geomcam,"Cal;NotValidFit","Pixels with not valid fit results");
  PixNoValid.SetCamContent(calcam, 29);
  PixNoValid.DrawClone("nonewpixelindex");
  c3->cd(4);
  MHCamera PixFFactor(geomcam,"Cal;FFactorValid","Pixels with valid F-Factor calibration");
  PixFFactor.SetCamContent(calcam, 35);
  PixFFactor.DrawClone("nonewpixelindex");
  c3->cd(5);
  MHCamera PixBlind(geomcam,"Cal;BlindPixelValid","Pixels with valid BlindPixel calibration");
  PixBlind.SetCamContent(calcam, 36);
  PixBlind.DrawClone("nonewpixelindex");
  c3->cd(6);
  MHCamera PixPINDiode(geomcam,"Cal;PINdiodeFFactorValid","Pixels with valid PINDiode calibration");
  PixPINDiode.SetCamContent(calcam, 37);
  PixPINDiode.DrawClone("nonewpixelindex");
  // ------------------------- Fitted mean charges ----------------------------
  if ((d = evtloop2.GetDisplay()))
    TCanvas &c4 = d.AddTab("FIT CHARGES");
  else
    TCanvas *c4 = new TCanvas();
  c4->Divide(1,2);
  MHCamera fitmeancharge(geomcam,"Cal;Charge","Fitted Mean Charges");
  fitmeancharge.SetCamContent(calcam, 0);
  fitmeancharge.SetCamError(calcam, 1);
  c4->cd(1);
  gPad->SetBorderMode(0);
  fitmeancharge.DrawClone("nonewEPhist");
  c4->cd(2);
  c = gPad;
  gPad->SetBorderMode(0);
  c->Divide(2,1);
  c->cd(1);
  gPad->SetBorderMode(0);
  fitmeancharge.DrawClone("nonewpixelindex");
  c->cd(2);
  gPad->SetBorderMode(0);
  hist = (TH1*)fitmeancharge->Projection("Charge");
  hist->SetXTitle("Charge");
  hist->DrawCopy("hist");
  // ----------------------- Reduced Sigma per Charge -------------------------
  if ((d = evtloop2.GetDisplay()))
    TCanvas &c5 = d.AddTab("RED SIGMA");
  else
    TCanvas *c5 = new TCanvas();
  c5->Divide(1,2);
  MHCamera redsigmacharge(geomcam,"Cal;RSigmaCharge","Reduced Sigma per Charge");
  redsigmacharge.SetCamContent(calcam, 7);
  redsigmacharge.SetCamError(calcam, 8);
  c5->cd(1);
  gPad->SetBorderMode(0);
  redsigmacharge.DrawClone("nonewEPhist");
  c5->cd(2);
  c = gPad;
  gPad->SetBorderMode(0);
  c->Divide(2,1);
  c->cd(1);
  gPad->SetBorderMode(0);
  redsigmacharge.DrawClone("nonewpixelindex");
  c->cd(2);
  gPad->SetBorderMode(0);
  hist = redsigmacharge->Projection("RSigmaCharge");
  hist->SetXTitle("Reduced Sigma per Charge");
  hist->Draw("hist");
  // ----------------------- Absolute arrival times -------------------------
  if ((d = evtloop2.GetDisplay()))
    TCanvas &c6 = d.AddTab("ABS TIMES");
  else
    TCanvas *c6 = new TCanvas();
  c6->Divide(1,2);
  MHCamera abstime(geomcam,"Cal;AbsTimeMean","Absolute Arrival Times");
  abstime.SetCamContent(calcam, 42);
  abstime.SetCamError(calcam, 43);
  c6->cd(1);
  gPad->SetBorderMode(0);
  abstime.DrawClone("nonewEPhist");
  c6->cd(2);
  c = gPad;
  gPad->SetBorderMode(0);
  c->Divide(2,1);
  c->cd(1);
  gPad->SetBorderMode(0);
  abstime.DrawClone("nonewpixelindex");
  c->cd(2);
  gPad->SetBorderMode(0);
  hist = (TH1*)abstime->Projection("AbsTimeMean");
  hist->SetXTitle("Absolute arrival time (time slice)");
  hist->DrawCopy("hist");

  //
  // Save canvas/display into a postcript file
  //

  if ((d = evtloop2.GetDisplay()))
    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
  else
    {
      c2->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
      c3->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
      c4->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
      c5->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
      c6->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
    }

  if(calfiles->GetEntries()!=0)
    c6->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps]");

  // ********************************************************************
  // ***************************** DATA *********************************
  // ********************************************************************

  if(datafiles->GetEntries()==0)
    {
      cout << "Warning, no entries found in data run(s)!!!"<< endl;
      break;
    }
  //
  // Create an empty Parameter List and an empty Task List
  //
  MParList  plist3;
  MTaskList tlist3;
  plist3.AddToList(&tlist3);
  
//    datafiles->Print("all");
  //
  // Now setup the new tasks and tasklist for the data
  // ---------------------------------------------------
  //
  
  datafiles->DisableAutoScheme();
  tlist3.AddToList(datafiles);  
  
  // Task containers
  MExtractedSignalCam sigcam;
  //
  // Get the previously created MPedestalCam and MGeomCamMagic
  // into the new Parameter List
  //
  plist3.AddToList(&geomcam);
  plist3.AddToList(&pedcam);

  plist3.AddToList(&sigcam);

  // Display containers
  MHTimeDiffTime difftime("DiffTime","Differential time between events");
  plist3.AddToList(&difftime);
  MBinning binsdifftime("BinningTimeDiff");
  binsdifftime.SetEdges(200, 0., 0.2);
  plist3.AddToList(&binsdifftime);
  MBinning binstime("BinningTime");
  binstime.SetEdges(10000, 0., 1e10);
  plist3.AddToList(&binstime);

  // Data tasks
  MExtractSignal2 signal2;

  MFillH fillDiffTime(&difftime,"MRawEvtData");
  fillDiffTime.SetBit(MFillH::kDoNotDisplay);
  tlist3.AddToList(&fillDiffTime);

  //
  // Get the previously created MGeomApply into the new Task List
  //
  tlist3.AddToList(&geomapl);
  tlist3.AddToList(&signal2);

  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop3;
  evtloop3.SetParList(&plist3);
  if(!batchmode)
    evtloop3.SetDisplay(d);
  
  //
  // Execute third analysis
  //
  cout << "********************" << endl;
  cout << "** COMPUTING DATA **" << endl;
  cout << "********************" << endl;
  if (!evtloop3.Eventloop())
    return;
  
  tlist3.PrintStatistics();

  TH1 *hist1;
  TH2 *hist2;

  //
  // Display the data checking canvas
  //
  // -------------------- Diffential time of events ---------------------
  if ((d = evtloop3.GetDisplay()))
    TCanvas &c9 = d.AddTab("TIME");
  else
    TCanvas *c9 = new TCanvas();
  gStyle->SetPadGridX(kTRUE);
  gStyle->SetPadGridY(kTRUE);
  gStyle->SetOptFit();
  c9->Divide(2,2);
  hist2 = difftime.GetHist();
  c9->cd(1);
  gPad->SetLogy();
  hist1 = hist2->ProjectionX("ProjX_sumtime", -1, 9999, "E");  
  hist1->SetTitle("Distribution of \\Delta t [s]");
  hist1->SetXTitle("\\Delta t [s]");
  hist1->SetYTitle("Counts");
  hist1->GetXaxis()->SetRange(0,hist1->GetXaxis()->GetNbins());
  TF1 *f1 = new TF1("f1","[0]*exp(x*[1])",hist1->GetMaximumBin(),hist1->GetXaxis()->GetXmax());
  hist1->Fit(f1,"");
  hist1->DrawCopy();
  c9->cd(2);
  gPad->SetLogy();
  hist1 = hist2->ProjectionX("ProjX_sumtime", -1, 9999, "E");  
  hist1->SetTitle("Distribution of \\Delta t [s] (deadtime)");
  hist1->SetXTitle("\\Delta t [s]");
  hist1->SetYTitle("Counts");
//    cout << hist1->GetXaxis()->GetBinCenter(hist1->GetMaximumBin())<< endl;
  hist1->GetXaxis()->SetRange(0,hist1->GetMaximumBin());
  hist1->DrawCopy();
  c9->cd(3);
  hist1 = difftime.GetHist()->ProjectionY("ProjY_sumtime", -1, 9999, "E");  
  hist1->SetTitle("Distribution of time [s]");
  hist1->SetXTitle("Time t [s]");
  hist1->SetYTitle("Counts");
  hist1->DrawCopy();
  c9->cd(4);

  //
  // Save canvas/display into a postcript file
  //

  if ((d = evtloop3.GetDisplay()))
    d->SaveAsPS(directory+source+"_"+ped+"-"+cal+"-"+data+".ps");
  else
    {
      c9->Print(directory+source+"_"+ped+"-"+cal+"-"+data+".ps)");
    }
  
}

void DAQDataCheck(const TString directory="/remote/data2/data/Magic-DATA/Period014/rootdata/2004_02_10/")
{

  MDirIter Next;
  Next.AddDirectory(directory,"*_E.root",-1);

//    Next.Print("all");

  TString fName="file.root";

  MReadMarsFile *pedfiles;
  MReadMarsFile *calfiles;
  MReadMarsFile *datafiles;

  TString source="";

  fName = Next();  // absolut path

  while(!fName.IsNull())
    {

      TString file = fName(fName.Last('/')+1,fName.Length()); // root file name
      file = file(0,file.Last('_'));

//        TString run = file(file.First('_')+1,5);
//        TString type = file(file.Index('_',file.First('_')+1)+1,1);
      source = file(file.Last('_')+1,file.Length());

      // Pedestal runs
      if(fName.Contains("_P_"))
	pedfiles = new MReadMarsFile("Events");
      else
	pedfiles->Rewind();

      while(fName.Contains("_P_")&&fName.Contains(source))
	{
	  pedfiles->AddFile(fName);
	  fName = Next();  
	}

      // Calibration runs
      if(fName.Contains("_C_")||(!fName.Contains(source)))
	calfiles = new MReadMarsFile("Events");
      else
	calfiles->Rewind();

      while(fName.Contains("_C_")&&fName.Contains(source))
	{
	  calfiles->AddFile(fName);
	  fName = Next();  
	}

      // Data runs
      if(fName.Contains("_D_")||(!fName.Contains(source)))
	datafiles = new MReadMarsFile("Events");
//        else
//  	break;
      while(fName.Contains("_D_")&&fName.Contains(source))
	{
	  datafiles->AddFile(fName);
	  fName = Next();  
	}

//        pedfiles->Print();
//        calfiles->Print("all");
//        datafiles->Print("all");

      DataAnalysis(pedfiles,calfiles,datafiles,directory);

      cout << "----------------------------------------------" << endl;

    }

}