Changeset 8911 for trunk/MagicSoft/Mars/mhcalib

Timestamp:

06/02/08 10:36:33 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhcalib

Files:

• MHCalibrationChargeBlindPix.cc (modified) (3 diffs)
• MHCalibrationChargeBlindPix.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.cc

@@ -18 +18 @@
 !   Author(s): Markus Gaug 02/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2008
 !
 !
@@ -432 +432 @@
     {
     case kEPoisson4:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6);
+      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto4,rmin,rmax,6);
       rmin += 6.5;
       break;
     case kEPoisson5:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6);
+      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto5,rmin,rmax,6);
       rmin = 0.;
       break;
     case kEPoisson6:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6);
+      fSinglePheFit = new TF1("SinglePheFit",&PoissonKto6,rmin,rmax,6);
       break;
     case kEPolya:
-      fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8);
+      fSinglePheFit = new TF1("SinglePheFit",&Polya,rmin,rmax,8);
       break;
     case kEMichele:
-      fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9);
+      fSinglePheFit = new TF1("SinglePheFit",&FitFuncMichele,rmin,rmax,9);
       break;
     default:
@@ -964 +964 @@
 }
 
-
-
-
-
-
-
-
-
-
-
+Double_t MHCalibrationChargeBlindPix::FitFuncMichele(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.;
+}
+
+Double_t MHCalibrationChargeBlindPix::PoissonKto4(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;
+}
+
+Double_t MHCalibrationChargeBlindPix::PoissonKto5(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;
+}
+
+Double_t MHCalibrationChargeBlindPix::PoissonKto6(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;
+}
+
+Double_t MHCalibrationChargeBlindPix::Polya(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;
+}

trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeBlindPix.h

@@ -5 +5 @@
 #include "MHCalibrationPix.h"
 #endif
+
+#include <TMath.h>
 
 #ifndef ROOT_TF1
@@ -59 +61 @@
   Byte_t    fFlags;                          // Bit-field for the flags
   enum { kSinglePheFitOK, kPedestalFitOK };  // Possible bits to be set
+
+  static Double_t FitFuncMichele(Double_t *x, Double_t *par);
+  static Double_t PoissonKto4(Double_t *x, Double_t *par);
+  static Double_t PoissonKto5(Double_t *x, Double_t *par);
+  static Double_t PoissonKto6(Double_t *x, Double_t *par);
+  static Double_t Polya(Double_t *x, Double_t *par);
 
 public:
@@ -137 +145 @@
                            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
 };