/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//const TString pedfile = "/remote/home/pc2/operator/Crab20040214/20040215_16743_P_CrabOn_E.root";
const TString pedfile = "../20040215_16770_P_OffCrab4_E.root";

//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";

void pedestalstudies(TString pedname=pedfile)
{

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

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

    //
    // Create a empty Parameter List and an empty Task List
    // The tasklist is identified in the eventloop by its name
    //
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);

    //
    // Now setup the tasks and tasklist for the pedestals:
    // ---------------------------------------------------
    //

    MReadMarsFile read("Events", pedname);
    read.DisableAutoScheme();

    MGeomApply      geomapl;
    MExtractSignal  sigcalc;

    //
    // Set the extraction range higher:
    //		
    //sigcalc.SetRange(1,14,1,14);

    MPedCalcPedRun pedcalc;

    //
    // Additionally to calculating the pedestals, 
    // you can fill histograms and look at them
    //
    MFillH fill("MHPedestalCam", "MExtractedSignalCam");

    tlist.AddToList(&read);
    tlist.AddToList(&geomapl);
    tlist.AddToList(&sigcalc);
    tlist.AddToList(&pedcalc);
    tlist.AddToList(&fill);

    MGeomCamMagic  geomcam;
    MPedestalCam   pedcam;
    MHPedestalCam  hpedcam;
    plist.AddToList(&geomcam);
    plist.AddToList(&pedcam);
    plist.AddToList(&hpedcam);

    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;

    evtloop.SetParList(&plist);
    evtloop.SetDisplay(display);
 
    //
    // Execute first analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist.PrintStatistics();

    // 
    // Look at one specific pixel, after all the histogram manipulations:
    //
//    hpedcam[9].DrawClone("fourierevents");


    MHCamera dispped0  (geomcam, "Ped;Pedestal",               "Mean per Slice");
    MHCamera dispped1  (geomcam, "Ped;PedestalErr",            "Mean Error per Slice");
    MHCamera dispped2  (geomcam, "Ped;PedestalRms",            "RMS per Slice");
    MHCamera dispped3  (geomcam, "Ped;PedestalRmsErr",         "RMS Error per Slice");

    MHCamera dispped4  (geomcam, "Ped;Mean",                   "Fitted Mean per Slice");
    MHCamera dispped5  (geomcam, "Ped;MeanErr",                "Fitted Error of Mean per Slice");
    MHCamera dispped6  (geomcam, "Ped;Sigma",                  "Fitted Sigma per Slice");
    MHCamera dispped7  (geomcam, "Ped;SigmaErr",               "Fitted Error of Sigma per Slice");
    MHCamera dispped8  (geomcam, "Ped;Prob",                   "Probability of Fit");
    MHCamera dispped9  (geomcam, "Ped;DeltaPedestalMean",      "Rel. Diff. Mean per Slice (Calc.-Fitte)");
    MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
    MHCamera dispped11  (geomcam, "Ped;DeltaRmsSigma",         "Rel. Diff. RMS per Slice (Calc.-Fitted)");
    MHCamera dispped12  (geomcam, "Ped;DeltaRmsSigmaError",    "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
    MHCamera dispped13  (geomcam, "Ped;FitOK",                 "Gaus Fit not OK");
    MHCamera dispped14  (geomcam, "Ped;FourierOK",             "Fourier Analysis not OK");

    dispped0.SetCamContent(  pedcam, 0);
    dispped0.SetCamError(    pedcam, 1);
    dispped1.SetCamContent(  pedcam, 1);
    dispped2.SetCamContent(  pedcam, 2);
    dispped2.SetCamError(    pedcam, 3);
    dispped3.SetCamContent(  pedcam, 3);

    dispped4.SetCamContent( hpedcam, 0);
    dispped4.SetCamError(   hpedcam, 1);
    dispped5.SetCamContent( hpedcam, 1);
    dispped6.SetCamContent( hpedcam, 2);
    dispped6.SetCamError(   hpedcam, 3);
    dispped7.SetCamContent( hpedcam, 3);
    dispped8.SetCamContent( hpedcam, 4);
    dispped9.SetCamContent( hpedcam, 5);
    dispped9.SetCamError(   hpedcam, 6);
    dispped10.SetCamContent(hpedcam, 7);
    dispped11.SetCamContent(hpedcam, 8);
    dispped11.SetCamError(  hpedcam, 9);
    dispped12.SetCamContent(hpedcam, 10);
    dispped13.SetCamContent(hpedcam, 11);
    dispped14.SetCamContent(hpedcam, 12);

    dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
    dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
    dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
    dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
    dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
    dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
    dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
    dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
    dispped8.SetYTitle("Fit Probability [1]");
    dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
    dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
    dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
    dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
    dispped13.SetYTitle("[1]");
    dispped14.SetYTitle("[1]");
    
    // Histogram values
    TCanvas &b1 = display->AddTab("Ped.Calc."); 
    b1.Divide(4,3);

    CamDraw(b1,dispped0,pedcam,1,4,1);
    CamDraw(b1,dispped1,pedcam,2,4,2);
    CamDraw(b1,dispped2,pedcam,3,4,2);
    CamDraw(b1,dispped3,pedcam,4,4,2);

    // Fitted values 
    TCanvas &b2 = display->AddTab("Ped.Fit"); 
    b2.Divide(4,3);

    CamDraw(b2,dispped4,hpedcam,1,4,1);
    CamDraw(b2,dispped5,hpedcam,2,4,2);
    CamDraw(b2,dispped6,hpedcam,3,4,2);
    CamDraw(b2,dispped7,hpedcam,4,4,2);


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

    CamDraw(b3,dispped8,hpedcam,1,1,3);

    // Differences
    TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
    c4.Divide(4,3);

    CamDraw(c4,dispped9,hpedcam,1,4,1);
    CamDraw(c4,dispped10,hpedcam,2,4,1);
    CamDraw(c4,dispped11,hpedcam,3,4,1);
    CamDraw(c4,dispped12,hpedcam,4,4,1);

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

    CamDraw(c5,dispped13,hpedcam,1,2,0);
    CamDraw(c5,dispped14,hpedcam,2,2,0);

}

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();
      
//      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;
        }
      
      gPad->Modified();
      gPad->Update();
      
    }
}
