Changeset 3636 for trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc

Timestamp:

04/03/04 17:27:50 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (modified) (22 diffs)

trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc

@@ -40 +40 @@
 //     It is created with a default name and title and resides in directory NULL.
 //     - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist
-//       (e.g. with the function TH1F::SetBins(..) )
+//       (e.g. by setting the variables fNbins, fFirst, fLast and calling the function 
+//       InitBins() or by directly calling fHGausHist.SetBins(..) )
 //     - The histogram is filled with the functions FillHist() or FillHistAndArray(). 
 //     - The histogram can be fitted with the function FitGaus(). This involves stripping 
@@ -48 +49 @@
 //       The NDF is compared to fNDFLimit and a check is made whether results are NaNs. 
 //       Accordingly, a flag IsGausFitOK() is set.
+//     - One can repeat the fit within a given amount of sigmas from the previous mean 
+//       (specified by the variable fPickupLimit) with the function RepeatFit()
+//     - One can completely skip the fitting to set mean, sigma and its errors directly 
+//       from the histograms with the function BypassFit()
 // 
 //  2) The TArrayF fEvents:
@@ -96 +101 @@
 const Float_t  MHGausEvents::fgProbLimit            = 0.001;
 const Int_t    MHGausEvents::fgNDFLimit             = 2;
+const Float_t  MHGausEvents::fgPickupLimit          = 5.;
 const Int_t    MHGausEvents::fgPowerProbabilityBins = 20;
 const Int_t    MHGausEvents::fgBinsAfterStripping   = 40;
@@ -105 +111 @@
 // - the default Probability Limit for fProbLimit  (fgProbLimit)
 // - the default NDF Limit for fNDFLimit  (fgNDFLimit)
+// - the default number for fPickupLimit  (fgPickupLimit)
 // - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping)
 // - the default name of the fHGausHist   ("HGausHist")
 // - the default title of the fHGausHist  ("Histogram of Events with Gaussian Distribution")
 // - the default directory of the fHGausHist (NULL)
+// - the default for fNbins (100)
+// - the default for fFirst (0.)
+// - the default for fLast  (100.)
 //
 // Initializes:
@@ -121 +131 @@
       fPowerSpectrum(NULL),
       fGraphEvents(NULL), fGraphPowerSpectrum(NULL),
+      fNbins(100), fFirst(0.), fLast(100.), fPixId(-1),
       fHGausHist(), fEvents(0),
       fFGausFit(NULL), fFExpFit(NULL)
@@ -133 +144 @@
   SetProbLimit();
   SetNDFLimit();
+  SetPickupLimit();
   SetBinsAfterStripping();
 
@@ -183 +195 @@
 // --------------------------------------------------------------------------
 //
+// Default Clear(), can be overloaded.
+//
 // Sets:
 // - all other pointers to NULL
-// - all variables to 0., except fEventFrequency
+// - all variables to 0., except fPixId to -1 and keep fEventFrequency
 // - all flags to kFALSE
 // 
@@ -197 +211 @@
   SetExpFitOK         ( kFALSE );
   SetFourierSpectrumOK( kFALSE );
+  SetExcluded         ( kFALSE );
 
   fMean              = 0.;
@@ -205 +220 @@
   
   fCurrentSize       = 0;
+  fPixId             = -1;
 
   if (fHPowerProbability)
@@ -247 +263 @@
 
 // --------------------------------------------------------------------------
+//
+// Default Reset(), can be overloaded.
 //
 // Executes:
@@ -264 +282 @@
 // --------------------------------------------------------------------------
 //
+// Default InitBins, can be overloaded.
+//
+// Executes:
+// - fHGausHist.SetBins(fNbins,fFirst,fLast)
+//
+void MHGausEvents::InitBins()
+{
+  fHGausHist.SetBins(fNbins,fFirst,fLast);
+}
+
+// --------------------------------------------------------------------------
+//
 // Executes:
 // - FillArray()
@@ -303 +333 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Set Excluded bit from outside
+//
+void MHGausEvents::SetExcluded(const Bool_t b)
+{
+    b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
+}
+
 
 const Double_t MHGausEvents::GetChiSquare()  const 
@@ -366 +405 @@
 {
   return TESTBIT(fFlags,kExpFitOK);
+}
+
+const Bool_t MHGausEvents::IsExcluded() const
+{
+  return TESTBIT(fFlags,kExcluded);
 }
 
@@ -409 +453 @@
   if (fFExpFit)
     return;
+
+  if (fEvents.GetSize() < 8)
+    {
+      *fLog << warn << "Cannot create Fourier spectrum in pixel: " << fPixId 
+            << ". Number of events smaller than 8 " << endl;
+      return;
+    }
+  
 
   //
@@ -424 +476 @@
   // This cuts only the non-used zero's, but MFFT will later cut the rest
   MArray::StripZeros(fEvents);
+
+  if (fEvents.GetSize() < 8)
+    {
+      *fLog << warn << "Cannot create Fourier spectrum. " << endl;
+      *fLog << warn << "Number of events (after stripping of zeros) is smaller than 8 " 
+            << "in pixel: " << fPixId << endl;
+      return;
+    }
 
   MFFT fourier;
@@ -506 +566 @@
   if (!fFGausFit) 
     {
-    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl;
+    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " 
+          << "in pixel: " << fPixId << endl;
     return kFALSE;
     }
@@ -528 +589 @@
   // The fit result is accepted under condition:
   // 1) The results are not nan's
-  // 2) The NDF is not smaller than fNDFLimit (5)
-  // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
+  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
+  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
   //
   if (   TMath::IsNaN(fMean) 
@@ -545 +606 @@
 
 // -----------------------------------------------------------------------------
+//
+// If flag IsGausFitOK() is set (histogram already successfully fitted), 
+// returns kTRUE
+// 
+// If both fMean and fSigma are still zero, call FitGaus() 
+// 
+// Repeats the Gauss fit in a smaller range, defined by: 
+// 
+// min = GetMean() - fPickupLimit * GetSigma();
+// max = GetMean() + fPickupLimit * GetSigma();
+//
+// The fit results are retrieved and stored in class-own variables.  
+//
+// A flag IsGausFitOK() is set according to whether the fit probability 
+// is smaller or bigger than fProbLimit, whether the NDF is bigger than 
+// fNDFLimit and whether results are NaNs.
+//
+Bool_t MHGausEvents::RepeatFit(const Option_t *option)
+{
+
+  if (IsGausFitOK())
+    return kTRUE;
+
+  if ((fMean == 0.) && (fSigma == 0.))
+    return FitGaus();
+
+  //
+  // Get new fitting ranges
+  //
+  Axis_t rmin = fMean - fPickupLimit * fSigma;
+  Axis_t rmax = fMean + fPickupLimit * fSigma;
+
+  Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
+  Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
+
+  fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
+
+  fHGausHist.Fit(fFGausFit,option);
+
+  fMean     = fFGausFit->GetParameter(1);
+  fMeanErr  = fFGausFit->GetParameter(2);
+  fSigma    = fFGausFit->GetParError(1) ; 
+  fSigmaErr = fFGausFit->GetParError(2) ; 
+  fProb     = fFGausFit->GetProb()      ;      
+
+  //
+  // The fit result is accepted under condition:
+  // 1) The results are not nan's
+  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
+  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
+  //
+  if (   TMath::IsNaN ( fMean     ) 
+      || TMath::IsNaN ( fMeanErr  )
+      || TMath::IsNaN ( fProb     )    
+      || TMath::IsNaN ( fSigma    )
+      || TMath::IsNaN ( fSigmaErr ) 
+      || fFGausFit->GetNDF() < fNDFLimit 
+      || fProb < fProbLimit )
+    return kFALSE;
+  
+  SetGausFitOK(kTRUE);
+  return kTRUE;
+
+}
+
+// -----------------------------------------------------------------------------
 // 
 // Bypasses the Gauss fit by taking mean and RMS from the histogram
@@ -559 +686 @@
   if (entries <= 0.)
     {
-      *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0" << endl;
+      *fLog << warn << GetDescriptor() 
+            << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
       return;
     }
   
-  SetMean     ( fHGausHist.GetMean() );
-  SetMeanErr  ( fHGausHist.GetRMS() / TMath::Sqrt(entries) );
-  SetSigma    ( fHGausHist.GetRMS() );
-  SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2. );
+  fMean     = fHGausHist.GetMean();
+  fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
+  fSigma    = fHGausHist.GetRMS() ;
+  fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
 }
 
@@ -577 +705 @@
   
   *fLog << all                                                        << endl;
-  *fLog << all << "Results of the Gauss Fit: "                        << endl;
+  *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId     << endl;
   *fLog << all << "Mean: "             << GetMean()                   << endl;
   *fLog << all << "Sigma: "            << GetSigma()                  << endl;
@@ -585 +713 @@
   *fLog << all                                                        << endl;
   
+}
+
+// --------------------------------------------------------------------------
+//
+// - Set fPixId to id
+//
+// Add id to names and titles of:
+// - fHGausHist
+//
+void MHGausEvents::ChangeHistId(Int_t id)
+{
+
+  fPixId = id;
+
+  fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
+  fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
+
+  fName  = Form("%s%d", fName.Data(),  id);
+  fTitle = Form("%s%d", fTitle.Data(), id);
+
 }
 
@@ -766 +914 @@
 }
 
-
+const Double_t MHGausEvents::GetPickup() const 
+{
+  
+  if ((fMean == 0.) && (fSigma == 0.))
+    return -1.;
+
+  const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
+  const Int_t last  = fHGausHist.GetXaxis()->GetLast();
+
+  if (first >= last)
+    return 0.;
+  
+  return fHGausHist.Integral(first, last, "width");
+}
+
+