#ifndef MARS_MHGausEvents
#define MARS_MHGausEvents

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

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

#ifndef MARS_MArrayF
#include "MArrayF.h"
#endif

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

  const static Int_t    fgNDFLimit;             //! Default for fNDFLimit             (now set to: 2)
  const static Float_t  fgProbLimit;            //! Default for fProbLimit            (now set to: 0.001)
  const static Int_t    fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20)

  Float_t *CreateEventXaxis(Int_t n);  // Create an x-axis for the Event TGraphs
  Float_t *CreatePSDXaxis(Int_t n);    // Create an x-axis for the PSD TGraphs
  
protected:

  Float_t  fEventFrequency;              // Event frequency in Hertz (to be set)

  Int_t    fBinsAfterStripping;        // Bins for the Gauss Histogram after stripping off the zeros at both ends
  UInt_t   fCurrentSize;               // Current size of the array fEvents
  Byte_t   fFlags;                     // Bit field for the fit result bits
  Int_t    fPowerProbabilityBins;      // Bins for the projected power spectrum
  
  TH1I    *fHPowerProbability;         // Fourier transform of fEvents projected on y-axis
  MArrayF *fPowerSpectrum;             // Fourier transform of fEvents

  enum  { kGausFitOK,
          kExpFitOK,
          kFourierSpectrumOK,
          kExcluded };                 // Bits for information about fit results 
  
  MArrayF  fEvents;                    // Array which holds the entries of GausHist
  TF1     *fFGausFit;                  // Gauss fit for fHGausHist
  TF1     *fFExpFit;                   // Exponential fit for FHPowerProbability
  Axis_t   fFirst;                     // Lower histogram edge  for fHGausHist (used by InitBins()) 
  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!!)
  TH1F     fHGausHist;                 // Histogram to hold the Gaussian distribution
  Axis_t   fLast;                      // Upper histogram edge  for fHGausHist (used by InitBins()) 
  Double_t fMean;                      // Mean of the Gauss fit
  Double_t fMeanErr;                   // Error of the mean of the Gauss fit
  Int_t    fNbins;                     // Number histogram bins for fHGausHist (used by InitBins())
  Int_t    fNDFLimit;                  // NDF limit for judgement if fit is OK
  Double_t fSigma;                     // Sigma of the Gauss fit
  Double_t fSigmaErr;                  // Error of the sigma of the Gauss fit
  Double_t fProb;                      // Probability of the Gauss fit 
  Float_t  fProbLimit;                 // Probability limit for judgement if fit is OK 

  void DrawPowerSpectrum(TVirtualPad &pad, Int_t i);  // Draw graph of fPowerSpectrum and fHPowerProbability

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

public:

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

  TObject *Clone(const char* name="") const;
  
  void Clear(Option_t *o="");
  void Reset();  

  // Draws
  void Draw(Option_t *option="");       // Default Draw 
  void DrawEvents(Option_t *option="");   // Draw graph of fEvents

  // Inits
  virtual void InitBins();
  
  // Getters
  const Double_t GetChiSquare()          const;
  const Double_t GetExpChiSquare()       const;
  const Int_t    GetExpNdf()             const;
  const Double_t GetExpProb()            const;
        MArrayF *GetEvents()                   { return &fEvents;            }  
  const MArrayF *GetEvents()             const { return &fEvents;            }
  const Float_t  GetEventFrequency () const { return fEventFrequency; }
  TF1     *GetFExpFit()                  { return fFExpFit;            }
  const TF1     *GetFExpFit()            const { return fFExpFit;            } 
        TF1     *GetFGausFit()                 { return fFGausFit;           }
  const TF1     *GetFGausFit()           const { return fFGausFit;           } 
        TGraph  *GetGraphEvents()              { return fGraphEvents;        }
  const Double_t GetFirst()              const { return fFirst;              }
  const Double_t GetLast ()              const { return fLast ;              }  
  const TGraph  *GetGraphEvents()        const { return fGraphEvents;        }
        TGraph  *GetGraphPowerSpectrum()       { return fGraphPowerSpectrum; }
  const TGraph  *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum; }
        TH1F    *GetHGausHist()                { return &fHGausHist;         }
  const TH1F    *GetHGausHist()          const { return &fHGausHist;         } 
        TH1I    *GetHPowerProbability()        { return fHPowerProbability;  }
  const TH1I    *GetHPowerProbability()  const { return fHPowerProbability;  } 
  const Double_t GetMean()               const { return fMean;               }
  const Double_t GetMeanErr()            const { return fMeanErr;            }
  const Int_t    GetNdf()                const;
  const Int_t    GetNbins()               const { return fNbins;            }
  const Double_t GetOffset()             const;
        MArrayF *GetPowerSpectrum()            { return fPowerSpectrum;      }  
  const MArrayF *GetPowerSpectrum()      const { return fPowerSpectrum;      }
  const Double_t GetProb()               const { return fProb;               }
  const Double_t GetSigma()              const { return fSigma;              }
  const Double_t GetSigmaErr()           const { return fSigmaErr;           }
  const Double_t GetSlope()              const;

  const Bool_t IsExcluded()              const;
  const Bool_t IsExpFitOK()              const; 
  const Bool_t IsEmpty()                 const;
  const Bool_t IsFourierSpectrumOK()     const;
  const Bool_t IsGausFitOK()             const; 
  const Bool_t IsOnlyOverflow()          const;
  const Bool_t IsOnlyUnderflow()         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",
                   const Double_t xmin=0., 
	           const Double_t xmax=0.);       // Fit the histogram HGausHist with a Gaussian
  
  // Prints
  virtual void Print(const Option_t *o="") const; // Default Print
  
  // Setters
  void  SetEventFrequency   ( const Float_t  f                   ) { fEventFrequency = f;   }
  void  SetExcluded         ( const Bool_t   b=kTRUE             );  
  void  SetExpFitOK         ( const Bool_t   b=kTRUE             );
  void  SetFourierSpectrumOK( const Bool_t   b=kTRUE             );
  void  SetGausFitOK        ( const Bool_t   b=kTRUE             );
  void  SetLast             ( const Double_t d                   ) { fLast           = d;   }
  void  SetFirst            ( const Double_t d                   ) { fFirst          = d;   }
  void  SetMean             ( const Double_t d                   ) { fMean           = d;   }
  void  SetMeanErr          ( const Double_t d                   ) { fMeanErr        = d;   }
  void  SetNbins            ( const Int_t    i                   ) { fNbins          = i;   }  
  void  SetNDFLimit         ( const Int_t    lim=fgNDFLimit      ) { fNDFLimit       = lim; }  
  void  SetProb             ( const Double_t d                   ) { fProb           = d;   }
  void  SetProbLimit        ( const Float_t  lim=fgProbLimit     ) { fProbLimit      = lim; }
  void  SetSigma            ( const Double_t d                   ) { fSigma          = d;   }
  void  SetSigmaErr         ( const Double_t d                   ) { fSigmaErr       = d;   }

  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
