/* ======================================================================== *\
!
! *
! * 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, 8/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Wolfgang Wittek, 1/2002 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MHEffectiveOnTime
//
//  Fills a 2-D histogram Time-difference vs. Theta
//
//  From this histogram the effective on-time is determined by a fit.
//  The result of the fit (see Fit()) and the fit-parameters (like chi^2)
//  are stored in corresponding histograms
//
//////////////////////////////////////////////////////////////////////////////

#include "MHEffectiveOnTime.h"

#include <TF1.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TMinuit.h>
#include <TRandom.h>

#include "MTime.h"
#include "MParameters.h"
#include "MPointingPos.h"

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MHEffectiveOnTime);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor. It sets name and title of the histogram.
//
MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
    : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60),
    fNameProjDeltaT(Form("DeltaT_%p", this)), fNameProjTheta(Form("Theta_%p", this))
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHEffectiveOnTime";
    fTitle = title ? title : "2-D histogram in Theta and time difference";

    // Main histogram
    fHDeltaT.SetName("DeltaT");
    fHDeltaT.SetXTitle("\\Delta t [s]");
    fHDeltaT.SetYTitle("\\Theta [\\circ]");
    fHDeltaT.UseCurrentStyle();
    fHDeltaT.SetDirectory(NULL);

    // effective on time versus theta
    fHEffOnTheta.SetName("EffOnTime");
    fHEffOnTheta.SetTitle("Effective On Time T_{eff}");
    fHEffOnTheta.SetXTitle("\\Theta [\\circ]");
    fHEffOnTheta.SetYTitle("T_{eff} [s]");
    fHEffOnTheta.UseCurrentStyle();
    fHEffOnTheta.SetDirectory(NULL);
    fHEffOnTheta.GetYaxis()->SetTitleOffset(1.2);
    //fHEffOn.Sumw2();

    // effective on time versus time
    fHEffOnTime.SetName("EffOnTime");
    fHEffOnTime.SetTitle("Effective On Time T_{eff}");
    fHEffOnTime.SetXTitle("Time");
    fHEffOnTime.SetYTitle("T_{eff} [s]");
    fHEffOnTime.UseCurrentStyle();
    fHEffOnTime.SetDirectory(NULL);
    fHEffOnTime.GetYaxis()->SetTitleOffset(1.2);
    fHEffOnTime.GetXaxis()->SetLabelSize(0.033);
    fHEffOnTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    fHEffOnTime.GetXaxis()->SetTimeDisplay(1);
    fHEffOnTime.Sumw2();

    // chi2/NDF versus theta
    fHChi2.SetName("Chi2/NDF");
    fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit");
    fHChi2.SetXTitle("\\Theta [\\circ]");
    fHChi2.SetYTitle("\\chi^{2}/NDF");
    fHChi2.UseCurrentStyle();
    fHChi2.SetDirectory(NULL);
    fHChi2.GetYaxis()->SetTitleOffset(1.2);
    //fHChi2.Sumw2();

    // chi2 probability
    fHProb.SetName("Prob");
    fHProb.SetTitle("\\chi^{2} Probability of Fit");
    fHProb.SetXTitle("\\Theta [\\circ]");
    fHProb.SetYTitle("p [%]");
    fHProb.UseCurrentStyle();
    fHProb.SetDirectory(NULL);
    fHProb.GetYaxis()->SetTitleOffset(1.2);
    fHProb.SetMaximum(101);
    //fHChi2.Sumw2();

    // lambda versus theta
    fHLambda.SetName("lambda");
    fHLambda.SetTitle("lambda of Effective On Time Fit");
    fHLambda.SetXTitle("\\Theta [\\circ]");
    fHLambda.SetYTitle("\\lambda [Hz]");
    fHLambda.UseCurrentStyle();
    fHLambda.SetDirectory(NULL);
    // fHLambda.Sumw2();

    // N0 versus theta
    fHN0.SetName("N0");
    fHN0.SetTitle("Ideal number of events N_{0}");
    fHN0.SetXTitle("\\Theta [\\circ]");
    fHN0.SetYTitle("N_{0}");
    fHN0.UseCurrentStyle();
    fHN0.SetDirectory(NULL);
    //fHN0del.Sumw2();

    // setup binning
    MBinning btheta("BinningTheta");
    btheta.SetEdgesCos(101, 0, 60);

    MBinning btime("BinningDeltaT");
    btime.SetEdges(50, 0, 0.1);

    MH::SetBinning(&fHDeltaT, &btime, &btheta);

    btheta.Apply(fHEffOnTheta);
    btheta.Apply(fHChi2);
    btheta.Apply(fHLambda);
    btheta.Apply(fHN0);
    btheta.Apply(fHProb);
}

// --------------------------------------------------------------------------
//
// Set the binnings and prepare the filling of the histogram
//
Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist)
{
   fPointPos = (MPointingPos*)plist->FindObject("MPointingPos");
   if (!fPointPos)
   {
       *fLog << err << dbginf << "MPointingPos not found... aborting." << endl;
       return kFALSE;
   }

   // FIXME: Remove const-qualifier from base-class!
   fTime = (MTime*)const_cast<MParList*>(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime");
   if (!fTime)
       return kFALSE;
   fParam = (MParameterDerr*)const_cast<MParList*>(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime");
   if (!fParam)
       return kFALSE;

   const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningDeltaT");
   const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta");
   if (!binstheta || !binsdtime)
       *fLog << warn << dbginf << "At least one MBinning not found... ignored." << endl;
   else
   {
       SetBinning(&fHDeltaT, binsdtime, binstheta);

       binstheta->Apply(fHEffOnTheta);
       binstheta->Apply(fHChi2);
       binstheta->Apply(fHLambda);
       binstheta->Apply(fHN0);
       binstheta->Apply(fHProb);
   }

   return kTRUE;
}

Bool_t MHEffectiveOnTime::Finalize()
{
    FitThetaBins();
    Calc();

    fIsFinalized = kTRUE;

    return kTRUE;
}

void MHEffectiveOnTime::FitThetaBins()
{
    const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));

    // nbins = number of Theta bins
    const Int_t nbins = fHDeltaT.GetNbinsY();

    TH1D *h=0;
    for (int i=1; i<=nbins; i++)
    {
        //        TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
        h = fHDeltaT.ProjectionX(name, i, i, "E");

        Double_t res[7];
        if (!FitH(h, res))
            continue;

        // the effective on time is Nm/lambda
        fHEffOnTheta.SetBinContent(i, res[0]);
        fHEffOnTheta.SetBinError  (i, res[1]);

        // plot chi2-probability of fit
        fHProb.SetBinContent(i, res[2]);

        // plot chi2/NDF of fit
        fHChi2.SetBinContent(i, res[3]);

        // lambda of fit
        fHLambda.SetBinContent(i, res[4]);
        fHLambda.SetBinError  (i, res[5]);

        // N0 of fit
        fHN0.SetBinContent(i, res[6]);

        // Rdead (from fit) is the fraction from real time lost by the dead time
        //fHRdead.SetBinContent(i, Rdead);
        //fHRdead.SetBinError  (i,dRdead);
    }

    // Histogram is reused via gROOT->FindObject()
    // Need to be deleted only once
    if (h)
        delete h;
}

Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const
{
    const Double_t Nm = h->Integral();

    // FIXME: Do fit only if contents of bin has changed
    if (Nm<=0)
        return kFALSE;

    // determine range (yq[0], yq[1]) of time differences
    // where fit should be performed;
    // require a fraction >=xq[0] of all entries to lie below yq[0]
    //     and a fraction <=xq[1] of all entries to lie below yq[1];
    // within the range (yq[0], yq[1]) there must be no empty bin;
    // choose pedestrian approach as long as GetQuantiles is not available
    Double_t xq[2] = { 0.05, 0.95 };
    Double_t yq[2];
    h->GetQuantiles(2, yq, xq);

    // Nmdel = Nm * binwidth,  with Nm = number of observed events
    const Double_t Nmdel = h->Integral("width");

    //
    // Setup Poisson function for the fit:
    // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
    //
    // parameter 0 = lambda
    // parameter 1 = N0*del
    //
    TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
    func.SetParNames("lambda", "N0", "del");

    func.SetParameter(0, 100);       // Hz
    func.SetParameter(1, Nm);
    func.FixParameter(2, Nmdel/Nm);

    // options : 0  do not plot the function
    //           I  use integral of function in bin rather than value at bin center
    //           R  use the range specified in the function range
    //           Q  quiet mode
    h->Fit(&func, "0IRQ");

    const Double_t chi2   = func.GetChisquare();
    const Int_t    NDF    = func.GetNDF();

    // was fit successful ?
    if (NDF<=0 || chi2>=2.5*NDF)
        return kFALSE;

    const Double_t lambda = func.GetParameter(0);
    const Double_t N0     = func.GetParameter(1);
    const Double_t prob   = func.GetProb();

    /*
     *fLog << all << "Nm/lambda=" << Nm/lambda << "  chi2/NDF=";
     *fLog << (NDF ? chi2/NDF : 0.0) << "  lambda=";
     *fLog << lambda << "  N0=" << N0 << endl;
     */

    Double_t emat[2][2];
    gMinuit->mnemat((Double_t*)emat, 2);

    const Double_t dldl   = emat[0][0];
    //const Double_t dN0dN0 = emat[1][1];

    const Double_t teff   = Nm/lambda;
    const Double_t dteff  = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);

    const Double_t dl     = TMath::Sqrt(dldl);

    //const Double_t kappa  = Nm/N0;
    //const Double_t Rdead  = 1.0 - kappa;
    //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);

    // the effective on time is Nm/lambda
    res[0] = teff;
    res[1] = dteff;

    // plot chi2-probability of fit
    res[2] = prob*100;

    // plot chi2/NDF of fit
    res[3] = NDF ? chi2/NDF : 0.0;

    // lambda of fit
    res[4] = lambda;
    res[5] = dl;

    // N0 of fit
    res[6] = N0;

    // Rdead (from fit) is the fraction from real time lost by the dead time
    //fHRdead.SetBinContent(i, Rdead);
    //fHRdead.SetBinError  (i,dRdead);

    return kTRUE;
}

void MHEffectiveOnTime::Paint(Option_t *opt)
{
    TVirtualPad *padsave = gPad;

    TString o(opt);
    if (o==(TString)"fit")
    {
        TH1D *h0=0;

        padsave->GetPad(1)->cd(1);
        if ((h0 = (TH1D*)gPad->FindObject(fNameProjDeltaT)))
        {
            h0 = fHDeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
            if (h0->GetEntries()>0)
                gPad->SetLogy();
        }

        padsave->GetPad(2)->cd(1);
        if ((h0 = (TH1D*)gPad->FindObject(fNameProjTheta)))
            fHDeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");

        if (!fIsFinalized)
            FitThetaBins();
    }

    TH1D *h=0;
    if (o==(TString)"theta")
        h = &fHEffOnTheta;
    if (o==(TString)"time")
        h = &fHEffOnTime;

    if (h)
    {
        Double_t error = 0;
        for (int i=0; i<h->GetXaxis()->GetNbins(); i++)
            error += h->GetBinError(i);

        TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", h->Integral(), error));
        text.SetBit(TLatex::kTextNDC);
        text.SetTextSize(0.04);
        text.Paint();
    }
    gPad = padsave;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHEffectiveOnTime::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("fit");

    pad->Divide(2, 1, 0, 0);

    TH1 *h;

    pad->cd(1);
    gPad->SetBorderMode(0);
    gPad->Divide(1, 2, 0, 0);
    pad->GetPad(1)->cd(1);
    gPad->SetBorderMode(0);
    h = fHDeltaT.ProjectionX(fNameProjDeltaT, -1, 9999, "E");
    h->SetTitle("Distribution of \\Delta t [s]");
    h->SetXTitle("\\Delta t [s]");
    h->SetYTitle("Counts");
    h->SetDirectory(NULL);
    h->SetMarkerStyle(kFullDotMedium);
    h->SetBit(kCanDelete);
    h->Draw();

    pad->GetPad(1)->cd(2);
    gPad->SetBorderMode(0);
    fHEffOnTime.Draw();
    AppendPad("time");

    pad->cd(2);
    gPad->SetBorderMode(0);
    gPad->Divide(1, 3, 0, 0);

    pad->GetPad(2)->cd(1);
    gPad->SetBorderMode(0);
    h = fHDeltaT.ProjectionY(fNameProjTheta, -1, 9999, "E");
    h->SetTitle("Distribution of  \\Theta [\\circ]");
    h->SetXTitle("\\Theta [\\circ]");
    h->SetYTitle("Counts");
    h->SetDirectory(NULL);
    h->SetMarkerStyle(kFullDotMedium);
    h->SetBit(kCanDelete);
    h->GetYaxis()->SetTitleOffset(1.1);
    h->Draw();

    pad->GetPad(2)->cd(2);
    gPad->SetBorderMode(0);
    fHProb.Draw();

    pad->GetPad(2)->cd(3);
    gPad->SetBorderMode(0);
    fHEffOnTheta.Draw();
    AppendPad("theta");
}

void MHEffectiveOnTime::FillTimeBin(Double_t val, Double_t err)
{
    // Get x-axis
    TAxis &x = *fHEffOnTime.GetXaxis();

    // Get number of bins
    const Int_t n = x.GetNbins();

    // Fill last bin with new values
    fHEffOnTime.SetBinContent(n, fParam->GetVal());
    fHEffOnTime.SetBinError(n, fParam->GetErr());

    // Get time range of histogram, get bin width
    const Double_t lo = fHEffOnTime.GetXaxis()->GetXmin();
    const Double_t up = fHEffOnTime.GetXaxis()->GetXmax();
    const Double_t w  = fHEffOnTime.GetXaxis()->GetBinWidth(1);

    // Enhance histogram by one bin
    MBinning bins;
    bins.SetEdges(n+2, lo, up+w);
    bins.Apply(fHEffOnTime);

    // Transform overflow bin
    // fHEffOnTime.SetBinContent(n+2, fHEffOnTime.GetBinContent(n+1));
    // fHEffOnTime.SetBinError(n+2, fHEffOnTime.GetBinError(n+1));
}

void MHEffectiveOnTime::Calc()
{
    TH1D *h = fHDeltaT.ProjectionX("", -1, 99999, "E");
    h->SetDirectory(0);

    Double_t res[7];
    Bool_t rc = FitH(h, res);
    delete h;

    if (!rc)
        return;

    const Double_t val   = res[0];
    const Double_t error = res[1];

    fParam->SetVal(val>fEffOnTime0?val-fEffOnTime0:0, error-fEffOnErr0);
    fParam->SetReadyToSave();

    FillTimeBin(fParam->GetVal(), fParam->GetErr());

    fEffOnTime0 = val;
    fEffOnErr0  = error;

    MTime now(*fTime);
    now.AddMilliSeconds(fInterval*1000);

    *fLog <<all << now << " - ";// << fLastTime-fTime;
    *fLog << Form("T_{eff} = %.1fs \\pm %.1fs",
                  fParam->GetVal(), fParam->GetErr());
    *fLog << Form("   %.1f   %.1f",
                  val, error) << endl;

    fTime->AddMilliSeconds(fInterval*1000);
}

// --------------------------------------------------------------------------
//
//  Fill the histogram
//
Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w)
{
    const MTime *time = dynamic_cast<const MTime*>(par);
    if (!time)
    {
        *fLog << err << "ERROR - MHEffectiveOnTime::Fill without argument or container doesn't inherit from MTime... abort." << endl;
        return kFALSE;
    }

    if (*fTime==MTime())
    {
        *fLog << all << *time << " - " << *fTime << " " << TMath::Floor(*time) << endl;
        fEffOnTime0 = 0;
        fEffOnErr0  = 0;

        fParam->SetVal(0, 0);
        fParam->SetReadyToSave();

        *fTime = *time;

        MBinning bins;
        bins.SetEdges(2, fTime->GetAxisTime(), fTime->GetAxisTime()+fInterval);
        bins.Apply(fHEffOnTime);

        // Make this just a ns before the first event
        fTime->Minus1ns();
    }

    if (*fTime+fInterval<*time)
    {
        FitThetaBins();
        Calc();
    }

    fHDeltaT.Fill(*time-fLastTime, fPointPos->GetZd(), w);
    fLastTime = *time;

    return kTRUE;
}

