#ifndef MARS_MHGausEvents
#define MARS_MHGausEvents

#ifndef ROOT_TH1
#include <TH1.h>
#endif

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

class TVirtualPad;
class TGraph;
class TArrayF;
class TH1F;
class TH1I;
class TF1;
class MHGausEvents : public MH
{
private:

  const static Float_t  fgProbLimit;            // Default probability limit for judgement if fit is OK
  const static Int_t    fgNDFLimit;             // Default NDF limit for judgement if fit is OK
  const static Int_t    fgPowerProbabilityBins; // Default number of bins for the projected power spectrum 
  const static Int_t    fgBinsAfterStripping;   // Default number of bins for the Gauss Histogram after stripping off the zeros at both end

  Float_t  fProbLimit;                 // Probability limit for judgement if fit is OK 
  Int_t    fNDFLimit;                  // NDF limit for judgement if fit is OK
  Int_t    fPowerProbabilityBins;      // number of bins for the projected power spectrum
  Int_t    fBinsAfterStripping;        // number of bins for the Gauss Histogram after stripping off the zeros at both end
  Float_t  fEventFrequency;            // The event frequency in Hertz (to be set)
  
  TF1     *fFGausFit;                  //-> Gauss fit for fHGausHist
  TF1     *fFExpFit;                   //-> Exponential fit for FHPowerProbability
  
  TH1I    *fHPowerProbability;         //-> Fourier transform of fEvents projected on y-axis
  
  TArrayF *fPowerSpectrum;             //-> Fourier transform of fEvents

  TGraph  *fGraphEvents;               //!  TGraph to display the event array (will not be cloned!!)
  TGraph  *fGraphPowerSpectrum;        //!  TGraph to display the power spectrum array (will not be cloned!!)

  Float_t *CreateXaxis(Int_t n);       //   Create an x-axis for the TGraphs

  Double_t fMean;                      //   Mean of the Gauss fit
  Double_t fSigma;                     //   Sigma of the Gauss fit
  Double_t fMeanErr;                   //   Error of the mean of the Gauss fit
  Double_t fSigmaErr;                  //   Error of the sigma of the Gauss fit
  Double_t fProb;                      //   Probability of the Gauss fit (derived from Chi-Square and NDF

  enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK }; // Bits to hold information about fit results 
  
  Byte_t fFlags;                       //   Byte to hold the bits fit result bits
  
  UInt_t fCurrentSize;                 //   Current size of the array fEvents
  
protected:

  TH1F    fHGausHist;                  // Histogram which should hold the Gaussian distribution
  TArrayF fEvents;                     // Array which holds the entries of GausHist

  // Setters
  void  SetPowerProbabilityBins(const Int_t nbins=fgPowerProbabilityBins) {  fPowerProbabilityBins = nbins;  }
  void  SetBinsAfterStripping(const Int_t nbins=fgBinsAfterStripping)     { fBinsAfterStripping = nbins;  }

  void DrawEvents();                                  // Draw a graph of the array fEvents
  void DrawPowerSpectrum(TVirtualPad &pad, Int_t i);  // Draw a graph of the array fPowerSpectrum and the hist fHPowerProbability
  
public:

  MHGausEvents(const char* name=NULL, const char* title=NULL);
  ~MHGausEvents();

  virtual void Clear(Option_t *o="");
  virtual void Reset();  

  // Setters
  void  SetEventFrequency(const Float_t f=0)  {  fEventFrequency = f; }

  void  SetMean(     const Double_t d )   { fMean     = d;   }
  void  SetMeanErr(  const Double_t d )   { fMeanErr  = d;   }
  void  SetSigma(    const Double_t d )   { fSigma    = d;   }
  void  SetSigmaErr( const Double_t d )   { fSigmaErr = d;   }
  
  void  SetProbLimit( const Float_t lim=fgProbLimit ) {  fProbLimit = lim; }
  void  SetNDFLimit(  const Int_t   lim=fgNDFLimit  ) {  fNDFLimit = lim;  }  

  // Setters ONLY for MC:
  void  SetGausFitOK(         const Bool_t b );
  void  SetExpFitOK(          const Bool_t b );
  void  SetFourierSpectrumOK( const Bool_t b );

  // Getters
  const Double_t GetMean()          const { return fMean;      }
  const Double_t GetMeanErr()       const { return fMeanErr;   }
  const Double_t GetSigma()         const { return fSigma;     }
  const Double_t GetSigmaErr()      const { return fSigmaErr;  }
  const Double_t GetChiSquare()     const;
  const Double_t GetProb()          const { return fProb;      }
  const Int_t    GetNdf()           const;

  const Double_t GetSlope()        const;
  const Double_t GetOffset()       const;
  const Double_t GetExpChiSquare() const;
  const Double_t GetExpProb()      const;
  const Int_t    GetExpNdf()       const;

  TH1F *GetHGausHist()                     { return &fHGausHist;     }
  const TH1F *GetHGausHist()         const { return &fHGausHist;     } 

  TArrayF *GetEvents()                     { return &fEvents;        }  
  const TArrayF *GetEvents()         const { return &fEvents;        }
  
  TArrayF *GetPowerSpectrum()              { return fPowerSpectrum;  }  
  const TArrayF *GetPowerSpectrum()  const { return fPowerSpectrum;  }
  
  TF1 *GetFGausFit()                       { return fFGausFit;       }
  const TF1 *GetFGausFit()           const { return fFGausFit;       } 

  TH1I *GetHPowerProbability()             { return fHPowerProbability; }
  const TH1I *GetHPowerProbability() const { return fHPowerProbability; } 

  TF1 *GetFExpFit()                        { return fFExpFit;        }
  const TF1 *GetFExpFit()            const { return fFExpFit;        } 

  TGraph *GetGraphEvents()                 { return fGraphEvents;    }
  const TGraph *GetGraphEvents()     const { return fGraphEvents;    }
  
  TGraph *GetGraphPowerSpectrum()             { return fGraphPowerSpectrum;  }
  const TGraph *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum;  }
  
  const Bool_t IsGausFitOK()          const; 
  const Bool_t IsExpFitOK()           const; 
  const Bool_t IsEmpty()              const;
  const Bool_t IsFourierSpectrumOK()  const;

  // Fill
  void FillArray(const Float_t f);            // Fill only the array fEvents 
  Bool_t FillHist(const Float_t f);           // Fill only the histogram HGausHist 
  Bool_t FillHistAndArray(const Float_t f);   // Fill bothe the array fEvents and the histogram HGausHist
  
  // Fits
  Bool_t FitGaus(Option_t *option="RQ0");     // Fit the histogram HGausHist with a Gaussian

  // Draws
  virtual void Draw(Option_t *option="");     // Default Draw 
  
  // Prints
  virtual void Print(const Option_t *o="") const; // Default Print
  
  // Miscelleaneous
  void CreateFourierSpectrum();             // Create the fourier spectrum out of fEvents
  void CreateGraphEvents();                 // Create the TGraph fGraphEvents of fEvents 
  void CreateGraphPowerSpectrum();          // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum
  
  ClassDef(MHGausEvents, 1)                 // Base class for events with Gaussian distributed values
};

#endif