#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 MExtractedSignalBlindPixel;
class MHCalibrationChargeBlindPix : public MHCalibrationChargePix
{
private:

  MCalibrationChargeBlindPix *fBlindPix;  //! Storage container of the results  
  MExtractedSignalBlindPixel *fSignal;    //! Storage container of the extracted signal
  MRawEvtData                *fRawEvt;    //! Storage container of the raw data

  static const Int_t    fgChargeNbins;
  static const Axis_t   fgChargeFirst;
  static const Axis_t   fgChargeLast;  

  static const Float_t  fgSinglePheCut;  

  static const Float_t  fgNumSinglePheLimit;

  static const Double_t gkElectronicAmp;
  static const Double_t gkElectronicAmpErr;

  TVector fASinglePheFADCSlices;           // Array containing the averaged FADC slice entries for the supposed single-phe events
  TVector fAPedestalFADCSlices;            // Array containing the averaged FADC slice entries for the supposed pedestal   events

  TF1 *fSinglePheFit;                      //-> Single Phe Fit (Gaussians convoluted with Poisson) (needs to be a pointer, initialized later)

  UInt_t  fNumSinglePhes;
  UInt_t  fNumPedestals;

  Float_t fSinglePheCut;  

  Float_t fNumSinglePheLimit;

  Double_t  fLambda; 
  Double_t  fMu0; 
  Double_t  fMu1; 
  Double_t  fSigma0; 
  Double_t  fSigma1; 

  Double_t  fLambdaErr; 
  Double_t  fMu0Err; 
  Double_t  fMu1Err; 
  Double_t  fSigma0Err; 
  Double_t  fSigma1Err; 

  Double_t  fLambdaCheck;
  Double_t  fLambdaCheckErr;

  Double_t  fChisquare;
  Int_t     fNDF;
  Double_t  fProb;

  Double_t  fMeanPedestal;
  Double_t  fMeanPedestalErr;
  Double_t  fSigmaPedestal;
  Double_t  fSigmaPedestalErr;

  Float_t   fExtractSlices;

  Byte_t    fFlags;

  enum { kSinglePheFitOK, kPedestalFitOK };

  TPaveText *fFitLegend;                        //!
  TH1F      *fHSinglePheFADCSlices;             //!
  TH1F      *fHPedestalFADCSlices;              //!


  // 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 Init();

  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();
  
  // Setters
  void SetChargeNbins       ( const Int_t  bins =fgChargeNbins       )    { fChargeNbins       = bins;     }
  void SetChargeFirst       ( const Axis_t first=fgChargeFirst       )    { fChargeFirst       = first;    }
  void SetChargeLast        ( const Axis_t last =fgChargeLast        )    { fChargeLast        = last;     }
  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);
  
  // 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;      }

  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;
  
  // Draws
  void Draw(Option_t *opt="");

private:
  void DrawLegend();
  
  // Fits
public:
  enum FitFunc_t  { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele };

private:
  FitFunc_t fFitFunc;

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

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

  const static Double_t fNoWay = 10000000000.0;

  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 (mu1cat    < mu0)
        return fNoWay;

      if (sigma1cat < sigma0)
        return fNoWay;

      // if (sigma1cat < sigma1dyn)
      // return NoWay;

      //if (mu1cat < mu1dyn)
      // return NoWay;

      //      if (lambda1cat < lambda1dyn)
      // return NoWay;

      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 fNoWay;

      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma1 < sigma0)
        return fNoWay;
      
      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 fNoWay;
      
      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma1 < sigma0)
        return fNoWay;
      
      
      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 fNoWay;
      
      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma1 < sigma0)
        return fNoWay;
      
      
      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)  // Histograms from the Calibration Blind Pixel
};

#endif  /* MARS_MHCalibrationChargeBlindPix */