/* ======================================================================== *\
!
! *
! * 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): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */
#include "MAGIC.h"


void calibration()
{

//    const char *inpath = "/remote/home/pc2/operator/NewCalBoxTest/";
//    const char *inpath = "/remote/home/pc2/operator/Crab20040214/";
    const char *inpath = "./";
    //const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/";
    //const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_01_26/";
    //const TString inpath = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/";

    MRunIter pruns;
    MRunIter cruns;

    pruns.AddRun(14400,inpath);
    cruns.AddRuns(14401,14409,inpath);
//    pruns.AddRuns(15487,15494,inpath);
//    cruns.AddRuns(15495,15508,inpath);
//    pruns.AddRun(19816,inpath);  // CalboxTest
//    cruns.AddRun(19820,inpath); // 3blue
//    cruns.AddRun(19820,inpath);   // 1blue

    gStyle->SetOptStat(1111);
    gStyle->SetOptFit();

    MStatusDisplay *display = new MStatusDisplay;
    display->SetUpdateTime(3000);
    display->Resize(850,700);

    MJPedestal pedloop;
    pedloop.SetInput(&pruns);
    pedloop.SetDisplay(display);

    if (!pedloop.Process())
        return;

    // 
    // Now the short version: 
    //
/*
    MJCalibration calloop;
    calloop.SetInput(&cruns);
    calloop.SetDisplay(display);
    if (!calloop.Process(pedloop.GetPedestalCam()))
        return;
#if 0
*/
    // 
    // The longer version:
    //

    //
    // Create a empty Parameter List and an empty Task List 
    //
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);
    plist.AddToList(&pedloop.GetPedestalCam());

    gLog << endl;;
    gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
    gLog << endl;

    MReadMarsFile read("Events");
    read.DisableAutoScheme();
    static_cast<MRead&>(read).AddFiles(cruns);

    MGeomCamMagic              geomcam;
    MExtractedSignalCam        sigcam;
    MArrivalTimeCam            timecam;
    MCalibrationChargeCam      calcam;
    MCalibrationChargePINDiode pindiode;   
    MCalibrationChargeBlindPix blindpix;   
    
    MHCalibrationRelTimeCam     histtime;
    MHCalibrationChargeCam      histcharge;
    MHCalibrationChargePINDiode histpin;
    MHCalibrationChargeBlindPix histblind;
//    histblind.SetSinglePheCut(80);
    //
    // 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 into the new Parameter List 
    //
    plist.AddToList(&geomcam);
    plist.AddToList(&sigcam);
    plist.AddToList(&timecam);
    plist.AddToList(&calcam);
    plist.AddToList(&histtime);
    plist.AddToList(&histcharge);
//    plist.AddToList(&histpin);
    plist.AddToList(&histblind);

    //
    // We saw that the signal jumps between slices, 
    // thus take the sliding window
    //		
    MExtractSignal2        sigcalc2;
    MExtractPINDiode       pincalc;
    MExtractBlindPixel     blindcalc;
//    blindcalc.SetRange(11,16);

    MArrivalTimeCalc2      timecalc;
    MCalibrationChargeCalc calcalc;
    MGeomApply             geomapl;
    
    MFillH filltime( "MHCalibrationRelTimeCam"    , "MArrivalTimeCam");
//   MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
    MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
    MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");

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

    tlist.AddToList(&read);
    tlist.AddToList(&geomapl);
    tlist.AddToList(&sigcalc2);
    tlist.AddToList(&blindcalc);
//    tlist.AddToList(&pincalc);
    //
    // In case, you want to skip the cosmics rejection, 
    // uncomment the next line
    //
    tlist.AddToList(&cont);
    //
    // In case, you want to skip the somewhat lengthy calculation
    // of the arrival times using a spline, uncomment the next two lines
    //
    tlist.AddToList(&timecalc);
    tlist.AddToList(&filltime);
//    tlist.AddToList(&fillpin);
    tlist.AddToList(&fillblind);
    tlist.AddToList(&fillcam);
    //
    tlist.AddToList(&calcalc);
    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(display);
    
    //
    // Execute second analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist.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);
    //
    // just one example how to get the plots of individual pixels
    //
    histblind.DrawClone("all");
    histcharge[5].DrawClone("all");
    histcharge(5).DrawClone("all");
    histtime[5].DrawClone("fourierevents");

    // Create histograms to display
    MHCamera disp1  (geomcam, "Cal;Charge",         "Fitted Mean Charges");
    MHCamera disp2  (geomcam, "Cal;SigmaCharge",    "Sigma of Fitted Charges");
    MHCamera disp3  (geomcam, "Cal;FitProb",        "Probability of Fit");
    MHCamera disp4  (geomcam, "Cal;RSigma",         "Reduced Sigmas");
    MHCamera disp5  (geomcam, "Cal;RSigma/Charge",  "Reduced Sigma per Charge");
    MHCamera disp6  (geomcam, "Cal;FFactorPh",      "Nr. of Photo-electrons (F-Factor Method)");
    MHCamera disp7  (geomcam, "Cal;FFactorConv",    "Conversion Factor to photons (F-Factor Method)");
    MHCamera disp8  (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
    MHCamera disp9  (geomcam, "Cal;BlindPixPh",     "Photon flux inside plexiglass (Blind Pixel Method)");
    MHCamera disp10 (geomcam, "Cal;BlindPixConv",   "Conversion Factor to photons (Blind Pixel Method)");
    MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
    MHCamera disp12 (geomcam, "Cal;PINDiodePh",     "Photon flux outside plexiglass (PIN Diode Method)");
    MHCamera disp13 (geomcam, "Cal;PINDiodeConv",   "Conversion Factor tp photons (PIN Diode Method)");
    MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
    MHCamera disp15 (geomcam, "Cal;Excluded",       "Pixels previously excluded");
    MHCamera disp16 (geomcam, "Cal;NotFitted",      "Pixels that could not be fitted");
    MHCamera disp17 (geomcam, "Cal;NotFitValid",    "Pixels with not valid fit results");
    MHCamera disp18 (geomcam, "Cal;HiGainOscillating",     "Oscillating Pixels HI Gain");
    MHCamera disp19 (geomcam, "Cal;LoGainOscillating",     "Oscillating Pixels LO Gain");
    MHCamera disp20 (geomcam, "Cal;HiGainPickup",          "Number Pickup events Hi Gain");
    MHCamera disp21 (geomcam, "Cal;LoGainPickup",          "Number Pickup events Lo Gain");
    MHCamera disp22 (geomcam, "Cal;Saturation",            "Pixels with saturated Hi Gain");
    MHCamera disp23 (geomcam, "Cal;FFactorValid",          "Pixels with valid F-Factor calibration");
    MHCamera disp24 (geomcam, "Cal;BlindPixelValid",       "Pixels with valid BlindPixel calibration");
    MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid",  "Pixels with valid PINDiode calibration");

    MHCamera disp26 (geomcam, "Cal;Ped",         "Pedestals");
    MHCamera disp27 (geomcam, "Cal;PedRms",      "Pedestal RMS");

    MHCamera disp28 (geomcam, "time;Time",        "Rel. Arrival Times");
    MHCamera disp29 (geomcam, "time;SigmaTime",   "Sigma of Rel. Arrival Times");
    MHCamera disp30 (geomcam, "time;TimeProb",    "Probability of Time Fit");
    MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
    MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels");

    MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
    MHCamera disp34 (geomcam, "Cal;AbsTimeRms",  "RMS of Arrival Times");

    // Fitted charge means and sigmas
    disp1.SetCamContent(calcam,  0);
    disp1.SetCamError(  calcam,  1);
    disp2.SetCamContent(calcam,  2);
    disp2.SetCamError(  calcam,  3);

    // Fit probabilities
    disp3.SetCamContent(calcam,  4);

    // Reduced Sigmas and reduced sigmas per charge
    disp4.SetCamContent(calcam,  5);
    disp4.SetCamError(  calcam,  6);
    disp5.SetCamContent(calcam,  7);
    disp5.SetCamError(  calcam,  8);

    // F-Factor Method
    disp6.SetCamContent(calcam,  9);
    disp6.SetCamError(  calcam, 10);
    disp7.SetCamContent(calcam, 11);
    disp7.SetCamError(  calcam, 12);
    disp8.SetCamContent(calcam, 13);
    disp8.SetCamError(  calcam, 14);

    /// Blind Pixel Method
    disp9.SetCamContent(calcam, 15);
    disp9.SetCamError(  calcam, 16);
    disp10.SetCamContent(calcam,17);
    disp10.SetCamError(  calcam,18);
    disp11.SetCamContent(calcam,19);
    disp11.SetCamError(  calcam,20);

    // PIN Diode Method
    disp12.SetCamContent(calcam,21);
    disp12.SetCamError(  calcam,22);
    disp13.SetCamContent(calcam,23);
    disp13.SetCamError(  calcam,24);
    disp14.SetCamContent(calcam,25);
    disp14.SetCamError(  calcam,26);

    // Pixels with defects
    disp15.SetCamContent(calcam,27);
    disp16.SetCamContent(calcam,28);
    disp17.SetCamContent(calcam,29);
    disp18.SetCamContent(calcam,30);
    disp19.SetCamContent(calcam,31);
    disp20.SetCamContent(calcam,32);
    disp21.SetCamContent(calcam,33);

    // Lo Gain calibration
    disp22.SetCamContent(calcam,34);

    // Valid flags
    disp23.SetCamContent(calcam,35);
    disp24.SetCamContent(calcam,36);
    disp25.SetCamContent(calcam,37);

    // Pedestals
    disp26.SetCamContent(calcam,38);
    disp26.SetCamError(  calcam,39);
    disp27.SetCamContent(calcam,40);
    disp27.SetCamError(  calcam,41);

    // Relative Times
    disp28.SetCamContent(timecam,0);
    disp28.SetCamError(  timecam,1);
    disp29.SetCamContent(timecam,2);
    disp29.SetCamError(  timecam,3);
    disp30.SetCamContent(timecam,4);
    disp31.SetCamContent(timecam,5);
    disp32.SetCamContent(timecam,6);

    // Absolute Times
    disp33.SetCamContent(calcam,42);
    disp33.SetCamError(  calcam,43);
    disp34.SetCamContent(calcam,43);

    disp1.SetYTitle("Charge [FADC units]");
    disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
    disp3.SetYTitle("P_{Charge} [1]");

    disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
    disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");

    disp6.SetYTitle("Nr. Photo-electrons [1]");
    disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
    disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");

    disp9.SetYTitle("Photon flux [ph/mm^2]");
    disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
    disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");

    disp12.SetYTitle("Photon flux [ph/mm^2]");
    disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
    disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");

    disp15.SetYTitle("[1]");
    disp16.SetYTitle("[1]");
    disp17.SetYTitle("[1]");
    disp18.SetYTitle("[1]");
    disp19.SetYTitle("[1]");
    disp20.SetYTitle("[1]");
    disp21.SetYTitle("[1]");
    disp22.SetYTitle("[1]");
    disp23.SetYTitle("[1]");
    disp24.SetYTitle("[1]");
    disp25.SetYTitle("[1]");

    disp26.SetYTitle("Ped [FADC Counts ]");
    disp27.SetYTitle("RMS_{Ped} [FADC Counts ]");

    disp28.SetYTitle("Time Offset [ns]");
    disp29.SetYTitle("Timing resolution [ns]");
    disp30.SetYTitle("P_{Time} [1]");

    disp31.SetYTitle("[1]");
    disp32.SetYTitle("[1]");

    disp33.SetYTitle("Mean Abs. Time [FADC slice]");
    disp34.SetYTitle("RMS Abs. Time [FADC slices]");

    gStyle->SetOptStat(1111);
    gStyle->SetOptFit();
    
    // Charges
    TCanvas &c1 = display->AddTab("Fit.Charge");
    c1.Divide(2, 3);

    CamDraw(c1, disp1,calcam,1, 2 , 2);
    CamDraw(c1, disp2,calcam,2, 2 , 2);

    // Fit Probability
    TCanvas &c2 = display->AddTab("Fit.Prob");
    c2.Divide(1,3);

    CamDraw(c2, disp3,calcam,1, 1 , 4);

    // Reduced Sigmas
    TCanvas &c3 = display->AddTab("Red.Sigma");
    c3.Divide(2,3);

    CamDraw(c3, disp4,calcam,1, 2 , 2);
    CamDraw(c3, disp5,calcam,2, 2 , 2);

    // F-Factor Method
    TCanvas &c4 = display->AddTab("F-Factor");
    c4.Divide(3,3);

    CamDraw(c4, disp6,calcam,1, 3 , 2);
    CamDraw(c4, disp7,calcam,2, 3 , 2);
    CamDraw(c4, disp8,calcam,3, 3 , 2);

    // Blind Pixel Method
    TCanvas &c5 = display->AddTab("BlindPix");
    c5.Divide(3, 3);

    CamDraw(c5, disp9,calcam,1, 3 ,  9);
    CamDraw(c5, disp10,calcam,2, 3 , 2);
    CamDraw(c5, disp11,calcam,3, 3 , 2);

    // PIN Diode Method
    TCanvas &c6 = display->AddTab("PINDiode");
    c6.Divide(3,3);

    CamDraw(c6, disp12,calcam,1, 3 , 9);
    CamDraw(c6, disp13,calcam,2, 3 , 2);
    CamDraw(c6, disp14,calcam,3, 3 , 2);

    // Defects
    TCanvas &c7 = display->AddTab("Defects");
    c7.Divide(7,2);

    CamDraw(c7, disp15,calcam,1,7, 0);
    CamDraw(c7, disp16,calcam,2,7, 0);
    CamDraw(c7, disp17,calcam,3,7, 0);
    CamDraw(c7, disp18,calcam,4,7, 0);
    CamDraw(c7, disp19,calcam,5,7, 0);
    CamDraw(c7, disp20,calcam,6,7, 0);
    CamDraw(c7, disp21,calcam,7,7, 0);

    // Valid flags
    TCanvas &c8 = display->AddTab("Validity");
    c8.Divide(4,2);

    CamDraw(c8, disp22,calcam,1,4,0);
    CamDraw(c8, disp23,calcam,2,4,0);
    CamDraw(c8, disp24,calcam,3,4,0);
    CamDraw(c8, disp25,calcam,4,4,0);


    // Pedestals
    TCanvas &c9 = display->AddTab("Pedestals");
    c9.Divide(2,3);

    CamDraw(c9,disp26,calcam,1,2,1);
    CamDraw(c9,disp27,calcam,2,2,2);

    // Rel. Times
    TCanvas &c10 = display->AddTab("Fitted Rel. Times");
    c10.Divide(3,3);

    CamDraw(c10,disp28,calcam,1,3,2);
    CamDraw(c10,disp29,calcam,2,3,2);
    CamDraw(c10,disp30,calcam,3,3,4);

    // Time Defects
    TCanvas &c11 = display->AddTab("Time Def.");
    c11.Divide(2,2);

    CamDraw(c11, disp31,calcam,1,2, 0);
    CamDraw(c11, disp32,calcam,2,2, 0);

    // Abs. Times
    TCanvas &c12 = display->AddTab("Abs. Times");
    c12.Divide(2,3);

    CamDraw(c12,disp33,calcam,1,2,2);
    CamDraw(c12,disp34,calcam,2,2,2);

#endif
}

void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
{

  c.cd(i);
  gPad->SetBorderMode(0);
  MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
  //  obj1->AddNotify(evt);
  
  c.cd(i+j);
  gPad->SetBorderMode(0);
  obj1->Draw();
  ((MHCamera*)obj1)->SetPrettyPalette();

  if (fit != 0)
    {
      c.cd(i+2*j);
      gPad->SetBorderMode(0);
      TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
      
//      obj2->Sumw2();
      obj2->Draw();
      obj2->SetBit(kCanDelete);

      const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
      const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
      const Double_t integ = obj2->Integral("width")/2.5066283;
      const Double_t mean  = obj2->GetMean();
      const Double_t rms   = obj2->GetRMS();
      const Double_t width = max-min;

      if (rms == 0. || width == 0. )
        return;
      
      switch (fit)
        {
        case 1:
          TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
          sgaus->SetBit(kCanDelete);
          sgaus->SetParNames("Area","#mu","#sigma");
          sgaus->SetParameters(integ/rms,mean,rms);
          sgaus->SetParLimits(0,0.,integ);
          sgaus->SetParLimits(1,min,max);
          sgaus->SetParLimits(2,0,width/1.5);
          obj2->Fit("sgaus","QLR");
          obj2->GetFunction("sgaus")->SetLineColor(kYellow);
          break;

        case 2:
          TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
          dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
          TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
          dgaus->SetBit(kCanDelete);
          dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
          dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
                               integ/width/2.,(max+mean)/2.,width/4.);
          // The left-sided Gauss 
          dgaus->SetParLimits(0,integ-1.5,integ+1.5);
          dgaus->SetParLimits(1,min+(width/10.),mean);
          dgaus->SetParLimits(2,0,width/2.);
          // The right-sided Gauss 
          dgaus->SetParLimits(3,0,integ);
          dgaus->SetParLimits(4,mean,max-(width/10.));
          dgaus->SetParLimits(5,0,width/2.);
          obj2->Fit("dgaus","QLRM");
          obj2->GetFunction("dgaus")->SetLineColor(kYellow);
          break;
          
        case 3:
          TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
          tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
          tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
          TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
          tgaus->SetBit(kCanDelete);
          tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
                             "A_{2}","#mu_{2}","#sigma_{2}",
                             "A_{3}","#mu_{3}","#sigma_{3}");
          tgaus->SetParameters(integ,(min+mean)/2,width/4.,
                               integ/width/3.,(max+mean)/2.,width/4.,
                               integ/width/3.,mean,width/2.);
          // The left-sided Gauss 
          tgaus->SetParLimits(0,integ-1.5,integ+1.5);
          tgaus->SetParLimits(1,min+(width/10.),mean);
          tgaus->SetParLimits(2,width/15.,width/2.);
          // The right-sided Gauss 
          tgaus->SetParLimits(3,0.,integ);
          tgaus->SetParLimits(4,mean,max-(width/10.));
          tgaus->SetParLimits(5,width/15.,width/2.);
          // The Gauss describing the outliers
          tgaus->SetParLimits(6,0.,integ);
          tgaus->SetParLimits(7,min,max);
          tgaus->SetParLimits(8,width/4.,width/1.5);
          obj2->Fit("tgaus","QLRM");
          obj2->GetFunction("tgaus")->SetLineColor(kYellow);
          break;
        case 4:
          obj2->Fit("pol0","Q");
          obj2->GetFunction("pol0")->SetLineColor(kYellow);
          break;
        case 9:
          break;
        default:
          obj2->Fit("gaus","Q");
          obj2->GetFunction("gaus")->SetLineColor(kYellow);
          break;
        }
      
        TArrayI s0(3);
        s0[0] = 6;
        s0[1] = 1;
        s0[2] = 2;

        TArrayI s1(3);
        s1[0] = 3;
        s1[1] = 4;
        s1[2] = 5;

        TArrayI inner(1);
        inner[0] = 0;

        // Just to get the right (maximum) binning
        TH1D *half[2];
        half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2");
        half[1] = obj1->ProjectionS(s0, inner, "Sector 3-4-5");

        for (int i=0; i<2; i++)
        {
            half[i]->SetLineColor(kRed+i);
            half[i]->SetDirectory(0);
            half[i]->SetBit(kCanDelete);
            half[i]->Draw("same");
        }

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


