#ifndef MARS_MHCalibrationBlindPixel
#define MARS_MHCalibrationBlindPixel

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

class TArrayF;
class TArrayI;
class TH1F;
class TH1I;
class TF1;
class TPaveText;

class TMath;
class MParList;
class MHCalibrationBlindPixel : public MH
{
private:

  static const Int_t    fgBlindPixelChargeNbins;
  static const Int_t    fgBlindPixelTimeNbins;
  static const Axis_t   fgBlindPixelTimeFirst;
  static const Axis_t   fgBlindPixelTimeLast;
  static const Double_t fgBlindPixelElectronicAmp;
  static const Double_t fgBlindPixelElectronicAmpError;

  static const Int_t    fPSDNbins;
  static const Int_t    fPulserFrequency;
  
  TH1I* fHBlindPixelCharge;        // Histogram with the single Phe spectrum
  TH1F* fHBlindPixelTime;          // Variance of summed FADC slices
  TH1F* fHBlindPixelPSD;           // Power spectrum density of fHBlindPixelChargevsN

  TH1I* fHSinglePheFADCSlices;
  TH1I* fHPedestalFADCSlices;
  
  TF1 *fSinglePheFit;   
  TF1 *fTimeGausFit;    
  TF1 *fSinglePhePedFit;

  TArrayF* fPSDHiGain;           //-> Power spectrum density of fHiGains
  TArrayF* fPSDLoGain;           //-> Power spectrum density of fLoGains

  TH1I* fHPSD;                   //-> 
  TF1*  fPSDExpFit;              //->
  
  TArrayF *fHiGains;             //->
  TArrayF *fLoGains;             //->
  TArrayF *fChargeXaxis;         //
  TArrayF *fPSDXaxis;            //

  Float_t  fPSDProb;

  Int_t fTotalEntries;           // Number of entries
  Int_t fCurrentSize;
  
  Axis_t  fBlindPixelChargefirst;
  Axis_t  fBlindPixelChargelast; 
  
  void DrawLegend();
  void CreateChargeXaxis(Int_t n);
  void CreatePSDXaxis(Int_t n);
  void CutArrayBorder(TArrayF *array) const;
  void CutArrayBorder(TArrayI *array) const;

  TPaveText *fFitLegend;                  
  
  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  fChisquare; 
  Double_t  fProb;      
  Int_t     fNdf;       

  Double_t  fMeanTime; 
  Double_t  fMeanTimeErr; 
  Double_t  fSigmaTime; 
  Double_t  fSigmaTimeErr; 
  
  Double_t  fLambdaCheck;
  Double_t  fLambdaCheckErr;

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

  Byte_t   fFlags;

  enum     { kFitOK, kOscillating };
  
public:

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

  void Clear(Option_t *o="");  
  void Reset();
  
  Bool_t FillBlindPixelCharge(const Int_t q);
  Bool_t FillBlindPixelTime(const Float_t t);
  Bool_t FillGraphs(const Int_t qhi, const Int_t qlo);
  
  // Setters
  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; }
  
  // 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 Double_t GetChiSquare()      const { return fChisquare;      }
  const Double_t GetProb()           const { return fProb;           }  
  const Int_t    GetNdf()            const { return fNdf;            }   

  const Double_t GetMeanTime()       const { return fMeanTime;       }
  const Double_t GetMeanTimeErr()    const { return fMeanTimeErr;    }
  const Double_t GetSigmaTime()      const { return fSigmaTime;      }
  const Double_t GetSigmaTimeErr()   const { return fSigmaTimeErr;   }

  TH1I *GetHSinglePheFADCSlices()    { return fHSinglePheFADCSlices;  }
  TH1I *GetHPedestalFADCSlices()    { return fHPedestalFADCSlices;  }  
  
  const Bool_t IsFitOK()        const;
  const Bool_t IsOscillating();
  
  const TH1F *GetHBlindPixelPSD()       const  { return fHBlindPixelPSD;       }
  
  // Draws
  TObject *DrawClone(Option_t *option="") const;
  void Draw(Option_t *option="");

  
  // Fits
  enum FitFunc_t  { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele };

private:
  FitFunc_t fFitFunc;

public:
  Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="RL0+Q");
  Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+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);
  
  // Others
  void CutAllEdges();
  Bool_t CheckOscillations();

private:

  const static Double_t fNoWay = 10000000000.0;

  Bool_t InitFit(Axis_t min, Axis_t max);
  void ExitFit(TF1 *f);  
  
  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(MHCalibrationBlindPixel, 1)  // Histograms from the Calibration Blind Pixel
};

#endif  /* MARS_MHCalibrationBlindPixel */