/* ======================================================================== *\
!
! *
! * 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): Josep  Flix,  01/2004 <mailto:jflix@ifae.es>
!              Javier Rico,  01/2004 <mailto:jrico@ifae.es>
!
!   (based on bootcampstandardanalysis.C by Javier López)
!
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

#include "MParList.h"
#include "MTaskList.h"
#include "MPedestalCam.h"
#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MPedCalcPedRun.h"
#include "MGeomCamMagic.h"
#include "MEvtLoop.h"
#include "MCalibrationCam.h"
#include "MExtractedSignalCam.h"
#include "MExtractSignal.h" 
#include "MCalibrationCalc.h"
#include "MCerPhotEvt.h"
#include "MCalibrate.h"
#include "MPedPhotCalc.h"
#include "MPedPhotCam.h"
#include "MHCamera.h"

#include "TTimer.h"
#include "TString.h"
#include "TCanvas.h"

#include <iostream>


const TString pedfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
const TString calfile="/disc02/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root"; 

void pedphotcalc(TString pedname=pedfile,TString calname=calfile)
{
  
  /*
    This macro is an example of the use of MPedPhotCalc. It computes
    the pedestal mean and rms from pedestal files undergoing the
    signal extraction + calibration chain, in units of photons. It
    produces plots of the relevant computed quantities.
  */
  
  /************************************/
  /* FIRST LOOP: PEDESTAL COMPUTATION */
  /************************************/
  
  MParList plist;
  MTaskList tlist;
  plist.AddToList(&tlist);
  
  // containers
  MPedestalCam   pedcam;
  plist.AddToList(&pedcam);
    
  //tasks
  MReadMarsFile read("Events", pedname);
  MGeomApply     geomapl;
  MPedCalcPedRun pedcalc;
  MGeomCamMagic  geomcam;

  read.DisableAutoScheme();
    
  tlist.AddToList(&read);
  tlist.AddToList(&geomapl);
  tlist.AddToList(&pedcalc);
  
  // Create and execute the event looper
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
  
  cout << "*************************" << endl;
  cout << "** COMPUTING PEDESTALS **" << endl;
  cout << "*************************" << endl;
  if (!evtloop.Eventloop())
    return;
  
  tlist.PrintStatistics();

  /****************************/
  /* SECOND LOOP: CALIBRATION */
  /****************************/
  
  MParList  plist2;
  MTaskList tlist2;
  plist2.AddToList(&tlist2);
  
  // containers  
  MCalibrationCam     calcam;
  MExtractedSignalCam sigcam;
  
  plist2.AddToList(&geomcam);
  plist2.AddToList(&pedcam);
  plist2.AddToList(&calcam);
  plist2.AddToList(&sigcam);
  
  //tasks
  MReadMarsFile read2("Events", calname);
  read2.DisableAutoScheme();
  
  MExtractSignal       sigcalc;
  MCalibrationCalc     calcalc;
  
  //
  // As long as we don't have digital modules, 
  // we have to set the color of the pulser LED by hand
  //
  calcalc.SetPulserColor(MCalibrationCalc::kECT1);
  
  tlist2.AddToList(&read2);
  tlist2.AddToList(&geomapl);
  tlist2.AddToList(&sigcalc);
  tlist2.AddToList(&calcalc);
  
  // Execute second loop (calibration)
  MEvtLoop evtloop2;
  evtloop2.SetParList(&plist2);
  
  cout << "***********************************" << endl;
  cout << "** COMPUTING CALIBRATION FACTORS **" << endl;
  cout << "***********************************" << endl;
  if (!evtloop2.Eventloop())
    return;
  
  tlist2.PrintStatistics();

  
  /************************************************************************/
  /* THIRD LOOP: PEDESTAL COMPUTATION USING EXTRACTED SIGNAL (IN PHOTONS) */
  /************************************************************************/
  
  // Create an empty Parameter List and an empty Task List
  MParList  plist3;  
  MTaskList tlist3;
  plist3.AddToList(&tlist3);

  // containers
  MCerPhotEvt    photcam;
  MPedPhotCam    pedphotcam;
  MExtractedSignalCam extedsig;
  
  plist3.AddToList(&geomcam);
  plist3.AddToList(&pedcam);
  plist3.AddToList(&calcam);
  plist3.AddToList(&photcam);
  plist3.AddToList(&extedsig);
  plist3.AddToList(&pedphotcam);
  
  //tasks
  MReadMarsFile read3("Events", pedname);
  MCalibrate      photcalc;
  MPedPhotCalc    pedphotcalc;  

  read3.DisableAutoScheme();
  
  tlist3.AddToList(&read3);
  tlist3.AddToList(&geomapl);
  tlist3.AddToList(&sigcalc);
  tlist3.AddToList(&photcalc);
  tlist3.AddToList(&pedphotcalc);
  
  // Create and execute eventloop
  MEvtLoop evtloop3;
  evtloop3.SetParList(&plist3);
    
  cout << "*************************************************************" << endl;
  cout << "** COMPUTING PEDESTALS USING EXTRACTED SIGNAL (IN PHOTONS) **" << endl;
  cout << "*************************************************************" << endl;
  
  if (!evtloop3.Eventloop())  
    return;  
  tlist3.PrintStatistics();


  /**********************/
  /* PRODUCE NICE PLOTS */
  /**********************/

  const UShort_t npix = 577;

  // declare histograms
  // pedestals
  TH1F* pedhist    = new TH1F("pedhist","Pedestal",npix,0,npix);
  TH1F* pedrmshist = new TH1F("pedrmshist","Pedestal RMS",npix,0,npix);
  TH1F* peddist    = new TH1F("peddist","Pedestal",100,0,20);
  TH1F* pedrmsdist = new TH1F("pedrmsdist","Pedestal RMS",100,0,15);
  
  // calibration factors
  TH1F* calhist    = new TH1F("calhist","Calibration factors",npix,0,npix);
  TH1F* caldist    = new TH1F("caldist","Calibration factors",100,0,1);

  // pedestal signals
  TH1F* pedphothist    = new TH1F("pedphothist","Pedestal Signal",npix,0,npix);
  TH1F* pedphotrmshist = new TH1F("pedphotrmshist","Pedestal Signal RMS",npix,0,npix);
  TH1F* pedphotdist    = new TH1F("pedphotdist","Pedestal Signal",100,-0.4,0.4);
  TH1F* pedphotrmsdist = new TH1F("pedphotrmsdist","Pedestal Signal RMS",100,0,25);

  // fill histograms
  for(int i=0;i<npix;i++)
    {
      const Float_t ped        = pedcam[i].GetPedestal();
      const Float_t pedrms     = pedcam[i].GetPedestalRms();
      const Float_t cal        = calcam[i].GetMeanConversionBlindPixelMethod();
      const Float_t calerr     = calcam[i].GetErrorConversionBlindPixelMethod();
      const Float_t pedphot    = pedphotcam[i].GetMean();
      const Float_t pedphotrms = pedphotcam[i].GetRms();

      pedhist->Fill(i,ped);
      peddist->Fill(ped);
      pedrmshist->Fill(i,pedrms);
      pedrmsdist->Fill(pedrms);

      calhist->Fill(i,cal);
      caldist->Fill(cal);

      pedphothist->Fill(i,pedphot);
      pedphotdist->Fill(pedphot);
      pedphotrmshist->Fill(i,pedphotrms);
      pedphotrmsdist->Fill(pedphotrms);
    }

  // Draw
  gROOT->Reset();
  gStyle->SetCanvasColor(0);
  TCanvas* canvas = new TCanvas("canvas","pedphotcalc.C", 0, 100, 650, 800);
  canvas->SetBorderMode(0);    // Delete the Canvas' border line     
  canvas->cd();
  canvas->SetBorderMode(0);
     
  canvas->Divide(2,4);

  // draw pedestal histo
  canvas->cd(1);
  gPad->cd();
  gPad->SetBorderMode(0);

  pedhist->SetStats(kFALSE);
  pedhist->GetXaxis()->SetTitle("Pixel SW Id");
  pedhist->GetYaxis()->SetTitle("Pedestal (ADC counts)");
  pedrmshist->SetStats(kFALSE);
  pedrmshist->SetLineColor(2);
  pedhist->Draw();
  pedrmshist->Draw("same");

  TLegend* leg1 = new TLegend(.14,.68,.39,.88);
  leg1->SetHeader("");
  leg1->AddEntry(pedhist,"Pedestal","L");
  leg1->AddEntry(pedrmshist,"Pedestal RMS","L");
  leg1->SetFillColor(0);
  leg1->SetLineColor(0);
  leg1->SetBorderSize(0);
  leg1->Draw();
     
  // draw pedestal distribution
  canvas->cd(2);
  gPad->cd();
  gPad->SetBorderMode(0);
  peddist->GetXaxis()->SetTitle("Pedestal (ADC counts)");
  pedrmsdist->SetLineColor(2);
  peddist->Draw();
  pedrmsdist->Draw("same");


  TLegend* leg2 = new TLegend(.14,.68,.39,.88);
  leg2->SetHeader("");
  leg2->AddEntry(pedhist,"Pedestal","L");
  leg2->AddEntry(pedrmshist,"Pedestal RMS","L");
  leg2->SetFillColor(0);
  leg2->SetLineColor(0);
  leg2->SetBorderSize(0);
  leg2->Draw();

  // draw calibration histo
  canvas->cd(3);
  gPad->cd();
  gPad->SetBorderMode(0);
  calhist->GetXaxis()->SetTitle("Pixel SW Id");
  calhist->SetMaximum(1);
  calhist->SetMinimum(0);
  calhist->GetYaxis()->SetTitle("Calibration factor (photon/ADC count)");
  calhist->SetStats(kFALSE);
  calhist->Draw();

  // draw calibration distribution
  canvas->cd(4);
  gPad->cd();
  gPad->SetBorderMode(0);
  caldist->GetXaxis()->SetTitle("Calibration factor (photon/ADC count)");
  caldist->Draw();

  // draw pedestal signal histo
  canvas->cd(5);
  gPad->cd();
  gPad->SetBorderMode(0);
  pedphothist->GetXaxis()->SetTitle("Pixel SW Id");
  pedphothist->GetYaxis()->SetTitle("Pedestal signal (photons)");
  pedphothist->SetStats(kFALSE);
  pedphothist->Draw();

  // draw pedestal signal distribution
  canvas->cd(6);
  gPad->cd();
  gPad->SetBorderMode(0);
  pedphotdist->GetXaxis()->SetTitle("Pedestal signal (photons)");
  pedphotdist->Draw();

  // draw pedestal signal rms histo
  canvas->cd(7);
  gPad->cd();
  gPad->SetBorderMode(0);
  pedphotrmshist->GetXaxis()->SetTitle("Pixel SW Id");
  pedphotrmshist->GetYaxis()->SetTitle("Pedestal signal rms (photons)");
  pedphotrmshist->SetStats(kFALSE);
  pedphotrmshist->Draw();

  // draw pedestal signal rms distribution
  canvas->cd(8);
  gPad->cd();
  gPad->SetBorderMode(0);
  pedphotrmsdist->GetXaxis()->SetTitle("Pedestal signal rms (photons)");
  pedphotrmsdist->Draw();

  return; 
}