Index: trunk/MagicSoft/Mars/Changelog
===================================================================
--- trunk/MagicSoft/Mars/Changelog	(revision 2827)
+++ trunk/MagicSoft/Mars/Changelog	(revision 2828)
@@ -4,4 +4,14 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2004/01/16: Markus Gaug
+
+   * mcalib/MCalibrationFits.h 
+     - removed and incorporated in MHCalibrationBlindPixel.h
+
+   * mcalib/MHCalibrationBlindPixel.[h,cc] 
+     - incorporate Fit functions
+     - fixed a bug due to which DrawClone crashed when class was used 
+       in a compiled macro
 
  2004/01/16: Abelardo Moralejo
Index: trunk/MagicSoft/Mars/mcalib/MCalibrationFits.h
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MCalibrationFits.h	(revision 2827)
+++ 	(revision )
@@ -1,368 +1,0 @@
-#ifndef MARS_MCalibrationFits
-#define MARS_MCalibrationFits
-
-#ifndef ROOT_TMath
-#include <TMath.h>
-#endif
-
-#define GIMMEABREAK     10000000000.0
-
-inline Double_t gfKto4(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 GIMMEABREAK;
-
-  Double_t sigma0 = par[3];
-  Double_t sigma1 = par[4];
-
-  if (sigma1 < sigma0)
-    return GIMMEABREAK;
-
-
-  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 Double_t gfKto5(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 GIMMEABREAK;
-
-  Double_t sigma0 = par[3];
-  Double_t sigma1 = par[4];
-
-  if (sigma1 < sigma0)
-    return GIMMEABREAK;
-
-
-  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 Double_t gfKto6(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 GIMMEABREAK;
-
-  Double_t sigma0 = par[3];
-  Double_t sigma1 = par[4];
-
-  if (sigma1 < sigma0)
-    return GIMMEABREAK;
-
-
-  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 Double_t gfKto7(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 GIMMEABREAK;
-
-  Double_t sigma0 = par[3];
-  Double_t sigma1 = par[4];
-
-  if (sigma1 < sigma0)
-    return GIMMEABREAK;
-
-
-  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 mu7 = (7.*mu1)-(6.*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 sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*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;
-  Double_t lambda7 = lambda6*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;
-  
-  // k=7:
-  arg = (x[0] - mu7)/sigma7;
-  sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
-  
-  return TMath::Exp(-1.*lambda)*par[5]*sum;
-
-};
-
-
-inline Double_t gfKto8(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 GIMMEABREAK;
-
-  Double_t sigma0 = par[3];
-  Double_t sigma1 = par[4];
-
-  if (sigma1 < sigma0)
-    return GIMMEABREAK;
-
-
-  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 mu7 = (7.*mu1)-(6.*mu0);
-  Double_t mu8 = (8.*mu1)-(7.*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 sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
-  Double_t sigma8 = TMath::Sqrt((8.*sigma1*sigma1) - (7.*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;
-  Double_t lambda7 = lambda6*lambda;
-  Double_t lambda8 = lambda7*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;
-  
-  // k=7:
-  arg = (x[0] - mu7)/sigma7;
-  sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
-  
-  // k=8:
-  arg = (x[0] - mu8)/sigma8;
-  sum += 0.000024801587315*lambda8*TMath::Exp(-0.5*arg*arg)/sigma7;
-  
-  return TMath::Exp(-1.*lambda)*par[5]*sum;
-
-};
-
-#endif  /* MARS_MCalibrationFits */
-
Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h
===================================================================
--- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h	(revision 2827)
+++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h	(revision 2828)
@@ -25,4 +25,6 @@
 #include "TPaveText.h"
 #endif
+
+const Double_t NoWay = 10000000000.0;
 
 class TMath;
@@ -50,7 +52,4 @@
   Bool_t fFitOK;  
   
-  BlindPixelFitFunc fgSinglePheFitFunc;     // In the beginning, 
-  Int_t     fgSinglePheFitNPar;             // we want to be flexible using different functions
-
   Double_t  fLambda; 
   Double_t  fMu0; 
@@ -81,4 +80,6 @@
   MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL);
   ~MHCalibrationBlindPixel();
+
+  typedef Double_t (*BlindPixelFitFunc)(Double_t *, Double_t *);
 
   Bool_t FillBlindPixelCharge(Float_t q)             { return fHBlindPixelCharge->Fill(q) > -1;  }  
@@ -118,5 +119,5 @@
   Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+Q");
 
-  void ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par=5);
+  void ChangeFitFunc(BlindPixelFitFunc fitfunc, Int_t par=6);
 
   void CutAllEdges();
@@ -127,4 +128,202 @@
   Bool_t IsFitOK() { return fFitOK; }
 
+private:
+  
+  BlindPixelFitFunc fgSinglePheFitFunc;     //! In the beginning, 
+  Int_t             fgSinglePheFitNPar;     //! we want to be flexible using different functions
+  
+  inline static Double_t gfKto4(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 gfKto5(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 gfKto6(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) 
 };