#ifndef MARS_MHCalibrationChargeBlindPix
#define MARS_MHCalibrationChargeBlindPix


#ifndef MARS_MHCalibrationChargePix
#include "MHCalibrationChargePix.h"
#endif

#ifndef ROOT_TMatrix
#include <TMatrix.h>
#endif

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

class TH1F;
class TF1;
class TPaveText;
class TText;
class MRawEvtData;
class MRawEvtPixelIter;
class MCalibrationChargeBlindPix;
class MExtractBlindPixel;
class MExtractedSignalBlindPixel;
class MHCalibrationChargeBlindPix : public MHGausEvents
{
private:

  static const Int_t    fgChargeNbins;       //! Default for fNBins        (now set to: 5300   )
  static const Axis_t   fgChargeFirst;       //! Default for fFirst        (now set to: -100.5 )
  static const Axis_t   fgChargeLast;        //! Default for fLast         (now set to: 5199.5 )
  static const Float_t  fgSinglePheCut;      //! Default for fSinglePheCut (now set to: 200    )
  static const Float_t  fgNumSinglePheLimit; //! Default for fNumSinglePheLimit (now set to: 50)
  static const Float_t  gkSignalInitializer; //! Signal initializer (-9999.)
  
  static const Double_t gkElectronicAmp;     // Electronic Amplification after the PMT (in FADC counts/N_e)
  static const Double_t gkElectronicAmpErr;  // Error of the electronic amplification

  Float_t fSinglePheCut;                     // Value of summed FADC slices upon which event considered as single-phe
  Float_t fNumSinglePheLimit;                // Minimum number of single-phe events 

  MCalibrationChargeBlindPix *fBlindPix;     //! Storage container results  
  MExtractedSignalBlindPixel *fSignal;       //! Storage container extracted signal
  MRawEvtData                *fRawEvt;       //! Storage container raw data
 
  TVector fASinglePheFADCSlices;             // Averaged FADC slice entries supposed single-phe events
  TVector fAPedestalFADCSlices;              // Averaged FADC slice entries supposed pedestal   events
 
  TF1 *fSinglePheFit;                        // Single Phe Fit (Gaussians convoluted with Poisson) 

  UInt_t  fNumSinglePhes;                    // Number of entries in fASinglePheFADCSlices 
  UInt_t  fNumPedestals;                     // Number of entries in fAPedestalFADCSlices  

  Double_t  fLambda;                         // Poisson mean from Single-phe fit 
  Double_t  fLambdaCheck;                    // Poisson mean from Pedestal fit alone
  Double_t  fMu0;                            // Mean of the pedestal
  Double_t  fMu1;                            // Mean of single-phe peak
  Double_t  fSigma0;                         // Sigma of the pedestal
  Double_t  fSigma1;                         // Sigma of single-phe peak
  Double_t  fLambdaErr;                      // Error of Poisson mean from Single-phe fit 
  Double_t  fLambdaCheckErr;                 // Error of Poisson mean from Pedestal fit alone 
  Double_t  fMu0Err;                         // Error of  Mean of the pedestal    
  Double_t  fMu1Err;                         // Error of  Mean of single-phe peak 
  Double_t  fSigma0Err;                      // Error of  Sigma of the pedestal   
  Double_t  fSigma1Err;                      // Error of  Sigma of single-phe peak
  Double_t  fChisquare;                      // Chisquare of single-phe fit 
  Int_t     fNDF;                            // Ndof of single-phe fit 
  Double_t  fProb;                           // Probability of singleo-phe fit
  Double_t  fMeanPedestal;                   // Mean pedestal from pedestal run
  Double_t  fSigmaPedestal;                  // Sigma pedestal from pedestal run
  Double_t  fMeanPedestalErr;                // Error of Mean pedestal from pedestal run 
  Double_t  fSigmaPedestalErr;               // Error of Sigma pedestal from pedestal run
				     
  Byte_t    fFlags;                          // Bit-field for the flags
  enum { kSinglePheFitOK, kPedestalFitOK };  // Possible bits to be set

  TPaveText *fFitLegend;                     //! Some legend to display the fit results
  TH1F      *fHSinglePheFADCSlices;          //! A histogram created and deleted only in Draw()
  TH1F      *fHPedestalFADCSlices;           //! A histogram created and deleted only in Draw()

  // Fill histos
  void  FillSinglePheFADCSlices(const MRawEvtPixelIter &iter);
  void  FillPedestalFADCSlices( const MRawEvtPixelIter &iter);

  // Fit
  Bool_t InitFit();
  void   ExitFit();  
  
public:

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

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

  Bool_t SetupFill(const MParList *pList);
  Bool_t ReInit   (      MParList *pList);
  Bool_t Fill     (const MParContainer *par, const Stat_t w=1);
  Bool_t Finalize();
  
  // Getters
  const Double_t GetLambda        ()  const { return fLambda;         }
  const Double_t GetLambdaCheck   ()  const { return fLambdaCheck;    }
  const Double_t GetMu0           ()  const { return fMu0;            }
  const Double_t GetMu1           ()  const { return fMu1;            }
  const Double_t GetSigma0        ()  const { return fSigma0;         }
  const Double_t GetSigma1        ()  const { return fSigma1;         }
  const Double_t GetLambdaErr     ()  const { return fLambdaErr;      }
  const Double_t GetLambdaCheckErr()  const { return fLambdaCheckErr; }
  const Double_t GetMu0Err        ()  const { return fMu0Err;         }
  const Double_t GetMu1Err        ()  const { return fMu1Err;         }
  const Double_t GetSigma0Err     ()  const { return fSigma0Err;      }
  const Double_t GetSigma1Err     ()  const { return fSigma1Err;      }
  const Float_t  GetSinglePheCut  ()  const { return fSinglePheCut;   }
 
  TVector &GetASinglePheFADCSlices()             { return fASinglePheFADCSlices;  }
  const TVector &GetASinglePheFADCSlices() const { return fASinglePheFADCSlices;  }

  TVector &GetAPedestalFADCSlices()              { return fAPedestalFADCSlices;  }  
  const TVector &GetAPedestalFADCSlices()  const { return fAPedestalFADCSlices;  }  

  const Bool_t  IsSinglePheFitOK()         const;
  const Bool_t  IsPedestalFitOK()          const;
  
  // Setters
  void SetCalibrationChargeBlindPix ( MCalibrationChargeBlindPix *pix)    { fBlindPix          = pix;  }
  void SetSinglePheCut      ( const Float_t cut =fgSinglePheCut      )    { fSinglePheCut      = cut;  }
  void SetNumSinglePheLimit ( const Float_t lim =fgNumSinglePheLimit )    { fNumSinglePheLimit = lim;  }

  void SetMeanPedestal      ( const Float_t f )   { fMeanPedestal     = f;  }
  void SetMeanPedestalErr   ( const Float_t f )   { fMeanPedestalErr  = f;  }
  void SetSigmaPedestal     ( const Float_t f )   { fSigmaPedestal    = f;  }
  void SetSigmaPedestalErr  ( const Float_t f )   { fSigmaPedestalErr = f;  }

  void SetSinglePheFitOK    ( const Bool_t b=kTRUE);
  void SetPedestalFitOK     ( const Bool_t b=kTRUE);
  
  // Draws
  void Draw(Option_t *opt="");

private:
  void DrawLegend(Option_t *opt="");
  
  // Fits
public:
  enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; // The possible fit functions

private:
  FitFunc_t fFitFunc;

public:
  Bool_t FitSinglePhe (Option_t *opt="RL0+Q");
  void   FitPedestal  (Option_t *opt="RL0+Q");

  void   ChangeFitFunc(const FitFunc_t func)  { fFitFunc = func;  }
  
  // Simulation
  Bool_t SimulateSinglePhe(const Double_t lambda,
                           const Double_t mu0,    const Double_t mu1,
                           const Double_t sigma0, const Double_t sigma1);
  
private:

  inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
    {
      
      Double_t lambda1cat = par[0];  
      Double_t lambda1dyn = par[1];
      Double_t mu0        = par[2];
      Double_t mu1cat     = par[3];
      Double_t mu1dyn     = par[4];
      Double_t sigma0     = par[5];
      Double_t sigma1cat  = par[6];
      Double_t sigma1dyn  = par[7];
      
      Double_t sumcat = 0.;
      Double_t sumdyn = 0.;
      Double_t arg    = 0.;
      
      if (lambda1cat < lambda1dyn)
        return FLT_MAX;

      if (mu1cat    < mu0)
        return FLT_MAX;

      if (mu1dyn    < mu0)
        return FLT_MAX;

      if (mu1cat < mu1dyn)
        return FLT_MAX;

      if (sigma0 < 0.0001)
        return FLT_MAX;
      
      if (sigma1cat < sigma0)
        return FLT_MAX;

      if (sigma1dyn < sigma0)
        return FLT_MAX;

      Double_t mu2cat = (2.*mu1cat)-mu0;  
      Double_t mu2dyn = (2.*mu1dyn)-mu0;  
      Double_t mu3cat = (3.*mu1cat)-(2.*mu0);
      Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0);
      
      Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0));  
      Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0));  
      Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0));
      Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0));
      
      Double_t lambda2cat = lambda1cat*lambda1cat;
      Double_t lambda2dyn = lambda1dyn*lambda1dyn;
      Double_t lambda3cat = lambda2cat*lambda1cat;
      Double_t lambda3dyn = lambda2dyn*lambda1dyn;

     // k=0:
      arg = (x[0] - mu0)/sigma0;
      sumcat = TMath::Exp(-0.5*arg*arg)/sigma0;
      sumdyn = sumcat;

      // k=1cat:
      arg = (x[0] - mu1cat)/sigma1cat;
      sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat;
      // k=1dyn:
      arg = (x[0] - mu1dyn)/sigma1dyn;
      sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn;
      
      // k=2cat:
      arg = (x[0] - mu2cat)/sigma2cat;
      sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat;
      // k=2dyn:
      arg = (x[0] - mu2dyn)/sigma2dyn;
      sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn;
  
     
      // k=3cat:
      arg = (x[0] - mu3cat)/sigma3cat;
      sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat;
      // k=3dyn:
      arg = (x[0] - mu3dyn)/sigma3dyn;
      sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn;  
    
      sumcat = TMath::Exp(-1.*lambda1cat)*sumcat;
      sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn;
      
      return par[8]*(sumcat+sumdyn)/2.;

    }
    
  inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
    {

      Double_t lambda = par[0];  
      
      Double_t sum = 0.;
      Double_t arg = 0.;
      
      Double_t mu0 = par[1];
      Double_t mu1 = par[2];
      
      if (mu1 < mu0)
        return FLT_MAX;

      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];

      if (sigma0 < 0.0001)
        return FLT_MAX;
      
      if (sigma1 < sigma0)
        return FLT_MAX;
      
      Double_t mu2 = (2.*mu1)-mu0;  
      Double_t mu3 = (3.*mu1)-(2.*mu0);
      Double_t mu4 = (4.*mu1)-(3.*mu0);
      
      Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
      Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
      Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
      
      Double_t lambda2 = lambda*lambda;
      Double_t lambda3 = lambda2*lambda;
      Double_t lambda4 = lambda3*lambda;
      
      // k=0:
      arg = (x[0] - mu0)/sigma0;
      sum = TMath::Exp(-0.5*arg*arg)/sigma0;
      
      // k=1:
      arg = (x[0] - mu1)/sigma1;
      sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
      
      // k=2:
      arg = (x[0] - mu2)/sigma2;
      sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
      
      // k=3:
      arg = (x[0] - mu3)/sigma3;
      sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
      
      // k=4:
      arg = (x[0] - mu4)/sigma4;
      sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
      
      return TMath::Exp(-1.*lambda)*par[5]*sum;
      
    } 

  
  inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
    {
      
      Double_t lambda = par[0];  
      
      Double_t sum = 0.;
      Double_t arg = 0.;
      
      Double_t mu0 = par[1];
      Double_t mu1 = par[2];
      
      if (mu1 < mu0)
        return FLT_MAX;
      
      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma0 < 0.0001)
        return FLT_MAX;
      
      if (sigma1 < sigma0)
        return FLT_MAX;
      
      
      Double_t mu2 = (2.*mu1)-mu0;  
      Double_t mu3 = (3.*mu1)-(2.*mu0);
      Double_t mu4 = (4.*mu1)-(3.*mu0);
      Double_t mu5 = (5.*mu1)-(4.*mu0);
      
      Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
      Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
      Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
      Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
      
      Double_t lambda2 = lambda*lambda;
      Double_t lambda3 = lambda2*lambda;
      Double_t lambda4 = lambda3*lambda;
      Double_t lambda5 = lambda4*lambda;
      
      // k=0:
      arg = (x[0] - mu0)/sigma0;
      sum = TMath::Exp(-0.5*arg*arg)/sigma0;
      
      // k=1:
      arg = (x[0] - mu1)/sigma1;
      sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
      
      // k=2:
      arg = (x[0] - mu2)/sigma2;
      sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
      
      // k=3:
      arg = (x[0] - mu3)/sigma3;
      sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
      
      // k=4:
      arg = (x[0] - mu4)/sigma4;
      sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
      
      // k=5:
      arg = (x[0] - mu5)/sigma5;
      sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
      
      return TMath::Exp(-1.*lambda)*par[5]*sum;
      
    }
  
  
  inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
    {
      
      Double_t lambda = par[0];  
      
      Double_t sum = 0.;
      Double_t arg = 0.;
      
      Double_t mu0 = par[1];
      Double_t mu1 = par[2];
      
      if (mu1 < mu0)
        return FLT_MAX;
      
      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma0 < 0.0001)
        return FLT_MAX;
      
      if (sigma1 < sigma0)
        return FLT_MAX;
      
      
      Double_t mu2 = (2.*mu1)-mu0;  
      Double_t mu3 = (3.*mu1)-(2.*mu0);
      Double_t mu4 = (4.*mu1)-(3.*mu0);
      Double_t mu5 = (5.*mu1)-(4.*mu0);
      Double_t mu6 = (6.*mu1)-(5.*mu0);
      
      Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
      Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
      Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
      Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
      Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
      
      Double_t lambda2 = lambda*lambda;
      Double_t lambda3 = lambda2*lambda;
      Double_t lambda4 = lambda3*lambda;
      Double_t lambda5 = lambda4*lambda;
      Double_t lambda6 = lambda5*lambda;
      
      // k=0:
      arg = (x[0] - mu0)/sigma0;
      sum = TMath::Exp(-0.5*arg*arg)/sigma0;
      
      // k=1:
      arg = (x[0] - mu1)/sigma1;
      sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
      
      // k=2:
      arg = (x[0] - mu2)/sigma2;
      sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
      
      // k=3:
      arg = (x[0] - mu3)/sigma3;
      sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
      
      // k=4:
      arg = (x[0] - mu4)/sigma4;
      sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
      
      // k=5:
      arg = (x[0] - mu5)/sigma5;
      sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
      
      // k=6:
      arg = (x[0] - mu6)/sigma6;
      sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
      
      return TMath::Exp(-1.*lambda)*par[5]*sum;
      
    }

  inline static Double_t fPolya(Double_t *x, Double_t *par)
    {

      const Double_t QEcat = 0.247;            // mean quantum efficiency
      const Double_t sqrt2 = 1.4142135623731;
      const Double_t sqrt3 = 1.7320508075689;
      const Double_t sqrt4 = 2.;
      
      const Double_t lambda = par[0];           // mean number of photons
      
      const Double_t excessPoisson = par[1];    // non-Poissonic noise contribution
      const Double_t delta1 = par[2];           // amplification first dynode
      const Double_t delta2 = par[3];           // amplification subsequent dynodes
      
      const Double_t electronicAmpl = par[4];   // electronic amplification and conversion to FADC charges

      const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2;  // total PMT gain
      const Double_t A = 1. + excessPoisson - QEcat                        
        + 1./delta1 
                + 1./delta1/delta2
        + 1./delta1/delta2/delta2;                                  // variance contributions from PMT and QE
      
      const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl;        // Total gain and conversion
      
      const Double_t mu0 = par[7];                                      // pedestal
      const Double_t mu1 = totAmpl;                                 // single phe position
      const Double_t mu2 = 2*totAmpl;                               // double phe position
      const Double_t mu3 = 3*totAmpl;                               // triple phe position
      const Double_t mu4 = 4*totAmpl;                               // quadruple phe position
      
      const Double_t sigma0 = par[5];
      const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A);
      const Double_t sigma2 = sqrt2*sigma1;
      const Double_t sigma3 = sqrt3*sigma1;
      const Double_t sigma4 = sqrt4*sigma1;
      
      const Double_t lambda2 = lambda*lambda;
      const Double_t lambda3 = lambda2*lambda;
      const Double_t lambda4 = lambda3*lambda;
      
      //-- calculate the area----
      Double_t arg = (x[0] - mu0)/sigma0;
      Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0;
      
     // k=1:
      arg = (x[0] - mu1)/sigma1;
      sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
      
      // k=2:
      arg = (x[0] - mu2)/sigma2;
      sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
      
      // k=3:
      arg = (x[0] - mu3)/sigma3;
      sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
      
      // k=4:
      arg = (x[0] - mu4)/sigma4;
      sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;      
      
      return TMath::Exp(-1.*lambda)*par[6]*sum;
    }
  
  ClassDef(MHCalibrationChargeBlindPix, 1)  // Histogram class for Charge Blind Pixel Calibration
};

#endif  /* MARS_MHCalibrationChargeBlindPix */