/* ======================================================================== *\
!
! *
! * 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, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHAlpha
//
// Create a single Alpha-Plot. The alpha-plot is fitted online. You can
// check the result when it is filles in the MStatusDisplay
// For more information see MHFalseSource::FitSignificance
//
// For convinience (fit) the output significance is stored in a
// container in the parlisrt
//
// PRELIMINARY!
//
//////////////////////////////////////////////////////////////////////////////
#include "MHAlpha.h"

#include <TF1.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TPaveStats.h>
#include <TStopwatch.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MHillasSrc.h"
#include "MEnergyEst.h"
#include "MTime.h"
#include "MObservatory.h"
#include "MPointingPos.h"
#include "MAstroSky2Local.h"
#include "MStatusDisplay.h"
#include "MParameters.h"
#include "MHMatrix.h"

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

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

ClassImp(MHAlpha);
ClassImp(MAlphaFitter);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHAlpha::MHAlpha(const char *name, const char *title)
    : fResult(0), /*fExcess(0),*/ fEnergy(0), fPointPos(0), fTimeEffOn(0),
    fTime(0), fNameProjAlpha(Form("Alpha_%p", this)), fMatrix(0)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHAlpha";
    fTitle = title ? title : "Alpha plot";

    fHAlpha.SetName("Alpha");
    fHAlpha.SetTitle("Alpha");
    fHAlpha.SetXTitle("\\Theta [deg]");
    fHAlpha.SetYTitle("E_{est} [GeV]");
    fHAlpha.SetZTitle("|\\alpha| [\\circ]");
    fHAlpha.SetDirectory(NULL);
    fHAlpha.UseCurrentStyle();

    // Main histogram
    fHAlphaTime.SetName("Alpha");
    fHAlphaTime.SetXTitle("|\\alpha| [\\circ]");
    fHAlphaTime.SetYTitle("Counts");
    fHAlphaTime.UseCurrentStyle();
    fHAlphaTime.SetDirectory(NULL);


    fHEnergy.SetName("Energy");
    fHEnergy.SetTitle(" N_{exc} vs. E_{est} ");
    fHEnergy.SetXTitle("E_{est} [GeV]");
    fHEnergy.SetYTitle("N_{exc}");
    fHEnergy.SetDirectory(NULL);
    fHEnergy.UseCurrentStyle();

    fHTheta.SetName("Theta");
    fHTheta.SetTitle(" N_{exc} vs. \\Theta ");
    fHTheta.SetXTitle("\\Theta [\\circ]");
    fHTheta.SetYTitle("N_{exc}");
    fHTheta.SetDirectory(NULL);
    fHTheta.UseCurrentStyle();

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

    MBinning binsa, binse, binst;
    binsa.SetEdges(18, 0, 90);
    binse.SetEdgesLog(25, 10, 100000);
    binst.SetEdgesCos(50, 0, 60);
    binse.Apply(fHEnergy);
    binst.Apply(fHTheta);
    binsa.Apply(fHAlphaTime);

    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);
}

void MHAlpha::FitEnergySpec(Bool_t paint)
{
    /*
    TF1 f("Spectrum", "pow(x, [0])*exp(-x/[1])*[2]*x");//
    f.SetParameter(0,  -2.1);
    f.SetParameter(1,  1400); // 50
    f.SetParameter(2,  fHEnergy.Integral()); // 50
    f.SetLineWidth(2);

    fHEnergy.Fit(&f,  "0QI", "", 100, 2400); // Integral? 90, 2800

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

    // was fit successful ?
    const Bool_t ok = NDF>0 && chi2<2.5*NDF;
    f.SetLineColor(ok ? kGreen : kRed);

    if (paint)
    {
        f.Paint("same");

        if (ok)
        {
        TString s;
        s += Form("\\alpha=%.1f ",   f.GetParameter(0));
        s += Form("E_{c}=%.1f TeV ", f.GetParameter(1)/1000);
        s += Form("N=%.1e",          f.GetParameter(2));

        TLatex text(0.25, 0.94, s.Data());
        text.SetBit(TLatex::kTextNDC);
        text.SetTextSize(0.04);
        text.Paint();
        }
    }*/
}

void MHAlpha::FitEnergyBins(Bool_t paint)
{
    // Do not store the 'final' result in fFit
    MAlphaFitter fit(fFit);

    const Int_t n = fHAlpha.GetNbinsY();

    for (int i=1; i<=n; i++)
    {
        TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", -1, 9999, i, i, i==1?"E":"");
        if (fit.Fit(*h))
        {
            fHEnergy.SetBinContent(i, fit.GetEventsExcess());
            fHEnergy.SetBinError(i, fit.GetEventsExcess()*0.2);
        }
        delete h;
    }
}

void MHAlpha::FitThetaBins(Bool_t paint)
{
    // Do not store the 'final' result in fFit
    MAlphaFitter fit(fFit);

    const Int_t n = fHAlpha.GetNbinsX();

    for (int i=1; i<=n; i++)
    {
        TH1D *h = fHAlpha.ProjectionZ("Alpha_EE", i, i, -1, 9999, i==1?"E":"");
        if (fit.Fit(*h))
        {
            fHTheta.SetBinContent(i, fit.GetEventsExcess());
            fHTheta.SetBinError(i, fit.GetEventsExcess()*0.2);
        }
        delete h;
    }
}

Bool_t MHAlpha::SetupFill(const MParList *pl)
{
    fHAlpha.Reset();
    fHEnergy.Reset();
    fHTheta.Reset();

    fEnergy = (MEnergyEst*)pl->FindObject("MEnergyEst");
    if (!fEnergy)
        *fLog << warn << "MEnergyEst not found... ignored." << endl;

    fPointPos = (MPointingPos*)pl->FindObject("MPointingPos");
    if (!fPointPos)
        *fLog << warn << "MPointingPos not found... ignored." << endl;

    fTimeEffOn = (MTime*)pl->FindObject("MTimeEffectiveOnTime", "MTime");
    if (!fTimeEffOn)
        *fLog << warn << "MTimeEffectiveOnTime [MTime] not found... ignored." << endl;

    fTime = (MTime*)pl->FindObject("MTime");
    if (!fTime)
        *fLog << warn << "MTime not found... ignored." << endl;

    fResult = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "Significance");
    if (!fResult)
        return kFALSE;

    //fExcess = (MParameterD*)const_cast<MParList*>(pl)->FindCreateObj("MParameterD", "MExcess");
    //if (!fExcess)
    //    return kFALSE;

    fLastTime = MTime();

    MBinning binst, binse, binsa;
    binst.SetEdges(fHAlpha, 'x');
    binse.SetEdges(fHAlpha, 'y');
    binsa.SetEdges(fHAlpha, 'z');

    MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
    if (fPointPos && bins)
        binst.SetEdges(*bins);
    if (!fPointPos)
        binst.SetEdges(1, 0, 90);

    bins = (MBinning*)pl->FindObject("BinningEnergyEst", "MBinning");
    if (fEnergy && bins)
        binse.SetEdges(*bins);
    if (!fEnergy)
        binse.SetEdges(1, 10, 100000);

    bins = (MBinning*)pl->FindObject("BinningAlpha", "MBinning");
    if (bins)
        binsa.SetEdges(*bins);

    binse.Apply(fHEnergy);
    binst.Apply(fHTheta);
    binsa.Apply(fHAlphaTime);
    MH::SetBinning(&fHAlpha, &binst, &binse, &binsa);

    MAlphaFitter *fit = (MAlphaFitter*)pl->FindObject("MAlphaFitter");
    if (!fit)
        *fLog << warn << "MAlphaFitter not found... ignored." << endl;
    else
        fit->Copy(fFit);

    fFit.Print();

    return kTRUE;
}

void MHAlpha::InitAlphaTime(const MTime &t)
{
    //
    // If this is the first call we have to initialize the time-histogram
    //
    MBinning bins;
    bins.SetEdges(1, t.GetAxisTime()-60, t.GetAxisTime());
    bins.Apply(fHTime);

    fLastTime=t;
}

void MHAlpha::UpdateAlphaTime(Bool_t final)
{
    if (!fTimeEffOn)
        return;

    const Int_t steps = 10;

    static int rebin = steps;

    if (!final)
    {
        if (fLastTime==MTime() && *fTimeEffOn!=MTime())
        {
            InitAlphaTime(*fTimeEffOn);
            return;
        }

        if (fLastTime!=*fTimeEffOn)
        {
            fLastTime=*fTimeEffOn;
            rebin--;
        }

        if (rebin>0)
            return;
    }

    MAlphaFitter fit(fFit);
    if (!fit.Fit(fHAlphaTime))
        return;

    // Reset Histogram
    fHAlphaTime.Reset();

    // Get number of bins
    const Int_t n = fHTime.GetNbinsX();

    //
    // Prepare Histogram
    //

    MTime *time = final ? fTime : fTimeEffOn;
    if (final)
        time->Plus1ns();

    // Enhance binning
    MBinning bins;
    bins.SetEdges(fHTime, 'x');
    bins.AddEdge(time->GetAxisTime());
    bins.Apply(fHTime);

    //
    // Fill histogram
    //
    fHTime.SetBinContent(n+1, fit.GetEventsExcess());
    fHTime.SetBinError(n+1, fit.GetEventsExcess()*0.1);

    *fLog << all << *fTimeEffOn << ": " << fit.GetEventsExcess() << endl;

    rebin = steps;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Bool_t MHAlpha::Fill(const MParContainer *par, const Stat_t w)
{
    Double_t alpha, energy, theta;

    if (fMatrix)
    {
        alpha  = (*fMatrix)[fMap[0]];
        energy = 1000; 
        theta  =    0; 
    }
    else
    {
        const MHillasSrc *hil = dynamic_cast<const MHillasSrc*>(par);
        if (!par)
        {
            *fLog << err << dbginf << "MHillasSrc not found... abort." << endl;
            return kFALSE;
        }

        alpha  = hil->GetAlpha();
        energy = fEnergy   ? fEnergy->GetEnergy() : 1000;
        theta  = fPointPos ? fPointPos->GetZd()   : 0;
    }

    UpdateAlphaTime();

    fHAlpha.Fill(theta, energy, TMath::Abs(alpha), w);
    fHAlphaTime.Fill(TMath::Abs(alpha), w);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Paint the integral and the error on top of the histogram
//
void MHAlpha::PaintText(Double_t val, Double_t error) const
{
    TLatex text(0.45, 0.94, Form("N_{exc} = %.1f \\pm %.1f", val, error));
    text.SetBit(TLatex::kTextNDC);
    text.SetTextSize(0.04);
    text.Paint();
}

// --------------------------------------------------------------------------
//
// Update the projections
//
void MHAlpha::Paint(Option_t *opt)
{
    // Note: Any object cannot be added to a pad twice!
    //       The object is searched by gROOT->FindObject only in
    //       gPad itself!
    TVirtualPad *padsave = gPad;

    TH1D *h0=0;

    TString o(opt);
    if (o==(TString)"proj")
    {
        TPaveStats *st=0;
        for (int x=0; x<4; x++)
        {
            TVirtualPad *p=padsave->GetPad(x+1);
            if (!p || !(st = (TPaveStats*)p->GetPrimitive("stats")))
                continue;

            if (st->GetOptStat()==11)
                continue;

            const Double_t y1 = st->GetY1NDC();
            const Double_t y2 = st->GetY2NDC();
            const Double_t x1 = st->GetX1NDC();
            const Double_t x2 = st->GetX2NDC();

            st->SetY1NDC((y2-y1)/3+y1);
            st->SetX1NDC((x2-x1)/3+x1);
            st->SetOptStat(11);
        }

        padsave->cd(1);
        fHAlpha.ProjectionZ(fNameProjAlpha);
        FitEnergyBins();
        FitThetaBins();
    }

    if (o==(TString)"alpha")
        if ((h0 = (TH1D*)gPad->FindObject(fNameProjAlpha)))
        {
            // Do not store the 'final' result in fFit
            MAlphaFitter fit(fFit);
            fit.Fit(*h0, kTRUE);
            fit.PaintResult();
        }

    if (o==(TString)"time")
        PaintText(fHTime.Integral(), 0);

    if (o==(TString)"theta")
        PaintText(fHTheta.Integral(), 0);

    if (o==(TString)"energy")
    {
        if (fHEnergy.GetMaximum()>1)
        {
            gPad->SetLogx();
            gPad->SetLogy();
        }
        FitEnergySpec(kTRUE);

    }

    gPad = padsave;
}

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

    // Do the projection before painting the histograms into
    // the individual pads
    AppendPad("proj");

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

    TH1D *h=0;

    pad->cd(1);
    gPad->SetBorderMode(0);
    h = fHAlpha.ProjectionZ(fNameProjAlpha, -1, 9999, -1, 9999, "E");
    h->SetBit(TH1::kNoTitle);
    h->SetXTitle("\\alpha [\\circ]");
    h->SetYTitle("Counts");
    h->SetDirectory(NULL);
    h->SetMarkerStyle(kFullDotMedium);
    h->SetBit(kCanDelete);
    h->Draw();
    // After the Alpha-projection has been drawn. Fit the histogram
    // and paint the result into this pad
    AppendPad("alpha");

    if (fHEnergy.GetNbinsX()>1)
    {
        pad->cd(2);
        gPad->SetBorderMode(0);
        fHEnergy.Draw();
        AppendPad("energy");
    }
    else
        delete pad->GetPad(2);

    if (fTimeEffOn && fTime || fHTime.GetNbinsX()>1)
    {
        pad->cd(3);
        gPad->SetBorderMode(0);
        fHTime.Draw();
        AppendPad("time");
    }
    else
        delete pad->GetPad(3);

    if (fHTheta.GetNbinsX()>1)
    {
        pad->cd(4);
        gPad->SetBorderMode(0);
        fHTheta.Draw();
        AppendPad("theta");
    }
    else
        delete pad->GetPad(4);

}

Bool_t MHAlpha::Finalize()
{
    // Store the final result in fFit
    TH1D *h = fHAlpha.ProjectionZ("AlphaExc_px", -1, 9999, -1, 9999, "E");
    h->SetDirectory(0);
    fFit.Print();
    Bool_t rc = fFit.Fit(*h);
    delete h;
    if (!rc)
    {
        *fLog << warn << "Histogram empty." << endl;
        return kTRUE;
    }

    if (fResult)
        fResult->SetVal(fFit.GetSignificance());

    FitEnergyBins();
    FitThetaBins();
    UpdateAlphaTime(kTRUE);
    MH::RemoveFirstBin(fHTime);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of
// MMcEvt. This function adds all necessary columns to the
// given matrix. Afterward you should fill the matrix with the corresponding
// data (eg from a file by using MHMatrix::Fill). If you now loop
// through the matrix (eg using MMatrixLoop) MHHadronness::Fill
// will take the values from the matrix instead of the containers.
//
void MHAlpha::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
        return;

    fMatrix = mat;

    fMap[0] = fMatrix->AddColumn("MHillasSrc.fAlpha");
    //fMap[1] = fMatrix->AddColumn("MEnergyEst.fEnergy");
}

void MHAlpha::StopMapping()
{
    fMatrix = NULL; 
}