#ifndef MARS_MHCalibrationBlindPixel
#define MARS_MHCalibrationBlindPixel

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

#ifndef MARS_MHCalibrationConfig
#include "MHCalibrationConfig.h"
#endif

#ifndef ROOT_TH1
#include "TH1.h"
#endif

#ifndef ROOT_TH1F
#include "TH1F.h"
#endif

#ifndef ROOT_TF1
#include "TF1.h"
#endif

#ifndef ROOT_TPaveText
#include "TPaveText.h"
#endif

const Double_t NoWay = 10000000000.0;

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

  TH1F* fHBlindPixelCharge;        //-> Histogram with the single Phe spectrum
  TH1I* fHBlindPixelTime;          //-> Variance of summed FADC slices
  TH1I* fHBlindPixelChargevsN;     //-> Summed Charge vs. Event Nr. 
  
  TF1 fSinglePheFit;
  TF1 fTimeGausFit;  
  TF1 fSinglePhePedFit;

  Axis_t  fBlindPixelChargefirst;
  Axis_t  fBlindPixelChargelast;
  Int_t   fBlindPixelChargenbins;
  
  void ResetBin(Int_t i);
  void DrawLegend();

  TPaveText *fFitLegend;                  //!
  Bool_t fFitOK;  
  
  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;

public:

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

  Bool_t FillBlindPixelCharge(Float_t q)             { return fHBlindPixelCharge->Fill(q) > -1;  }  
  Bool_t FillBlindPixelTime(Int_t t)                 { return fHBlindPixelTime->Fill(t) > -1;  }
  Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1;  }  
  

  //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; }

  const Bool_t IsFitOK()                { return fFitOK; }

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

private:

  void InitFit(TF1& f, Axis_t min, Axis_t max);
  void ExitFit(TF1& f);  
  
  inline static Double_t fFitFuncMichele(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 NoWay;

      Double_t sigma0 = par[3];
      Double_t sigma1 = par[4];
      
      if (sigma1 < sigma0)
        return NoWay;
      
      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 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 NoWay;

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

#endif  /* MARS_MHCalibrationBlindPixel */