#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 Int_t    fgBinsAfterStripping;   //! Default for fBinsAfterStripping   (now set to: 40)
  const static Float_t  fgBlackoutLimit;        //! Default for fBlackoutLimit        (now set to: 5. )
  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 Float_t  fgPickupLimit;          //! Default for fPickupLimit          (now set to: 5. )
  const static Int_t    fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20)
  
  Int_t    fBinsAfterStripping;        // Bins for the Gauss Histogram after stripping off the zeros at both ends
  Int_t    fCurrentSize;               // Current size of the array fEvents
  Float_t  fEventFrequency;            // Event frequency in Hertz (to be set)
  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
  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!!)

  enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK, kExcluded }; // Bits for information about fit results 
  
protected:

  Float_t  fBlackoutLimit;             // Lower number sigmas from mean until event is considered blackout
  TArrayF  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()) 
  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
  Int_t    fSaturated;                 // Number of events classified as saturated
  Double_t fSigma;                     // Sigma of the Gauss fit
  Double_t fSigmaErr;                  // Error of the sigma of the Gauss fit
  Float_t  fPickupLimit;               // Upper number sigmas from mean until event is considered pickup
  Int_t    fPixId;                     // Pixel ID 
  Double_t fProb;                      // Probability of the Gauss fit 
  Float_t  fProbLimit;                 // Probability limit for judgement if fit is OK 

  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

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

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

 public:

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

  virtual void Clear(Option_t *o="");
  virtual void Reset();  
  virtual void InitBins();
  
  // Draws
  virtual void Draw(Option_t *option="");         // Default Draw 
  void DrawEvents();                              // Draw graph of fEvents
  
  // Getters
  const Double_t GetBlackout()           const;  
  const Double_t GetChiSquare()          const;
  const Double_t GetExpChiSquare()       const;
  const Int_t    GetExpNdf()             const;
  const Double_t GetExpProb()            const;
        TArrayF *GetEvents()                   { return &fEvents;            }  
  const TArrayF *GetEvents()             const { return &fEvents;            }
        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 Double_t GetOffset()             const;
  const Double_t GetPickup()             const;
  const Int_t    GetPixId()              const { return fPixId;              }
        TArrayF *GetPowerSpectrum()            { return fPowerSpectrum;      }  
  const TArrayF *GetPowerSpectrum()      const { return fPowerSpectrum;      }
  const Double_t GetProb()               const { return fProb;               }
  const Float_t  GetSaturated()          const { return fSaturated;          }
  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
  Bool_t RepeatFit(const Option_t *option="RQ0"); // Repeat fit within limits defined by fPickupLimit
  void   BypassFit();                             // Take mean and RMS from the histogram
  
  // Prints
  virtual void Print(const Option_t *o="") const; // Default Print
  
  // Setters
  void  SetBlackoutLimit    ( const Float_t  lim=fgBlackoutLimit ) { fBlackoutLimit  = lim; }
  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  SetPickupLimit      ( const Float_t  lim=fgPickupLimit   ) { fPickupLimit    = lim; }
  void  SetPixId            ( const Int_t    i                   ) { fPixId          = i;   }
  void  SetProb             ( const Double_t d                   ) { fProb           = d;   }
  void  SetProbLimit        ( const Float_t  lim=fgProbLimit     ) { fProbLimit      = lim; }
  void  SetSaturated        ( const Int_t    i                   ) { fSaturated     += i;   }
  void  SetSigma            ( const Double_t d                   ) { fSigma          = d;   }
  void  SetSigmaErr         ( const Double_t d                   ) { fSigmaErr       = d;   }

  // Miscelleaneous
  virtual void ChangeHistId(const Int_t id);      // Changes names and titles of the histogram
  virtual void Renorm();                          // Re-normalize the results 
  
  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
