/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MAlphaFitter
//
// 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 "MAlphaFitter.h"

#include <TF1.h>
#include <TH1.h>
#include <TH3.h>

#include <TRandom.h>

#include <TLatex.h>

#include "MMath.h"

#include "MLogManip.h"

ClassImp(MAlphaFitter);

using namespace std;

void MAlphaFitter::Clear(Option_t *o)
{
    fSignificance=0;
    fEventsExcess=0;
    fEventsSignal=0;
    fEventsBackground=0;

    fChiSqSignal=0;
    fChiSqBg=0;
    fIntegralMax=0;
    fScaleFactor=1;

    fCoefficients.Reset();
}

// --------------------------------------------------------------------------
//
// This is a preliminary implementation of a alpha-fit procedure for
// all possible source positions. It will be moved into its own
// more powerfull class soon.
//
// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
// or
//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
//
// Parameter [1] is fixed to 0 while the alpha peak should be
// symmetric around alpha=0.
//
// Parameter [4] is fixed to 0 because the first derivative at
// alpha=0 should be 0, too.
//
// In a first step the background is fitted between bgmin and bgmax,
// while the parameters [0]=0 and [2]=1 are fixed.
//
// In a second step the signal region (alpha<sigmax) is fittet using
// the whole function with parameters [1], [3], [4] and [5] fixed.
//
// The number of excess and background events are calculated as
//   s = int(hist,    0, 1.25*sigint)
//   b = int(pol2(3), 0, 1.25*sigint)
//
// The Significance is calculated using the Significance() member
// function.
//
Bool_t MAlphaFitter::Fit(TH1D &h, Bool_t paint)
{
    Clear();
    if (h.GetEntries()==0)
        return kFALSE;

    Double_t sigmax=fSigMax;
    Double_t bgmin =fBgMin;
    Double_t bgmax =fBgMax;

    //*fLog << inf << "Fit: " << sigmax << " " << fSigInt << "  " << bgmin << " " << bgmax << endl;

    //TF1 fFunc("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90);

    //fFunc->SetParameters(fCoefficients.GetArray());

    fFunc->FixParameter(1, 0);
    if (fPolynomOrder!=1)
        fFunc->FixParameter(4, 0);
    fFunc->SetParLimits(2, 0, 90);
    fFunc->SetParLimits(3, -1, 1);

    const Double_t alpha0 = h.GetBinContent(1);
    const Double_t alphaw = h.GetXaxis()->GetBinWidth(1);

    // Check for the regios which is not filled...
    if (alpha0==0)
        return kFALSE; //*fLog << warn << "Histogram empty." << endl;

    // First fit a polynom in the off region
    fFunc->FixParameter(0, 0);
    fFunc->FixParameter(2, 1);
    fFunc->ReleaseParameter(3);

    for (int i=5; i<fFunc->GetNpar(); i++)
        if (fFitBackground)
            fFunc->ReleaseParameter(i);
        else
            fFunc->SetParameter(i, 0);

    // options : N  do not store the function, do not draw
    //           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
    //           E  Perform better Errors estimation using Minos technique
    if (fFitBackground)
    {
        h.Fit(fFunc, "NQI", "", bgmin, bgmax);
        fChiSqBg = fFunc->GetChisquare()/fFunc->GetNDF();
    }


    // ------------------------------------
    if (paint && fFitBackground)
    {
        fFunc->SetRange(0, 90);
        fFunc->SetLineColor(kRed);
        fFunc->SetLineWidth(2);
        fFunc->Paint("same");
    }
    // ------------------------------------

    fFunc->ReleaseParameter(0);
    //func.ReleaseParameter(1);
    fFunc->ReleaseParameter(2);
    fFunc->FixParameter(3, fFunc->GetParameter(3));
    fFunc->FixParameter(4, fFunc->GetParameter(4));
    for (int i=5; i<fFunc->GetNpar(); i++)
        fFunc->FixParameter(i, fFunc->GetParameter(i));

    // Do not allow signals smaller than the background
    const Double_t A  = alpha0-fFunc->GetParameter(3);
    const Double_t dA = TMath::Abs(A);
    fFunc->SetParLimits(0, -dA*4, dA*4);
    fFunc->SetParLimits(2, 0, 90);

    // Now fit a gaus in the on region on top of the polynom
    fFunc->SetParameter(0, A);
    fFunc->SetParameter(2, sigmax*0.75);

    // options : N  do not store the function, do not draw
    //           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
    //           E  Perform better Errors estimation using Minos technique
    h.Fit(fFunc, "NQI", "", 0, sigmax);

    fChiSqSignal = fFunc->GetChisquare()/fFunc->GetNDF();
    fCoefficients.Set(fFunc->GetNpar(), fFunc->GetParameters());

    //const Bool_t ok = NDF>0 && chi2<2.5*NDF;

    // ------------------------------------
    if (paint)
    {
        fFunc->SetLineColor(kGreen);
        fFunc->SetLineWidth(2);
        fFunc->Paint("same");
    }
    // ------------------------------------

     //const Double_t s = fFunc->Integral(0, fSigInt)/alphaw;
     fFunc->SetParameter(0, 0);
     fFunc->SetParameter(2, 1);
     //const Double_t b = fFunc->Integral(0, fSigInt)/alphaw;
     //fSignificance = MMath::SignificanceLiMaSigned(s, b);

    const Int_t bin = h.GetXaxis()->FindFixBin(fSigInt*0.999);

    fIntegralMax      = h.GetBinLowEdge(bin+1);
    fEventsBackground = fFunc->Integral(0, fIntegralMax)/alphaw;
    fEventsSignal     = h.Integral(1, bin);
    fEventsExcess     = fEventsSignal-fEventsBackground;
    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground);

    if (TMath::IsNaN(fSignificance))
        fSignificance=0;

    if (fEventsExcess<0)
        fEventsExcess=0;

    return kTRUE;
}

Bool_t MAlphaFitter::Fit(const TH1D &hon, const TH1D &hof, Double_t alpha, Bool_t paint)
{
    TH1D h(hon);
    h.Add(&hof, -1); // substracts also number of entries!
    h.SetEntries(hon.GetEntries());

    MAlphaFitter fit(*this);
    fit.SetPolynomOrder(0);

    if (alpha<=0 || !fit.Fit(h, paint))
        return kFALSE;

    fChiSqSignal  = fit.GetChiSqSignal();
    fChiSqBg      = fit.GetChiSqBg();
    fCoefficients = fit.GetCoefficients();

    const Int_t bin = hon.GetXaxis()->FindFixBin(fSigInt*0.999);

    fIntegralMax      = hon.GetBinLowEdge(bin+1);
    fEventsBackground = hof.Integral(1, bin);
    fEventsSignal     = hon.Integral(1, bin);
    fEventsExcess     = fEventsSignal-fEventsBackground;
    fScaleFactor      = alpha;
    fSignificance     = MMath::SignificanceLiMaSigned(fEventsSignal, fEventsBackground/alpha, alpha);

    if (TMath::IsNaN(fSignificance))
        fSignificance=0;
    if (fEventsExcess<0)
        fEventsExcess=0;

    return kTRUE;
}

void MAlphaFitter::PaintResult(Float_t x, Float_t y, Float_t size) const
{
    const Double_t w  = GetGausSigma();
    const Double_t m  = fIntegralMax;

    const Int_t    l1 = w<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(w));
    const Int_t    l2 = m<=0 ? 0 : (Int_t)TMath::Ceil(-TMath::Log10(m));
    const TString fmt = Form("\\sigma_{L/M}=%%.1f  \\omega=%%.%df\\circ  E=%%d B=%%d  x<%%.%df  \\tilde\\chi_{b}=%%.1f  \\tilde\\chi_{s}=%%.1f  c=%%.1f  f=%%.2f",
                             l1<1?1:l1+1, l2<1?1:l2+1);

    TLatex text(x, y, Form(fmt.Data(), fSignificance, w, (int)fEventsExcess,
                           (int)fEventsBackground, m, fChiSqBg, fChiSqSignal,
                           fCoefficients[3], fScaleFactor));

    text.SetBit(TLatex::kTextNDC);
    text.SetTextSize(size);
    text.Paint();
}

void MAlphaFitter::Copy(TObject &o) const
{
    MAlphaFitter &f = static_cast<MAlphaFitter&>(o);

    // Setup
    f.fSigInt       = fSigInt;
    f.fSigMax       = fSigMax;
    f.fBgMin        = fBgMin;
    f.fBgMax        = fBgMax;
    f.fScaleMin     = fScaleMin;
    f.fScaleMax     = fScaleMax;
    f.fPolynomOrder = fPolynomOrder;
    f.fFitBackground= fFitBackground;
    f.fSignalFunc   = fSignalFunc;
    f.fScaleMode    = fScaleMode;
    f.fScaleUser    = fScaleUser;
    f.fStrategy     = fStrategy;
    f.fCoefficients.Set(fCoefficients.GetSize());
    f.fCoefficients.Reset();

    // Result
    f.fSignificance     = fSignificance;
    f.fEventsExcess     = fEventsExcess;
    f.fEventsSignal     = fEventsSignal;
    f.fEventsBackground = fEventsBackground;
    f.fChiSqSignal      = fChiSqSignal;
    f.fChiSqBg          = fChiSqBg;
    f.fIntegralMax      = fIntegralMax;
    f.fScaleFactor      = fScaleFactor;

    // Function
    TF1 *fcn = f.fFunc;
    f.fFunc = new TF1(*fFunc);
    gROOT->GetListOfFunctions()->Remove(f.fFunc);
    f.fFunc->SetName("Dummy");
    delete fcn;
}

void MAlphaFitter::Print(Option_t *o) const
{
    *fLog << GetDescriptor() << ": Fitting..." << endl;
    *fLog << " ...signal to " << fSigMax << " (integrate into bin at " << fSigInt << ")" << endl;
    *fLog << " ...signal function: ";
    switch (fSignalFunc)
    {
    case kGauss:       *fLog << "gauss(x)";        break;
    case kThetaSq:     *fLog << "gauss(sqrt(x))";  break;
    }
    *fLog << endl;
    if (!fFitBackground)
        *fLog << " ...no background." << endl;
    else
    {
        *fLog << " ...background from " << fBgMin << " to " << fBgMax << endl;
        *fLog << " ...polynom order " << fPolynomOrder << endl;
        *fLog << " ...scale mode: ";
        switch (fScaleMode)
        {
        case kNone:        *fLog << "none.";         break;
        case kEntries:     *fLog << "entries.";      break;
        case kIntegral:    *fLog << "integral.";     break;
        case kOffRegion:   *fLog << "off region (integral between " << fScaleMin << " and " << fScaleMax << ")"; break;
        case kBackground:  *fLog << "background (integral between " << fBgMin    << " and " << fBgMax    << ")"; break;
        case kLeastSquare: *fLog << "least square (N/A)"; break;
        case kUserScale:   *fLog << "user def (" << fScaleUser << ")"; break;
        }
        *fLog << endl;
    }

    if (TString(o).Contains("result"))
    {
        *fLog << "Result:" << endl;
        *fLog << " - Significance (Li/Ma)    " << fSignificance << endl;
        *fLog << " - Excess Events           " << fEventsExcess << endl;
        *fLog << " - Signal Events           " << fEventsSignal << endl;
        *fLog << " - Background Events       " << fEventsBackground << endl;
        *fLog << " - Chi^2/ndf (Signal)      " << fChiSqSignal << endl;
        *fLog << " - Chi^2/ndf (Background)  " << fChiSqBg << endl;
        *fLog << " - Signal integrated up to " << fIntegralMax << "°" << endl;
        *fLog << " - Scale Factor (Off)      " << fScaleFactor << endl;
    }
}

Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, UInt_t bin, Bool_t paint)
{
    const TString name(Form("TempAlphaEnergy%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, -1, 9999, bin, bin, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}

Bool_t MAlphaFitter::FitTheta(const TH3D &hon, UInt_t bin, Bool_t paint)
{
    const TString name(Form("TempAlphaTheta%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, bin, bin, -1, 9999, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}
/*
Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint)
{
    const TString name(Form("TempAlphaTime%06d", gRandom->Integer(1000000)));

    hon.GetZaxis()->SetRange(bin,bin);
    TH1D *h = (TH1D*)hon.Project3D("ye");
    hon.GetZaxis()->SetRange(-1,9999);

    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}
*/
Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint)
{
    const TString name(Form("TempAlpha%06d", gRandom->Integer(1000000)));
    TH1D *h = hon.ProjectionZ(name, -1, 9999, -1, 9999, "E");
    h->SetDirectory(0);

    const Bool_t rc = Fit(*h, paint);
    delete h;
    return rc;
}

Bool_t MAlphaFitter::FitEnergy(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));

    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    const Bool_t rc = ScaleAndFit(*h1, h0, paint);

    delete h0;
    delete h1;

    return rc;
}

Bool_t MAlphaFitter::FitTheta(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));

    TH1D *h1 = hon.ProjectionZ(name1, bin, bin, -1, 9999, "E");
    TH1D *h0 = hof.ProjectionZ(name0, bin, bin, -1, 9999, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    const Bool_t rc = ScaleAndFit(*h1, h0, paint);

    delete h0;
    delete h1;

    return rc;
}
/*
Bool_t MAlphaFitter::FitTime(const TH3D &hon, const TH3D &hof, UInt_t bin, Bool_t paint)
{
    const TString name1(Form("TempAlphaTime%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlphaTime%06d_off", gRandom->Integer(1000000)));

    hon.GetZaxis()->SetRange(bin,bin);
    TH1D *h1 = (TH1D*)hon.Project3D("ye");
    hon.GetZaxis()->SetRange(-1,9999);
    h1->SetDirectory(0);

    hof.GetZaxis()->SetRange(bin,bin);
    TH1D *h0 = (TH1D*)hof.Project3D("ye");
    hof.GetZaxis()->SetRange(-1,9999);
    h0->SetDirectory(0);

    const Bool_t rc = ScaleAndFit(*h1, h0, paint);

    delete h0;
    delete h1;

    return rc;
}
*/
Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, const TH3D &hof, Bool_t paint)
{
    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));

    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, -1, 9999, "E");
    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, -1, 9999, "E");
    h1->SetDirectory(0);
    h0->SetDirectory(0);

    const Bool_t rc = ScaleAndFit(*h1, h0, paint);

    delete h0;
    delete h1;

    return rc;
}

Double_t MAlphaFitter::Scale(TH1D &of, const TH1D &on) const
{
    Float_t scaleon = 1;
    Float_t scaleof = 1;
    switch (fScaleMode)
    {
    case kNone:
        return 1;

    case kEntries:
        scaleon = on.GetEntries();
        scaleof = of.GetEntries();
        break;

    case kIntegral:
        scaleon = on.Integral();
        scaleof = of.Integral();
        break;

    case kOffRegion:
        {
            const Int_t min = on.GetXaxis()->FindFixBin(fScaleMin);
            const Int_t max = on.GetXaxis()->FindFixBin(fScaleMax);
            scaleon = on.Integral(min, max);
            scaleof = of.Integral(min, max);
        }
        break;

    case kBackground:
        {
            const Int_t min = on.GetXaxis()->FindFixBin(fBgMin);
            const Int_t max = on.GetXaxis()->FindFixBin(fBgMax);
            scaleon = on.Integral(min, max);
            scaleof = of.Integral(min, max);
        }
        break;

    case kUserScale:
        scaleon = fScaleUser;
        break;

        // This is just to make some compiler happy
    default:
        return 1;
    }

    if (scaleof!=0)
    {
        of.Scale(scaleon/scaleof);
        return scaleon/scaleof;
    }
    else
    {
        of.Reset();
        return 0;
    }
}

Double_t MAlphaFitter::GetMinimizationValue() const
{
    switch (fStrategy)
    {
    case kSignificance:
        return -GetSignificance();
    case kSignificanceChi2:
        return -GetSignificance()/GetChiSqSignal();
    case kSignificanceLogExcess:
        if (GetEventsExcess()<1)
            return 0;
        return -GetSignificance()*TMath::Log10(GetEventsExcess());
    case kSignificanceExcess:
        return -GetSignificance()*GetEventsExcess();
    case kExcess:
        return -GetEventsExcess();
    case kGaussSigma:
        return GetGausSigma();
    }
    return 0;
}

Int_t MAlphaFitter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;

    //void SetScaleUser(Float_t scale)       { fScaleUser = scale; fScaleMode=kUserScale; }
    //void SetScaleMode(ScaleMode_t mode)    { fScaleMode    = mode; }

    if (IsEnvDefined(env, prefix, "SignalIntegralMax", print))
    {
        SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalIntegralMax", fSigInt));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "SignalFitMax", print))
    {
        SetSignalIntegralMax(GetEnvValue(env, prefix, "SignalFitMax", fSigMax));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "BackgroundFitMax", print))
    {
        SetBackgroundFitMax(GetEnvValue(env, prefix, "BackgroundFitMax", fBgMax));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "BackgroundFitMin", print))
    {
        SetBackgroundFitMin(GetEnvValue(env, prefix, "BackgroundFitMin", fBgMin));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "ScaleMin", print))
    {
        SetScaleMin(GetEnvValue(env, prefix, "ScaleMin", fScaleMin));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "ScaleMax", print))
    {
        SetScaleMax(GetEnvValue(env, prefix, "ScaleMax", fScaleMax));
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "PolynomOrder", print))
    {
        SetPolynomOrder(GetEnvValue(env, prefix, "PolynomOrder", fPolynomOrder));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "MinimizationStrategy", print))
    {
        TString txt = GetEnvValue(env, prefix, "MinimizationStrategy", "");
        txt = txt.Strip(TString::kBoth);
        txt.ToLower();
        if (txt==(TString)"significance")
            fStrategy = kSignificance;
        if (txt==(TString)"significancechi2")
            fStrategy = kSignificanceChi2;
        if (txt==(TString)"significanceexcess")
            fStrategy = kSignificanceExcess;
        if (txt==(TString)"excess")
            fStrategy = kExcess;
        if (txt==(TString)"gausssigma" || txt==(TString)"gaussigma")
            fStrategy = kGaussSigma;
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "Scale", print))
    {
        fScaleUser = GetEnvValue(env, prefix, "Scale", fScaleUser);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "ScaleMode", print))
    {
        TString txt = GetEnvValue(env, prefix, "ScaleMode", "");
        txt = txt.Strip(TString::kBoth);
        txt.ToLower();
        if (txt==(TString)"none")
            fScaleMode = kNone;
        if (txt==(TString)"entries")
            fScaleMode = kEntries;
        if (txt==(TString)"integral")
            fScaleMode = kIntegral;
        if (txt==(TString)"offregion")
            fScaleMode = kOffRegion;
        if (txt==(TString)"background")
            fScaleMode = kBackground;
        if (txt==(TString)"leastsquare")
            fScaleMode = kLeastSquare;
        if (txt==(TString)"userscale")
            fScaleMode = kUserScale;
        if (txt==(TString)"fixed")
        {
            fScaleMode = kUserScale;
            fScaleUser = fScaleFactor;
            cout << "---------> " << fScaleFactor << " <----------" << endl;
        }
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "SignalFunction", print))
    {
        TString txt = GetEnvValue(env, prefix, "SignalFunction", "");
        txt = txt.Strip(TString::kBoth);
        txt.ToLower();
        if (txt==(TString)"gauss" || txt==(TString)"gaus")
            SetSignalFunction(kGauss);
        if (txt==(TString)"thetasq")
            SetSignalFunction(kThetaSq);
        rc = kTRUE;
    }

    return rc;
}