/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Markus Gaug, 02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MJCalibration
//
/////////////////////////////////////////////////////////////////////////////
#include "MJCalibration.h"

#include <TF1.h>
#include <TFile.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TSystem.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MRunIter.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MEvtLoop.h"

#include "MHCamera.h"
#include "MGeomCam.h"

#include "MPedestalCam.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePINDiode.h"
#include "MCalibrationChargeBlindPix.h"
#include "MCalibrationChargeCalc.h"

#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MBadPixelsMerge.h"
#include "MBadPixelsCam.h"
#include "MExtractSignal.h"
#include "MExtractPINDiode.h"
#include "MExtractBlindPixel.h"
#include "MExtractSignal2.h"
#include "MFCosmics.h"
#include "MContinue.h"
#include "MFillH.h"

#include "MJCalibration.h"
#include "MStatusDisplay.h"

ClassImp(MJCalibration);
using namespace std;

MJCalibration::MJCalibration(const char *name, const char *title) : fRuns(0)
{
    fName  = name  ? name  : "MJCalibration";
    fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
}

void MJCalibration::DrawProjection(MHCamera *obj1, Int_t fit) const
{

  TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
  obj2->Draw();
  obj2->SetBit(kCanDelete);
  
  if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    {
      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;
      
      TArrayI outer(1);
      outer[0] = 1;
      
      // Just to get the right (maximum) binning
      TH1D *half[4];
      half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
      half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
      half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
      half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
      
      for (int i=0; i<4; i++)
        {
          half[i]->SetLineColor(kRed+i);
          half[i]->SetDirectory(0);
          half[i]->SetBit(kCanDelete);
          half[i]->Draw("same");
        }
    }
  
  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.5;
  const Double_t mean  = obj2->GetMean();
  const Double_t rms   = obj2->GetRMS();
  const Double_t width = max-min;
  
  const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";

  const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"
                               "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
  TF1 *f=0;
  switch (fit)
    {
    case 1:
      f = new TF1("sgaus", "gaus(0)", min, max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("Area", "#mu", "#sigma");
      f->SetParameters(integ/rms, mean, rms);
      f->SetParLimits(0, 0,   integ);
      f->SetParLimits(1, min, max);
      f->SetParLimits(2, 0,   width/1.5);
      
      obj2->Fit(f, "QLR");
      break;
      
    case 2:
      f = new TF1("dgaus",dgausformula.Data(),min,max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
      f->SetParameters(integ,(min+mean)/2.,width/4.,
                       integ/width/2.,(max+mean)/2.,width/4.);
      // The left-sided Gauss 
      f->SetParLimits(0,integ-1.5      , integ+1.5);
      f->SetParLimits(1,min+(width/10.), mean);
      f->SetParLimits(2,0              , width/2.);
      // The right-sided Gauss 
      f->SetParLimits(3,0   , integ);
      f->SetParLimits(4,mean, max-(width/10.));
      f->SetParLimits(5,0   , width/2.);
      obj2->Fit(f,"QLRM");
      break;
      
    case 3:
      f = new TF1("tgaus",tgausformula.Data(),min,max);
      f->SetLineColor(kYellow);
      f->SetBit(kCanDelete);
      f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
                     "A_{2}","#mu_{2}","#sigma_{2}",
                     "A_{3}","#mu_{3}","#sigma_{3}");
      f->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 
      f->SetParLimits(0,integ-1.5,integ+1.5);
      f->SetParLimits(1,min+(width/10.),mean);
      f->SetParLimits(2,width/15.,width/2.);
      // The right-sided Gauss 
      f->SetParLimits(3,0.,integ);
      f->SetParLimits(4,mean,max-(width/10.));
      f->SetParLimits(5,width/15.,width/2.);
      // The Gauss describing the outliers
      f->SetParLimits(6,0.,integ);
      f->SetParLimits(7,min,max);
      f->SetParLimits(8,width/4.,width/1.5);
      obj2->Fit(f,"QLRM");
      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;
    }
}

void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit)
{
    c.cd(x);
    gPad->SetBorderMode(0);
    MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
    obj1->AddNotify(&fCalibrationCam);

    c.cd(x+y);
    gPad->SetBorderMode(0);
    obj1->Draw();

    if (!fit)
        return;

    c.cd(x+2*y);
    gPad->SetBorderMode(0);
    DrawProjection(obj1, fit);
}


void MJCalibration::DisplayResult(MParList &plist)
{
    if (!fDisplay)
        return;

    //
    // Update display
    //
    TString title = fDisplay->GetTitle();
    title += "--  Calibration ";
    title += fRuns->GetRunsAsString();
    title += "  --";
    fDisplay->SetTitle(title);

    //
    // Get container from list
    //
    MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");

    // 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;FFactorPhe",     "Nr. of Phe's (F-Factor Method)");
    MHCamera disp7 (geomcam, "Cal;FFactorConv",    "Conversion Factor (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 (Blind Pixel Method)");
    MHCamera disp11(geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
    MHCamera disp12(geomcam, "Cal;PINDiodePh",     "Photons flux outside plexiglass (PIN Diode Method)");
    MHCamera disp13(geomcam, "Cal;PINDiodeConv",   "Conversion Factor (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;Oscillation",    "Oscillating Pixels");
    MHCamera disp19(geomcam, "Cal;Saturation",     "Pixels with saturated Hi Gain");

    // Fitted charge means and sigmas
    disp1.SetCamContent(fCalibrationCam,  0);
    disp1.SetCamError(  fCalibrationCam,  1);
    disp2.SetCamContent(fCalibrationCam,  2);
    disp2.SetCamError(  fCalibrationCam,  3);
    // Fit probabilities
    disp3.SetCamContent(fCalibrationCam,  4);

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

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

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

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

    // Pixels with defects
    disp15.SetCamContent(fCalibrationCam,27);
    disp16.SetCamContent(fCalibrationCam,28);
    disp17.SetCamContent(fBadPixels     , 9);
    disp18.SetCamContent(fBadPixels     , 8);

    // Lo Gain calibration
    disp19.SetCamContent(fCalibrationCam,31);

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

    disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Q} - RMS^{2}_{Ped}} [FADC Counts]");
    disp5.SetYTitle("Red.Sigma/<Q> [1]");

    disp6.SetYTitle("PhE [#]");
    disp7.SetYTitle("Conv.Factor [PhE/FADC units]");
    disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Q}/\\mu_{Q} [1]");

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

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

    disp15.SetYTitle("[1]");
    disp16.SetYTitle("[1]");
    disp17.SetYTitle("[1]");
    disp18.SetYTitle("[1]");

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

    // Charges
    TCanvas &c1 = fDisplay->AddTab("Fit.Charge");
    c1.Divide(2, 3);

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

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

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

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

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

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

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

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

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

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

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

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

    CamDraw(c7, 1, 4, disp15, 0);
    CamDraw(c7, 2, 4, disp16, 0);
    CamDraw(c7, 3, 4, disp17, 0);
    CamDraw(c7, 4, 4, disp18, 0);

    // Lo Gain Calibration
    TCanvas &c8 = fDisplay->AddTab("LowGain");
    c8.Divide(1,3);

    CamDraw(c8, 1, 1, disp19, 0);
}

Bool_t MJCalibration::WriteResult()
{
    if (fOutputPath.IsNull())
        return kTRUE;

    const TString oname(GetOutputFile());

    *fLog << inf << "Writing to file: " << oname << endl;

    TFile file(oname, "UPDATE");

    if (fDisplay && fDisplay->Write()<=0)
    {
        *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
        return kFALSE;
    }

    if (fCalibrationCam.Write()<=0)
    {
        *fLog << err << "Unable to write MCalibrationCam to " << oname << endl;
        return kFALSE;
    }

    if (fBadPixels.Write()<=0)
    {
        *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
        return kFALSE;
    }

    return kTRUE;

}

void MJCalibration::SetOutputPath(const char *path)
{
    fOutputPath = path;
    if (fOutputPath.EndsWith("/"))
        fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
}

Bool_t MJCalibration::Process(MPedestalCam &pedcam)
{
    if (!ReadCalibrationCam())
        return ProcessFile(pedcam);

    return kTRUE;
}

TString MJCalibration::GetOutputFile() const
{
    if (!fRuns)
        return "";

    return Form("%s/%s-F1.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
}

Bool_t MJCalibration::ReadCalibrationCam()
{
    const TString fname = GetOutputFile();

    if (gSystem->AccessPathName(fname, kFileExists))
    {
        *fLog << err << "Input file " << fname << " doesn't exist." << endl;
        return kFALSE;
    }

    *fLog << inf << "Reading from file: " << fname << endl;

    TFile file(fname, "READ");
    if (fCalibrationCam.Read()<=0)
    {
        *fLog << "Unable to read MCalibrationCam from " << fname << endl;
        return kFALSE;
    }

    if (file.FindKey("MBadPixelsCam"))
    {
        MBadPixelsCam bad;
        if (bad.Read()<=0)
        {
            *fLog << "Unable to read MBadPixelsCam from " << fname << endl;
            return kFALSE;
        }
        fBadPixels.Merge(bad);
    }

    if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
        fDisplay->Read();

    return kTRUE;
}


Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
{
    if (!fRuns)
    {
        *fLog << err << "No Runs choosen... abort." << endl;
        return kFALSE;
    }
    if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
    {
        *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
        return kFALSE;
    }

    *fLog << inf;
    fLog->Separator(GetDescriptor());
    *fLog << "Calculate MCalibrationCam from Runs " << fRuns->GetRunsAsString() << endl;
    *fLog << endl;

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

    //
    // As long, as we don't have digital modules, 
    // we have to set the color of the pulser LED by hand
    //
    MCalibrationChargePINDiode pindiode;
    pindiode.SetColor(kCT1);

    MCalibrationChargeBlindPix blindpix;
    blindpix.SetColor(kCT1);

    // Setup Tasklist
    MParList plist;
    plist.AddToList(&pedcam);
    plist.AddToList(&fCalibrationCam);
    plist.AddToList(&pindiode);
    plist.AddToList(&blindpix);

    MTaskList tlist;
    plist.AddToList(&tlist);

    MGeomApply               apply;
    MBadPixelsMerge          merge(&fBadPixels);
    MExtractPINDiode         pinext;
    MExtractBlindPixel       blindext;
    // MExtractSignal   extract; // Do not use this at the moment...
    MExtractSignal2          extract;
    MCalibrationChargeCalc   calcalc;

    MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
    MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
    MFillH fillcam  ("MHCalibrationChargeCam",      "MExtractedSignalCam");
    // 
    // Apply a filter against cosmics
    // (will have to be needed in the future
    // when the calibration hardware-trigger is working)
    // 
    MFCosmics cosmics;
    MContinue cont(&cosmics);

    //calcalc.SkipBlindPixelFit();

    tlist.AddToList(&read);
    tlist.AddToList(&merge);
    tlist.AddToList(&apply);
    tlist.AddToList(&extract);
    tlist.AddToList(&pinext);  
    tlist.AddToList(&blindext);
    tlist.AddToList(&cont);
    tlist.AddToList(&fillcam);
    tlist.AddToList(&fillpin);
    tlist.AddToList(&fillblind);
    tlist.AddToList(&calcalc);

    // Create and setup the eventloop
    MEvtLoop evtloop(fName);
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(fDisplay);
    evtloop.SetLogStream(fLog);

    // Execute first analysis
    if (!evtloop.Eventloop())
    {
        *fLog << err << GetDescriptor() << ": Failed." << endl;
        return kFALSE;
    }

    tlist.PrintStatistics();

    DisplayResult(plist);

    if (!WriteResult())
        return kFALSE;

    *fLog << inf << GetDescriptor() << ": Done." << endl;

    return kTRUE;
}