Changeset 2790

Timestamp:

01/13/04 17:29:27 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mcalib/MHCalibrationBlindPixel.cc (modified) (9 diffs)
• mhist/MHCamera.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -9 +9 @@
      - make Stat Window in option 'proj' bigger
      - add comment line at beginning of Draw()
+
+   * mcalib/MHCalibrationBlindPixel.[h.cc]
+     - add Function DrawClone()
+     - fix bug in calculation of integral (GetIntegral("width"))
+     - fit now with fixed normalization
+     - restrict parameters tighter
+     - declare fit as not valid, if less than 100 events 
+       in single photo-electron peak
 
 

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

@@ -27 +27 @@
 //  MHCalibrationBlindPixel                                                 //
 //                                                                          //
-//  Performs all the necessary fits to extract the mean number of photons   //
-//              out of the derived light flux                               //
+//  Performs all the Single Photo-Electron Fit to extract                   //
+//  the mean number of photons and to derive the light flux                 //
+//                                                                          //
+// The fit result is accepted under condition that:                         //
+// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)  //
+// 2) at least 100 events are in the single Photo-electron peak             // 
 //                                                                          //
 //////////////////////////////////////////////////////////////////////////////
@@ -83 +87 @@
     fHBlindPixelCharge->SetYTitle("Nr. of events");
     fHBlindPixelCharge->Sumw2();
+    fHBlindPixelCharge->SetDirectory(NULL);
 
     Axis_t tfirst = -0.5;
@@ -92 +97 @@
     fHBlindPixelTime->SetYTitle("Nr. of events");
     fHBlindPixelTime->Sumw2();
+    fHBlindPixelTime->SetDirectory(NULL);
 
     // We define a reasonable number and later enlarge it if necessary
@@ -101 +107 @@
     fHBlindPixelChargevsN->SetXTitle("Event Nr.");
     fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices");
-
-    fgSinglePheFitFunc = &gfKto8;
+    fHBlindPixelChargevsN->SetDirectory(NULL);
+
+    fgSinglePheFitFunc = &gfKto4;
     fgSinglePheFitNPar = 5;
 }
@@ -187 +194 @@
 
 }
+
+
+TObject *MHCalibrationBlindPixel::DrawClone(Option_t *opt) const
+{
+
+   //TVirtualPad *pad = gROOT->GetSelectedPad();
+   TVirtualPad *padsav = gPad;
+   //if (pad) pad->cd();
+
+   TObject *newobj = Clone();
+
+   if (!newobj)
+       return 0;
+
+   newobj->Draw(opt);
+
+   if (padsav)
+       padsav->cd();
+
+   return newobj;
+}
+  
 
 
@@ -317 +346 @@
   // otherwise the fit goes gaga because of high number of dimensions ...
   //
-  const Stat_t   entries      = fHBlindPixelCharge->Integral();
+  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 si_0_guess = 20.;
-  const Double_t mu_1_guess = mu_0_guess + 100.;
+  const Double_t mu_1_guess = mu_0_guess + 50.;
   const Double_t si_1_guess = si_0_guess + si_0_guess;
 
   fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1);
+  //  fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1);
+  //  fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess);
   fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries);
+  //  fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1");
   fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area");
-  fSinglePheFit->SetParLimits(0,0.,3.);
-  fSinglePheFit->SetParLimits(1,rmin,rmax);
-  fSinglePheFit->SetParLimits(2,rmin,rmax);
-  fSinglePheFit->SetParLimits(3,1.0,rmax-rmin);
-  fSinglePheFit->SetParLimits(4,1.7,rmax-rmin);
-  fSinglePheFit->SetParLimits(5,0.,1.5*entries);
+  fSinglePheFit->SetParLimits(0,0.,1.);
+  fSinglePheFit->SetParLimits(1,rmin,(rmax-rmin)/1.5);
+  fSinglePheFit->SetParLimits(2,(rmax-rmin)/2.,(rmax-0.05*(rmax-rmin)));
+  fSinglePheFit->SetParLimits(3,1.0,(rmax-rmin)/2.0);
+  fSinglePheFit->SetParLimits(4,1.0,(rmax-rmin)/2.5);
+  //  fSinglePheFit->SetParLimits(5,entries/gkSq2Pi,entries/gkSq2Pi);
+  //  fSinglePheFit->SetParLimits(5,0.,1.5*entries);
   //
   // Normalize the histogram to facilitate faster fitting of the area
-  // For speed reasons, FKto8 is normalized to Sqrt(2 pi).
+  // For speed reasons, gfKto4 is normalized to Sqrt(2 pi).
   // Therefore also normalize the histogram to that value 
   //
-  // ROOT gives us another nice example of user-unfriendly behavior:
-  // Although the normalization of the function fSinglePhe and the 
-  // Histogram fHBlindPixelCharge agree (!!), the fit does not normalize correctly INTERNALLY
-  // in the fitting procedure !!!
-  // 
-  // This has to do with the fact that the internal function histogramming 
-  // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does
-  // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST)
-  // 
-  // So, WE have to adapt to that internal flaw of ROOT:
-  //
-  //  const Int_t  npx     = fSinglePheFit->GetNpx();
-  //  const Int_t  bins    = fHBlindPixelCharge->GetXaxis()->GetLast()-fHBlindPixelCharge->GetXaxis()->GetFirst();
   //  fHBlindPixelCharge->Scale(gkSq2Pi*(float)bins/npx/entries);
-
-  // 
-  // we need this, otherwise, ROOT does not calculate the area correctly
-  // don't ask me why it does not behave correctly, it's one of the nasty
-  // mysteries of ROOT which takes you a whole day to find out :-)
-  //
-  //  fSinglePheFit->SetNpx(fChargenbins);  
-
-  fHBlindPixelCharge->Fit(fSinglePheFit,opt);
+  // fHBlindPixelCharge->Scale(gkSq2Pi/entries);
+  Float_t  norm   = entries/gkSq2Pi;
+  //           norm  *= (Float_t)fSinglePheFit->GetNpx()/(Float_t)fBlindPixelChargenbins;
+  fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1);
+
   fHBlindPixelCharge->Fit(fSinglePheFit,opt);
 
@@ -384 +400 @@
 
   Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2);
-  cout << "Parameter0: " << fSinglePhePedFit->GetParameter(0) << endl;
-  cout << "Parameter2: " << fSinglePhePedFit->GetParameter(2) << endl;
-  cout << "Pedarea: " << pedarea << endl;
-  cout << "entries: " << entries << endl;
-  fLambdaCheck     = TMath::Log((double)entries/pedarea);
+  fLambdaCheck     = TMath::Log(entries/pedarea);
   fLambdaCheckErr  = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0)
                      + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2);
@@ -396 +408 @@
   *fLog << "DoF: " << fNdf << endl;
   *fLog << "Probability: " << fProb << endl;
-  *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax);
-
-  //
-  // The fit result is accepted under condition
-  // The fit result is accepted under condition
-  // The Probability is greater than gkProbLimit (default 0.001 == 99.7%)
+
+  //
+  // The fit result is accepted under condition that: 
+  // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%)
+  // 2) at least 100 events are in the single Photo-electron peak
   //
   if (fProb < gkProbLimit) 
     {
-      *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl;
+      *fLog << err << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl;
       fFitOK = kFALSE;
       return kFALSE;
     }
-  
+
+  Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries;
+  
+  if (contSinglePhe < 100.) 
+    {
+      *fLog << err << "Statistics is too low: Only " << contSinglePhe 
+            << " in the first electron peak " << endl;
+      fFitOK = kFALSE;
+      return kFALSE;
+    } 
+  else
+    *fLog << GetDescriptor() << ": " << contSinglePhe 
+          << " in first electron peak " << endl;
   
   fFitOK = kTRUE;
@@ -420 +443 @@
 {
 
-  Int_t nbins = 50;
+  Int_t nbins = 30;
 
   CutEdges(fHBlindPixelCharge,nbins);

trunk/MagicSoft/Mars/mhist/MHCamera.cc

@@ -375 +375 @@
 //   'content'     Display the relative content aligned to GetMaximum() and
 //                 GeMinimum() ((val-min)/(max-min))
+//   'proj'        Display the y-projection of the histogram
 //
 void MHCamera::Draw(Option_t *option)
@@ -620 +621 @@
       CreateProjection();
       opt.ReplaceAll("proj", "");
+      Float_t he = gStyle->GetStatH();
+      Float_t wi = gStyle->GetStatH();
+      gStyle->SetStatH(0.4);
+      gStyle->SetStatW(0.25);
       fYProj->Paint(opt);
+      gStyle->SetStatH(he);      
+      gStyle->SetStatW(wi);      
       return;
     }