#ifndef MARS_MHFindSignificanceONOFF
#define MARS_MHFindSignificanceONOFF

#ifndef MARS_MH
#include "MH.h"
#endif

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

#ifndef ROOT_TPostScript
#include <TPostScript.h>
#endif


class TF1;
class TH1;
class TCanvas;

class MHFindSignificanceONOFF : public MH
{
private:

    TH1  *fHistOrig;  // original plot of |alpha| (0.0 to 90.0 degrees)
    TH1  *fHist;      // copy of fHistOrig or rebinned histogram


    
   TH1  *fHistOrigOFF;  // original plot of |alpha| (0.0 to 90.0 degrees)
   TH1  *fHistOFF;      // copy of fHistOrig or rebinned histogram
   TH1  *fHistOFFNormalized;   //  fHistOFF normalized (contents and errors) with 
                               //   with fNormFactor; useful when plotting OFF 
                               //  Data together with ON data

 
  // Postscript file where, for each theta bin,  
   // alpha ON and OFF distributions 
  //  after cuts (and hence, Nex and SigmaLiMa computations) will be stored
   // This TPostScript object has to be initialized and closed finally 
   // outside this class.
 
   // NOT WORKING WELL !!!
   
   TPostScript* fPsFilename;
   

   TString fPsFilenameString; // Name of the file where plots will be stored while the 
  // TPostScript option is not working 



    Double_t fAlphamin;  // requested lower limit of fit range
    Double_t fAlphammm;  // center of fit range
    Double_t fAlphamax;  // requested lower limit of fit range

    Double_t fAlphaminOFF; // requested lower limit of fit range for OFF data
    Double_t fAlphamaxOFF;  // requested lower limit of fit range for OFF data


    Double_t fAlphami;  // actual lower limit of fit range 
    Double_t fAlphamm;  // actual center of fit range 
    Double_t fAlphama;  // actual upper limit of fit range

    Double_t fAlphasig; // requested signal range
    Double_t fAlphasi;  // actual signal range
    Double_t fAlphasiOFF;  // actual signal range for OFF data


    Double_t fAlphalow; // requested lower edge of signal range
    Double_t fAlphalo;  // actual lower edge of signal range

    Double_t fAlphahig; // requested upper edge of background range
    Double_t fAlphahi;  // actual upper edge of background range
 
    // number of events in signal region
    Double_t fNon;     // total number of events in signal region
    Double_t fNbg;     // number of background events in signal region
    Double_t fNex;     // number of excess events in signal region

    Double_t fdNon;
    Double_t fdNbg;
    Double_t fdNex;

    // number of events in background region
    Double_t fNbgtot;  // total number of events in background region
    Double_t fNbgtotFitted;  // fitted total no. of events in background region
    Double_t fdNbgtotFitted; // fitted error of this number

    // Variables related to OFF data. 
    // Many variables are defined for similarity with 
    // existing program. Many of them might be removed in future...


    Double_t fNexONOFF;     // number of excess events (ON-OFF) 
                            // in signal region
    Double_t fNexONOFFFitted;     // number of excess events (ON-OFF fitted) 
                                  // in signal region

    Double_t fdNexONOFF;     // error in number of excess events (ON-OFF) 
                            // in signal region
    Double_t fdNexONOFFFitted;     // error in number of excess events (ON-OFF fitted) 
                                  // in signal region


    
    Double_t fNoffTot;       // Total number of OFF events
    Double_t fNoffSig;    // Number of OFF events in "signal" region 
    Double_t fNoffBg;     // Number of OFF events in "background" region
    
    // errors for previous quantities

    Double_t fdNoffTot;       
    Double_t fdNoffSig;    
    Double_t fdNoffBg;     

    // Fitted quantities for OFF alpha distribution

    Double_t fNoffTotFitted;       // Total number of OFF events
    Double_t fNoffSigFitted;    // Number of OFF events in "signal" region 
    Double_t fNoffBgFitted;     // Number of OFF events in "background" region
    
    // errors for previous quantities

    Double_t fdNoffTotFitted;       
    Double_t fdNoffSigFitted;    
    Double_t fdNoffBgFitted;     

    
    // Boolean variable used to disable the usage ("serious" usage) of the 
    // quantities computed from fits. This will be useful in those cases 
    // where there is too few events to perform a decent fit to the 
    // alpha histograms.
    
    Bool_t fUseFittedQuantities;



    // Bool variable used to decide wether to print or not the results 
    // of the fit, significance, Nex... onto the final alpha plot. 
    // for the time being, this variable is set in the constructor. 
    // At some point, I might make it such it can be set externally...

    Bool_t fPrintResultsOntoAlphaPlot;


    // effective number of background events
    Double_t fNoff; // fNoff = fNoffSig * fNormFactor
    Double_t fGamma;   // 
    Double_t fNormFactor; // = TotalNon/TotalNoff)

    Double_t fSigLiMa; // significance of gamma signal according to Li & Ma
                       // Computed using quantities fGama and Noff
                       // and formula (17) 


    Double_t fSigLiMa2; // significance of gamma signal according to Li & Ma
                        // Computed using quantities  fNormFactor and fNoffSig
                        // and formula (17)

    Double_t fSigLiMa3; // significance of gamma signal according to Li & Ma
                        // Computed using quantities  fNormFactor and fNoffSig
                        // and formula (5)

                       


    const static Double_t fEps = 1.e-4;  // tolerance for floating point comparisons

    Bool_t fDraw;          // if true : draw plots
    Bool_t fSavePlots;     // if true : Save plots in Psfile
    Bool_t fFitGauss;      // if true : do the (polynomial+Gauss fit)
    Bool_t fRebin;         // if true : allow rebinning of the alpha plot    
    Bool_t fReduceDegree;  // if true : allow reducing of the order of the polynomial

    Bool_t fConstantBackg; // if set true if background fit is not possible
                           // due to low statistics

    TCanvas  *fCanvas;

    Double_t fNexGauss;    // total number of excess events 
                           // (from fitted Gauss function)
    Double_t fdNexGauss;   // error of the total number of excess events
 
    Double_t fSigmaGauss;  // sigma of fitted Gauss function
    Double_t fdSigmaGauss; // error of this sigma

    //--------------------
    TF1      *fPoly;   // polynomial function
    Int_t    fFitBad;  // if != 0 fit failed
    Int_t    fDegree;  // degree of polynomial to be fitted to the background
    Int_t    fNdf;     // number of degrees of freedom of polynomial fit
    Double_t fChisq;   // chi squared of polynomial fit
    Double_t fProb;    // chi squared probability ofg polynomial fit 

    TArrayD fValues;
    TArrayD fErrors;

    const static Int_t    fNdim = 6;
    Double_t fEmat[fNdim][fNdim];
    Double_t fEma [fNdim][fNdim];
    Double_t fCorr[fNdim][fNdim];

    Int_t  fMbins;     // number of bins in the fit range
    Int_t  fMlow;      // number of bins in the fit range with too few entries
    Int_t  fNzero;     // number of bins in the fit range with zero entry
    Int_t  fIstat;

    //--------------------

    // Fit to OFF data

    
    TF1      *fPolyOFF;   // polynomial function for OFF data
    
    TF1      *fPolyOFFNormalized;  //  polynomial function for OFF data normalized with fNormFactor*(BinWidthON/BinWidthOFF)

    Int_t    fFitBadOFF;  // if != 0 fit failed
    Int_t    fDegreeOFF;  // degree of polynomial to be fitted to the background
    Int_t    fNdfOFF;     // number of degrees of freedom of polynomial fit
    Double_t fChisqOFF;   // chi squared of polynomial fit
    Double_t fProbOFF;    // chi squared probability ofg polynomial fit 

    TArrayD fValuesOFF;
    TArrayD fErrorsOFF;
    TArrayD fValuesOFFNormalized;
    TArrayD fErrorsOFFNormalized;

    const static Int_t    fNdimOFF = 6;
    Double_t fEmatOFF[fNdim][fNdim];
    Double_t fEmaOFF [fNdim][fNdim];
    Double_t fCorrOFF[fNdim][fNdim];

    Int_t  fMbinsOFF;     // number of bins in the fit range
    Int_t  fMlowOFF;      // number of bins in the fit range with too few entries
    Int_t  fNzeroOFF;     // number of bins in the fit range with zero entry
    Int_t  fIstatOFF;


    //--------------------
    TF1      *fGPoly;   // (Gauss+polynomial) function
    TF1      *fGBackg;  // polynomial part of (Gauss+polynomial) function
    Int_t    fGFitBad;  // if != 0 fit failed
    Int_t    fGDegree;  // degree of polynomial to be fitted to the background
    Int_t    fGNdf;     // number of degrees of freedom of polynomial fit
    Double_t fGChisq;   // chi squared of polynomial fit
    Double_t fGProb;    // chi squared probability ofg polynomial fit 

    TArrayD fGValues;
    TArrayD fGErrors;

    const static Int_t    fGNdim = 9;
    Double_t fGEmat[fGNdim][fGNdim];
    Double_t fGEma[fGNdim][fGNdim];
    Double_t fGCorr[fGNdim][fGNdim];

    Int_t  fGMbins;     // number of bins in the fit range
    Int_t  fGNzero;     // numnber of bins in the fit range with zero entry
    Int_t  fGIstat;

    //--------------------

    static const TString gsDefName;  //! Default Name
    static const TString gsDefTitle; //! Default Title

    Bool_t DetExcess(); 
    Bool_t DetExcessONOFF(); 
    Bool_t FitPolynomial();
    Bool_t FitPolynomialOFF();
    Bool_t FitGaussPoly();
    //Bool_t FitGaussPolyONOFF();
    Bool_t RebinHistogram(Double_t x0, Int_t nrebin);
    Bool_t RebinHistogramOFF(Double_t x0, Int_t nrebin);

 public:
    MHFindSignificanceONOFF(const char *name=NULL, const char *title=NULL);
    ~MHFindSignificanceONOFF();

    
    // fhistOFF is the alpha distribution of OFF data (after cuts), 
    // BUT NOT normalized. Normalization factor is also one of the 
    // arguments. Off alpha distribution will be normalized 
    // at the beginning of the function FindSigmaONOFF. 

    Bool_t FindSigmaONOFF(TH1 *fhistON, TH1 *fhistOFF, 
			  Double_t NormFactor,
			  Double_t alphamin, Double_t alphamax,
			  Int_t degreeON,  Int_t degreeOFF, 
			  Double_t alphasig, 
			  Bool_t drawpoly, 
			  Bool_t fitgauss, 
			  Bool_t print, Bool_t saveplots, 
			  //TPostScript* PsFile
			  const TString psfilename);

    // Set limits for the polynomial fit of the OFF data
    // It will give the user the possibility of using, 
    // for the OFF data substraction, a fit in the region 
    // of interest. Might be useful in case of "pathological"  
    // alpha distribution for the OFF data.

    void SetLimitsForOFFDataFit(Double_t Low, Double_t Up)
	{ fAlphaminOFF = Low; fAlphamaxOFF = Up;  }

    void SetAlphaSig (Double_t alphasig)
	{ fAlphasig = alphasig; }


    void SetUseFittedQuantities (Bool_t b)
	{fUseFittedQuantities = b;}
	    

    Bool_t SigmaLiMa(Double_t non, Double_t noff, Double_t gamma,
                     Double_t *siglima);

    Bool_t SigmaLiMaForm5(Double_t non, Double_t noff, Double_t gamma,
			  Double_t *siglima);

   
    Bool_t SigmaVsAlphaONOFF(TH1 *fhistON, TH1 *fhistOFF,  
			     Double_t alphamin, Double_t alphamax, 
			     Int_t degree, Bool_t print);
    


    // Function that returns one of the 3 LiMa sigmas. 
    // The returned value is the one used in the optimization 
    // and final alpha plots. 
    

    Double_t GetSignificance();//   { return fSigLiMa; }

    Bool_t GetUseFittedQuantities() {return fUseFittedQuantities;}



    // Following function computes a clone of fHistOFF and normalizes 
    // contents, errors and fPolyOFF (if exists) with the fNormFactor. 
    // This normalized OFF hist will be used when plotting OFF data 
    // together with ON data.
 
    Bool_t ComputeHistOFFNormalized(); 


    Bool_t DrawFit(Option_t *opt=NULL);

    Bool_t DrawHistOFF();
    Bool_t DrawHistOFFNormalized();

    Float_t GetDegree()     const { return fDegree;  }
    Float_t GetProb()       const { return fProb;    }
    Float_t GetNdf()        const { return fNdf;     }
    Float_t GetGamma()      const { return fGamma;   }
    Float_t GetNon()        const { return fNon;     }
    Float_t GetNex()        const { return fNex;     }
    Float_t GetNbg()        const { return fNbg;     }
    Float_t GetSigLiMa1()    const { return fSigLiMa; }
    Float_t GetSigLiMa2()    const { return fSigLiMa2; }
    Float_t GetSigLiMa3()    const { return fSigLiMa3; }
       
    Float_t GetMbins()      const { return fMbins;   }
    Float_t GetAlphasi()    const { return fAlphasi; }
    Float_t GetNexONOFF() const {return fNexONOFF;}
    Float_t GetNexONOFFFitted() const {return fNexONOFFFitted;}
    
    

    void SetRebin(Bool_t b=kTRUE);
    void SetReduceDegree(Bool_t b=kTRUE);

    void SetPsFilename (TPostScript* PsFile) {fPsFilename = PsFile;}

    void PrintPoly(Option_t *opt=NULL);
    void PrintPolyOFF(Option_t *opt=NULL);
    void PrintPolyGauss(Option_t *opt=NULL);

    ClassDef(MHFindSignificanceONOFF, 1) // Determine significance from alpha plot
};

#endif




















