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

//////////////////////////////////////////////////////////////////////////////
//
// MHDisp
//
// Create a false source plot using disp.
//
// Currently the use of this class requires to be after MFMagicCuts
// in the tasklist. Switching of the M3L cut in MFMagicCuts is recommended.
//
//////////////////////////////////////////////////////////////////////////////
#include "MHDisp.h"

#include <TStyle.h>
#include <TCanvas.h>

#include <TF1.h>
#include <TF2.h>
#include <TProfile.h>

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

#include "MMath.h"

#include "MParList.h"
#include "MParameters.h"

#include "MHillasExt.h"
#include "MHillasSrc.h"
#include "MSrcPosCam.h"
#include "MPointingPos.h"

ClassImp(MHDisp);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHDisp::MHDisp(const char *name, const char *title)
    : fHilExt(0), fDisp(0)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHDisp";
    fTitle = title ? title : "3D-plot using Disp vs x, y";

    fHist.SetDirectory(NULL);

    fHist.SetName("Alpha");
    fHist.SetTitle("3D-plot of ThetaSq vs x, y");
    fHist.SetXTitle("x [\\circ]");
    fHist.SetYTitle("y [\\circ]");
    fHist.SetZTitle("\\vartheta [deg^{2}]");
}

// --------------------------------------------------------------------------
//
// Set binnings (takes BinningFalseSource) and prepare filling of the
// histogram.
//
// Also search for MTime, MObservatory and MPointingPos
//
Bool_t MHDisp::SetupFill(const MParList *plist)
{
    if (!MHFalseSource::SetupFill(plist))
        return kFALSE;

    fDisp = (MParameterD*)plist->FindObject("Disp", "MParameterD");
    if (!fDisp)
    {
        *fLog << err << "Disp [MParameterD] not found... abort." << endl;
        return kFALSE;
    }

    MParameterD *p = (MParameterD*)plist->FindObject("M3lCut", "MParameterD");
    if (!p)
    {
        *fLog << err << "M3lCut [MParameterD] not found... abort." << endl;
        return kFALSE;
    }
    fM3lCut = TMath::Abs(p->GetVal());

    p = (MParameterD*)plist->FindObject("DispXi", "MParameterD");
    if (!p)
    {
        *fLog << err << "DispXi [MParameterD] not found... abort." << endl;
        return kFALSE;
    }
    fXi = p->GetVal();

    p = (MParameterD*)plist->FindObject("DispXiTheta", "MParameterD");
    if (!p)
    {
        *fLog << err << "DispXiTheta [MParameterD] not found... abort." << endl;
        return kFALSE;
    }
    fXiTheta = p->GetVal();

    fHilExt = (MHillasExt*)plist->FindObject("MHillasExt");
    if (!fHilExt)
    {
        *fLog << err << "MHillasExt not found... aborting." << endl;
        return kFALSE;
    }

    // Initialize all bins with a small (=0) value otherwise
    // these bins are not displayed
    for (int x=0; x<fHist.GetNbinsX(); x++)
        for (int y=0; y<fHist.GetNbinsX(); y++)
            fHist.Fill(fHist.GetXaxis()->GetBinCenter(x+1),
                       fHist.GetYaxis()->GetBinCenter(y+1),
                       0.0, 1e-10);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
{
    const MHillas *hil = dynamic_cast<const MHillas*>(par);
    if (!hil)
    {
        *fLog << err << "MHDisp::Fill: No container specified!" << endl;
        return kFALSE;
    }

    // Get camera rotation angle
    Double_t rho = 0;
    if (fTime && fObservatory && fPointPos)
        rho = fPointPos->RotationAngle(*fObservatory, *fTime);

    // Get Disp from Parlist
    const Double_t disp = fDisp->GetVal();

    // Calculate where disp is pointing
    TVector2 pos1(hil->GetCosDelta(), hil->GetSinDelta());
    TVector2 pos2(hil->GetCosDelta(), hil->GetSinDelta());
    pos1 *=  disp;
    pos2 *= -disp;

    pos1 += hil->GetMean()*fMm2Deg;
    pos2 += hil->GetMean()*fMm2Deg;

    // gammaweight: If we couldn't decide which position makes the
    // event a gamma, both position are assigned 'half of a gamma'
    Double_t gweight = 0.5;

    // Check whether our g/h-separation allows to asign definitly
    // to one unique position. Therefor we requeire that the event
    // would be considered a gamma for one, but not for the other
    // position. This can only be the case if the third moment
    // has a value higher than the absolute cut value.
    if (TMath::Abs(fHilExt->GetM3Long()) > fM3lCut)
    {
        // Because at one position the event is considered a gamma
        // we have to find out which position it is...
        MSrcPosCam src;
        MHillasSrc hsrc;
        hsrc.SetSrcPos(&src);

        // Calculate the sign for one of the desired source positions
        // The other position must have the opposite sign
        src.SetXY(pos1/fMm2Deg);
        if (hsrc.Calc(*hil)>0)
        {
            *fLog << warn << "Calculation of MHillasSrc failed." << endl;
            return kFALSE;
        }
        const Double_t m3l = fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, hsrc.GetCosDeltaAlpha());

        gweight = m3l>fM3lCut ? 1 : 0;
    }

    // Now we can safly derotate both position...
    TVector2 srcp;
    if (fSrcPos)
        srcp = fSrcPos->GetXY();

    if (rho!=0)
    {
        pos1=pos1.Rotate(-rho);
        pos2=pos2.Rotate(-rho);
        srcp=srcp.Rotate(-rho);
    }

    pos1 -= srcp*fMm2Deg;
    pos2 -= srcp*fMm2Deg;

    fHist.Fill(pos1.X(), pos1.Y(), 0.0, w*gweight);
    fHist.Fill(pos2.X(), pos2.Y(), 0.0, w*(1-gweight));

    return kTRUE;
}
/*
static Double_t FcnGauss2d(Double_t *x, Double_t *par)
{
    TVector2 v = TVector2(x[0], x[1]).Rotate(par[5]*TMath::DegToRad());

    const Double_t g0 = TMath::Gaus(v.X(), par[1], par[2]);
    const Double_t g1 = TMath::Gaus(v.Y(), par[3], par[4]);

    //Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE);
    return par[0]*g0*g1  + par[6]*(v.X()*v.X() + v.Y()*v.Y()) +  par[7];
}*/

void MHDisp::Paint(Option_t *o)
{
    TVirtualPad *pad = gPad;

    pad->cd(1);

    fHist.GetZaxis()->SetRange(0,9999);
    TH1 *h1=fHist.Project3D("yx_on");
    gStyle->SetPalette(1, 0);

    Int_t ix, iy, iz;
    TF1 func("fcn", "gaus + [3]*x*x + [4]");

    pad->cd(3);
    TH1 *h2 = (TH1*)gPad->FindObject("RadProf");
    /*
    if (!h2)
    {
        const Double_t maxr = TMath::Hypot(h1->GetXaxis()->GetXmax(), h1->GetYaxis()->GetXmax());
        const Int_t    nbin = (h1->GetNbinsX()+h1->GetNbinsY())/2;
        TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
        //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
        h->Sumw2();
        h->SetXTitle("\\vartheta [\\circ]");
        h->SetYTitle("<cts>/\\Delta R");
        h->SetBit(kCanDelete);
        h->Draw();
    }*/
    if (h1 && h2)
    {
        h2->Reset();

        h1->GetMaximumBin(ix, iy, iz);

        const Double_t x0 = h1->GetXaxis()->GetBinCenter(ix);
        const Double_t y0 = h1->GetYaxis()->GetBinCenter(iy);
        const Double_t w0 = h1->GetXaxis()->GetBinWidth(1);

        for (int x=0; x<h1->GetNbinsX(); x++)
            for (int y=0; y<h1->GetNbinsY(); y++)
            {
                const Double_t r = TMath::Hypot(h1->GetXaxis()->GetBinCenter(x+1)-x0,
                                                h1->GetYaxis()->GetBinCenter(y+1)-y0);
                h2->Fill(r, h1->GetBinContent(x+1, y+1));
            }

        // Replace with MAlphaFitter?
        func.SetLineWidth(1);
        func.SetLineColor(kBlue);

        func.SetParLimits(2, h2->GetBinWidth(1), 1.0);

        func.SetParameter(0, h2->GetBinContent(1));
        func.FixParameter(1, 0);
        func.SetParameter(2, 0.15);
        func.SetParameter(4, h2->GetBinContent(10));
        h2->Fit(&func, "IMQ", "", 0, 1.0);

        // No wintegrate the function f(x) per Delta Area
        // which is f(x)/(pi*delta r*(2*r+delta r))
        TF1 func2("fcn2", Form("(gaus + [3]*x*x + [4])/(2*x+%.5f)", w0));
        for (int i=0; i<5; i++)
            func2.SetParameter(i, func.GetParameter(i));

        const Double_t r0 = 2*func2.GetParameter(2);
        const Double_t e  = func2.Integral(0, r0)/(w0*TMath::Pi());
        func2.SetParameter(0, 0);
        const Double_t b  = func2.Integral(0, r0)/(w0*TMath::Pi());
        const Double_t s  = MMath::SignificanceLiMa(e, b);

        h2->SetTitle(Form("P=(%.2f\\circ/%.2f\\circ) \\omega=%.2f\\circ \\sigma=%.1f E=%.0f B=%.0f", x0, y0, func.GetParameter(2), s, e-b, b));
    }
   /*
    if (h1)
    {
        const Double_t maxr = 0.9*TMath::Abs(fHist.GetXaxis()->GetXmax());

        TF2 f2d("Gaus2D", FcnGauss2d, -maxr, maxr, -maxr, maxr, 8);
        f2d.SetLineWidth(1);

        f2d.SetParameter(0, h1->GetMaximum()*5); // A
        f2d.SetParLimits(1, h1->GetXaxis()->GetBinCenter(ix)-h1->GetXaxis()->GetBinWidth(ix)*5, h1->GetXaxis()->GetBinCenter(ix)+h1->GetXaxis()->GetBinWidth(ix));  // mu_1
        f2d.SetParLimits(3, h1->GetYaxis()->GetBinCenter(iy)-h1->GetYaxis()->GetBinWidth(iy)*5, h1->GetYaxis()->GetBinCenter(iy)+h1->GetYaxis()->GetBinWidth(iy));  // mu_2
        f2d.SetParLimits(2, 0, func.GetParameter(2)*5);          // sigma_1
        f2d.SetParLimits(4, 0, func.GetParameter(2)*5);          // sigma_2
        f2d.SetParLimits(5, 0, 45);            // phi
        f2d.SetParLimits(6, 0, func.GetParameter(3)*5);
        f2d.SetParLimits(7, 0, func.GetParameter(4)*5);

        f2d.SetParameter(0, h1->GetMaximum()); // A
        f2d.SetParameter(1, h1->GetXaxis()->GetBinCenter(ix)); // mu_1
        f2d.SetParameter(2, func.GetParameter(2));             // sigma_1
        f2d.SetParameter(3, h1->GetYaxis()->GetBinCenter(iy)); // mu_2
        f2d.SetParameter(4, func.GetParameter(2));             // sigma_2
        f2d.FixParameter(5, 0);                                // phi
        f2d.SetParameter(6, func.GetParameter(3));
        f2d.SetParameter(7, func.GetParameter(4));
        h1->Fit(&f2d, "Q", "cont2");
        //f2d.DrawCopy("cont2same");
    }*/
}

void MHDisp::Draw(Option_t *o)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    const Int_t col = pad->GetFillColor();
    pad->SetBorderMode(0);

    AppendPad("");

    TString name = Form("%s_1", pad->GetName());
    TPad *p = new TPad(name,name, 0.005, 0.005, 0.65, 0.995, col, 0, 0);
    p->SetNumber(1);
    p->Draw();
    p->cd();

    TH1 *h3 = fHist.Project3D("yx_on");
    h3->SetTitle("Distribution of equivalent events");
    h3->SetDirectory(NULL);
    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
    h3->SetMinimum(0);
    h3->Draw("colz");
    h3->SetBit(kCanDelete);
//    catalog->Draw("mirror same *");

    pad->cd();
    name = Form("%s_2", pad->GetName());
    p = new TPad(name,name, 0.66, 0.005, 0.995, 0.5, col, 0, 0);
    p->SetNumber(2);
    p->Draw();
    p->cd();
    h3->Draw("surf3");

    pad->cd();
    name = Form("%s_3", pad->GetName());
    p = new TPad(name,name, 0.66, 0.5, 0.995, 0.995, col, 0, 0);
    p->SetNumber(3);
    p->Draw();
    p->cd();

    const Double_t maxr = TMath::Hypot(h3->GetXaxis()->GetXmax(), h3->GetYaxis()->GetXmax());
    const Int_t    nbin = (h3->GetNbinsX()+h3->GetNbinsY())/2;
    TProfile *h = new TProfile("RadProf", "Radial Profile", nbin, 0, maxr);
    h->SetDirectory(0);
    //TH1F *h = new TH1F("RadProf", "Radial Profile", nbin, 0, maxr);
    //h->Sumw2();
    h->SetXTitle("\\vartheta [\\circ]");
    h->SetYTitle("<cts>/\\Delta R");
    h->SetBit(kCanDelete);
    h->Draw();
}