/* ======================================================================== *\
!
! *
! * 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>
!
!   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 "MPedestalCam.h"
#include "MCalibrationCam.h"

#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MExtractSignal.h"
#include "MCalibrationCalc.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();
    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.5;
    const Double_t mean  = obj2->GetMean();
    const Double_t rms   = obj2->GetRMS();
    const Double_t width = max-min;

    TF1 *f=0;
    switch (fit)
    {
    case 0:
        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 1:
        f = new TF1("dgaus", "gaus(0)+gaus(3)", min, max);
        f->SetLineColor(kYellow);
        f->SetBit(kCanDelete);
        f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2");
        f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4);
        f->SetParLimits(0, 0,   integ);
        f->SetParLimits(1, min, max);
        f->SetParLimits(2, 0,   width/3);
        f->SetParLimits(3, 0,   integ);
        f->SetParLimits(4, min, max);
        f->SetParLimits(5, 0,   width/3);

        obj2->Fit(f, "QLR");
        break;

    case 2:
        f = new TF1("tgaus", "gaus(0)+gaus(3)+gaus(6)", min, max);
        f->SetLineColor(kYellow);
        f->SetBit(kCanDelete);
        f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2", "A3", "#mu3", "#sigma3");
        f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4, integ/width,min+width/6, width/2);
        f->SetParLimits(0, 0,   integ);
        f->SetParLimits(1, min, max);
        f->SetParLimits(2, 0,   width/4);
        f->SetParLimits(3, 0,   integ);
        f->SetParLimits(4, min, max);
        f->SetParLimits(5, 0,   width/4);
        f->SetParLimits(6, 0,   integ);
        f->SetParLimits(7, min, max);
        f->SetParLimits(8, 0,   width/2);

        obj2->Fit(f, "QLR");
        break;

    case 3:
        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, Int_t x, Int_t y, MHCamera &cam1, 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();

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

    // tbretz: Using gStyle in this context is totally useless!
    //const Float_t he = gStyle->GetStatH();
    //const Float_t wi = gStyle->GetStatH();
    //gStyle->SetStatH(0.4);
    //gStyle->SetStatW(0.25);
    //gStyle->SetStatH(he);
    //gStyle->SetStatW(wi);
}

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 disp3 (geomcam, "Cal;SigmaCharge",   "Sigma of Fitted Charges");
    MHCamera disp5 (geomcam, "Cal;FitProb",       "Probability of Fit");
    /*
     MHCamera disp6 (geomcam, "Cal;Time",          "Arrival Times");
     MHCamera disp7 (geomcam, "Cal;SigmaTime",     "Sigma of Arrival Times");
     MHCamera disp8 (geomcam, "Cal;TimeChiSquare", "Chi Square of Time Fit");
     */
//    MHCamera disp9 (geomcam, "Cal;Ped",           "Pedestals");
//    MHCamera disp10(geomcam, "Cal;PedRms",        "Pedestal RMS");
//    MHCamera disp11(geomcam, "Cal;RSigma",         "Reduced Sigmas");
    MHCamera disp12(geomcam, "Cal;FFactorPhe",     "Nr. of Phe's (F-Factor Method)");
    MHCamera disp13(geomcam, "Cal;FFactorConv",    "Conversion Factor (F-Factor Method)");
    MHCamera disp14(geomcam, "Cal;BlindPixPhe",    "Nr. of Photons inside plexiglass (Blind Pixel Method)");
    MHCamera disp15(geomcam, "Cal;BlindPixConv",   "Conversion Factor (Blind Pixel Method)");
//    MHCamera disp16(geomcam, "Cal;RSigma/Charge",  "Reduced Sigma per Charge");

    disp1.SetCamContent(fCalibrationCam, 0);
    disp1.SetCamError(fCalibrationCam, 1);

    disp3.SetCamContent(fCalibrationCam, 2);
    disp3.SetCamError(fCalibrationCam, 3);

    disp5.SetCamContent(fCalibrationCam, 4);

    /*
    disp6.SetCamContent(fCalibrationCam, 5);
    disp6.SetCamError(fCalibrationCam, 6);
    disp7.SetCamContent(fCalibrationCam, 6);
    disp8.SetCamContent(fCalibrationCam, 7);
    */
/*
    disp9.SetCamContent(calcam, 8);
    disp9.SetCamError(calcam, 9);
    disp10.SetCamContent(calcam, 9);
    disp11.SetCamContent(fCalibrationCam, 10);
 */
    disp12.SetCamContent(fCalibrationCam, 11);
    disp12.SetCamError(fCalibrationCam, 12);
    disp13.SetCamContent(fCalibrationCam, 13);
    disp13.SetCamError(fCalibrationCam, 14);

    disp14.SetCamContent(fCalibrationCam, 15);
    disp15.SetCamContent(fCalibrationCam, 16);
 //   disp16.SetCamContent(fCalibrationCam, 17);

    disp1.SetYTitle("Charge [FADC units]");
    disp3.SetYTitle("\\sigma_{Charge} [FADC units]");
    disp5.SetYTitle("P_{Charge} [1]");
    /*
     disp6.SetYTitle("Arr. Time [FADC slice Nr.]");
     disp7.SetYTitle("\\sigma_{Time} [FADC slices]");
     disp8.SetYTitle("\\chi^{2}_{Time} [1]");
    disp9.SetYTitle("Ped [FADC Counts ]");
    disp10.SetYTitle("RMS_{Ped} [FADC Counts ]");
    disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC units]");
     */
    disp12.SetYTitle("Nr. Photo-Electrons [1]");
    disp13.SetYTitle("Conversion Factor [PhE/FADC Count]");
    disp14.SetYTitle("Nr. Photons [1]");
    disp15.SetYTitle("Conversion Factor [Phot/FADC Count]");
//    disp16.SetYTitle("Reduced Sigma / Charge [1]");

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

    TCanvas &c1 = fDisplay->AddTab("FitCharge");
    c1.Divide(2, 3);

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

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

    CamDraw(c2, 1, 1, disp5, 3);
/*
    // Times
    TCanvas &c3 = fDisplay->AddTab("FitTimes");
    c3.Divide(3,3);

    CamDraw(c3, 1, 3, disp6, 1);
    CamDraw(c3, 2, 3, disp7, 0);
    CamDraw(c3, 3, 3, disp8, 0);

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

    CamDraw(c4, 1, 2, disp9,  0);
    CamDraw(c4, 2, 2, disp10, 1);

    // Reduced Sigmas
    TCanvas &c5 = fDisplay->AddTab("RedSigma");
    c5.Divide(2,3);

    CamDraw(c5, 1, 2, disp11, 2);
    CamDraw(c5, 2, 2, disp16, 2);
    */

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

    CamDraw(c6, 1, 2, disp12, 1);
    CamDraw(c6, 2, 2, disp13, 1);

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

    CamDraw(c7, 1, 2, disp14, 9);
    CamDraw(c7, 2, 2, disp15, 1);

}

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;
    }
    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 (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 metch 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);

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

    MTaskList tlist;
    plist.AddToList(&tlist);

    MGeomApply       apply;
    MExtractSignal   extract;
    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);
    //calcalc.SkipBlindPixelFit();

    tlist.AddToList(&apply);
    tlist.AddToList(&read);
    tlist.AddToList(&extract);
    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;
}
