Changeset 2913 for trunk/MagicSoft/Mars

Timestamp:

01/26/04 15:09:56 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mcalib/MHCalibrationBlindPixel.cc (modified) (3 diffs)
• mcalib/MHCalibrationBlindPixel.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -14 +14 @@
      - TSpline5 now on stack. 
 
+   * mcalib/MHCalibrationBlindPixel.[h,cc]:
+     - force mu_0 in Blind Pixel Fit to be around 0 in fKPoisson4
 
  2004/01/26: Thomas Bretz

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

@@ -347 +347 @@
   const Stat_t   entries      = fHBlindPixelCharge->Integral("width");
   const Double_t lambda_guess = 0.5;
-  const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
+  const Double_t maximum_bin  = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin());
+  const Double_t norm         = entries/gkSq2Pi;
+
+  //
+  // Initialize the fit function
+  //
+  switch (fFitFunc)
+    {
+    case kEPoisson4:
+      //      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,5);
+      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
+      break;
+    case kEPoisson5:
+      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
+      break;
+    case kEPoisson6:
+      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
+      break;
+    case kEPolya:
+      break;
+
+    case kEMichele:
+      break;
+
+    default:
+      *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl;
+      return kFALSE;
+      break;
+    }
+  
+  const Double_t mu_0_guess = maximum_bin;
   const Double_t si_0_guess = 20.;
   const Double_t mu_1_guess = mu_0_guess + 50.;
   const Double_t si_1_guess = si_0_guess + si_0_guess;
-  const Double_t norm       = entries/gkSq2Pi;
-
-  //
-  // Initialize
+  //
+  // Initialize boundaries and start parameters
   //
   switch (fFitFunc)
@@ -360 +388 @@
       
     case kEPoisson4:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6);
+      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
+      fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
+      fSinglePheFit->SetParLimits(0,0.,1.);
+      fSinglePheFit->SetParLimits(1,-2.,2.);
+      fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
+      fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
+      fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
+      fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5);
+      break;
+    case kEPoisson5:
+    case kEPoisson6:
       fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
       fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
@@ -369 +407 @@
       fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
       fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
-      break;
-      
-    case kEPoisson5:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6);
-      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
-      fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
-      fSinglePheFit->SetParLimits(0,0.,1.);
-      fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
-      fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
-      fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
-      fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
-      fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
-      break;
-
-    case kEPoisson6:
-      fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6);
-      fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm);
-      fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area");
-      fSinglePheFit->SetParLimits(0,0.,1.);
-      fSinglePheFit->SetParLimits(1,min,(max-min)/1.5);
-      fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min)));
-      fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0);
-      fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5);
-      fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
-      break;
-
-    case kEPoisson7:
       break;
 

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h

@@ -132 +132 @@
     {
 
+      const Double_t qe = 0.247;
+
       Double_t lambda = par[0];  
+
+      Double_t b1 = par[1];
+      Double_t b2 = par[2];
       
       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;
@@ -193 +154 @@
       Double_t arg = 0.;
       
+      //      Double_t mu0 = 0.;
       Double_t mu0 = par[1];
       Double_t mu1 = par[2];