#ifndef MARS_MHCalibrationChargeBlindPix #define MARS_MHCalibrationChargeBlindPix #ifndef MARS_MHCalibrationChargePix #include "MHCalibrationChargePix.h" #endif #ifndef ROOT_TMatrix #include #endif #ifndef ROOT_TF1 #include #endif class TH1F; class TF1; class TPaveText; class TText; class MRawEvtData; class MRawEvtPixelIter; class MCalibrationChargeBlindPix; class MExtractedSignalBlindPixel; class MHCalibrationChargeBlindPix : public MHCalibrationChargePix { 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 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=""); 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 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 }; // 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 offset = par[9]; 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. + offset; } 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 */