#ifndef MARS_MAlphaFitter
#define MARS_MAlphaFitter

#ifndef MARS_MParContainer
#include "MParContainer.h"
#endif

#ifndef ROOT_TArrayD
#include <TArrayD.h>
#endif

#ifndef ROOT_TF1
#include <TF1.h>
#endif

class TH1D;
class TH3D;

class MAlphaFitter : public MParContainer
{
public:
    enum ScaleMode_t {
        kNone,
        kEntries,
        kIntegral,
        kOffRegion,
        kLeastSquare,
        kUserScale
    };
private:
    Float_t fSigInt;
    Float_t fSigMax;
    Float_t fBgMin;
    Float_t fBgMax;
    Float_t fScaleMin;
    Float_t fScaleMax;
    Int_t   fPolynomOrder;

    Double_t fSignificance;
    Double_t fEventsExcess;
    Double_t fEventsSignal;
    Double_t fEventsBackground; // fScaleFactor already applied

    Double_t fChiSqSignal;
    Double_t fChiSqBg;
    Double_t fIntegralMax;
    Double_t fScaleFactor;  // Scale factor for off-data

    TArrayD fCoefficients;

    TF1 *fFunc;

    ScaleMode_t fScaleMode;
    Double_t fScaleUser;

public:
    // Implementing the function yourself is only about 5% faster
    MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kOffRegion), fScaleUser(1)
    {
        fName  = name  ? name  : "MAlphaFitter";
        fTitle = title ? title : "Fit alpha";

        gROOT->GetListOfFunctions()->Remove(fFunc);
        fFunc->SetName("Dummy");

        Clear();
    }

    MAlphaFitter(const MAlphaFitter &f) : fFunc(0)
    {
        f.Copy(*this);
    }
    ~MAlphaFitter()
    {
        delete fFunc;
    }

    void Clear(Option_t *o="");
    void Print(Option_t *o="") const;
    void Copy(TObject &o) const;
    /*
    TObject *Clone(const char *name) const
    {
        return new MAlphaFitter(*this);
    }
    */

    void SetScaleUser(Float_t scale)       { fScaleUser = scale; fScaleMode=kUserScale; }
    void SetScaleMode(ScaleMode_t mode)    { fScaleMode    = mode; }
    void SetSignalIntegralMax(Float_t s)   { fSigInt       = s; }
    void SetSignalFitMax(Float_t s)        { fSigMax       = s; }
    void SetBackgroundFitMin(Float_t s)    { fBgMin        = s; }
    void SetBackgroundFitMax(Float_t s)    { fBgMax        = s; }
    void SetScaleMin(Float_t s)            { fScaleMin     = s; }
    void SetScaleMax(Float_t s)            { fScaleMax     = s; }
    void SetPolynomOrder(Int_t s)          { fPolynomOrder = s; delete fFunc; fFunc=new TF1 ("", Form("gaus(0) + pol%d(3)", s));
        gROOT->GetListOfFunctions()->Remove(fFunc);
        fFunc->SetName("Dummy");
        fCoefficients.Set(3+s+1); fCoefficients.Reset(); }

    Double_t GetEventsExcess() const       { return fEventsExcess; }
    Double_t GetEventsSignal() const       { return fEventsSignal; }
    Double_t GetEventsBackground() const   { return fEventsBackground; }

    Double_t GetSignificance() const       { return fSignificance; }
    Double_t GetChiSqSignal() const        { return fChiSqSignal; }
    Double_t GetChiSqBg() const            { return fChiSqBg; }
    Double_t GetScaleFactor() const        { return fScaleFactor; }

    Double_t GetGausSigma() const          { return fCoefficients[2]; }
    Double_t GetGausMu() const             { return fCoefficients[1]; }
    Double_t GetGausA() const              { return fCoefficients[0]; }
    Double_t GetCoefficient(Int_t i) const { return fCoefficients[i]; }
    const TArrayD &GetCoefficients() const { return fCoefficients; }

    void PaintResult(Float_t x=0.04, Float_t y=0.94, Float_t size=0.035) const;

    Bool_t Fit(TH1D &h, Bool_t paint=kFALSE);
    Bool_t Fit(TH1D &on, TH1D &off, Double_t alpha, Bool_t paint=kFALSE);
    Bool_t Fit(TH1D &on, TH1D *off, Double_t alpha, Bool_t paint=kFALSE)
    {
        return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
    }
    Bool_t Fit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
    {
        return off ? Fit(on, *off, 1, paint) : Fit(on, paint);
    }
    Bool_t ScaleAndFit(TH1D &on, TH1D *off, Bool_t paint=kFALSE)
    {
        const Double_t alpha = off ? Scale(*off, on) : 1;
        return off ? Fit(on, *off, alpha, paint) : Fit(on, paint);
    }

    Bool_t FitAlpha(const TH3D &h, Bool_t paint=kFALSE);
    Bool_t FitEnergy(const TH3D &h, UInt_t bin, Bool_t paint=kFALSE);
    Bool_t FitTheta(const TH3D &h,  UInt_t bin, Bool_t paint=kFALSE);

    Bool_t FitAlpha(const TH3D &on, const TH3D &off, Bool_t paint=kFALSE);
    Bool_t FitEnergy(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);
    Bool_t FitTheta(const TH3D &on, const TH3D &off, UInt_t bin, Bool_t paint=kFALSE);

    Bool_t FitAlpha(const TH3D &on, const TH3D *off, Bool_t paint=kFALSE)
    {
        return off ? FitAlpha(on, *off, paint) : FitAlpha(on, paint);
    }
    Bool_t FitEnergy(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
    {
        return off ? FitEnergy(on, *off, bin, paint) : FitEnergy(on, bin, paint);
    }
    Bool_t FitTheta(const TH3D &on, const TH3D *off, UInt_t bin, Bool_t paint=kFALSE)
    {
        return off ? FitTheta(on, *off, bin, paint) : FitTheta(on, bin, paint);
    }

    Double_t Scale(TH1D &off, const TH1D &on) const;

    ClassDef(MAlphaFitter, 1)
};

#endif