Changeset 3766


Ignore:
Timestamp:
04/16/04 14:14:01 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc

    r3741 r3766  
    267267}
    268268
    269 // --------------------------------------------------------------------------
    270 //
    271 // Default Reset(), can be overloaded.
    272 //
    273 // Executes:
    274 // - Clear()
    275 // - fHGausHist.Reset()
    276 // - fEvents.Set(0)
    277 //
    278 void MHGausEvents::Reset()
    279 {
    280 
    281   Clear();
    282   fHGausHist.Reset();
    283   fEvents.Set(0);
    284 
    285 }
    286 
    287 // --------------------------------------------------------------------------
    288 //
    289 // Default InitBins, can be overloaded.
    290 //
    291 // Executes:
    292 // - fHGausHist.SetBins(fNbins,fFirst,fLast)
    293 //
    294 void MHGausEvents::InitBins()
    295 {
    296   fHGausHist.SetBins(fNbins,fFirst,fLast);
    297 }
    298 
    299 // --------------------------------------------------------------------------
    300 //
    301 // Executes:
    302 // - FillArray()
    303 // - FillHist()
    304 //
    305 Bool_t MHGausEvents::FillHistAndArray(const Float_t f)
    306 {
    307 
    308   FillArray(f);
    309   return FillHist(f);
    310 }
    311 
    312 // --------------------------------------------------------------------------
    313 //
    314 // Fills fHGausHist with f
    315 // Returns kFALSE, if overflow or underflow occurred, else kTRUE
    316 //
    317 Bool_t MHGausEvents::FillHist(const Float_t f)
    318 {
    319   return fHGausHist.Fill(f) > -1;
    320 }
    321 
    322 // --------------------------------------------------------------------------
    323 //
    324 // Fill fEvents with f
    325 // If size of fEvents is 0, initializes it to 512
    326 // If size of fEvents is smaller than fCurrentSize, double the size
    327 // Increase fCurrentSize by 1
    328 //
    329 void MHGausEvents::FillArray(const Float_t f)
    330 {
    331   if (fEvents.GetSize() == 0)
    332     fEvents.Set(512);
    333 
    334   if (fCurrentSize >= fEvents.GetSize())
    335     fEvents.Set(fEvents.GetSize()*2);
    336  
    337   fEvents.AddAt(f,fCurrentSize++);
    338 }
    339 
    340 // --------------------------------------------------------------------------
    341 //
    342 // Set Excluded bit from outside
    343 //
    344 void MHGausEvents::SetExcluded(const Bool_t b)
    345 {
    346     b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
    347 }
    348 
    349 
    350 const Double_t MHGausEvents::GetChiSquare()  const
    351 {
    352   return ( fFGausFit ? fFGausFit->GetChisquare() : 0.);
    353 }
    354 
    355 const Int_t MHGausEvents::GetNdf() const
    356 {
    357   return ( fFGausFit ? fFGausFit->GetNDF() : 0);
    358 }
    359 
    360 
    361 const Double_t MHGausEvents::GetExpChiSquare()  const
    362 {
    363   return ( fFExpFit ? fFExpFit->GetChisquare() : 0.);
    364 }
    365 
    366 
    367 const Int_t MHGausEvents::GetExpNdf()  const
    368 {
    369   return ( fFExpFit ? fFExpFit->GetNDF() : 0);
    370 }
    371 
    372 
    373 const Double_t MHGausEvents::GetExpProb()  const
    374 {
    375   return ( fFExpFit ? fFExpFit->GetProb() : 0.);
    376 }
    377 
    378 
    379 const Double_t MHGausEvents::GetOffset()  const
    380 {
    381   return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.);
    382 }
    383 
    384 
    385 const Double_t MHGausEvents::GetSlope()  const
    386 {
    387   return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.);
    388 }
    389 
    390 
    391 const Bool_t MHGausEvents::IsEmpty() const
    392 {
    393     return !(fHGausHist.GetEntries());
    394 }
    395 
    396 
    397 const Bool_t MHGausEvents::IsFourierSpectrumOK() const
    398 {
    399   return TESTBIT(fFlags,kFourierSpectrumOK);
    400 }
    401 
    402 
    403 const Bool_t MHGausEvents::IsGausFitOK() const
    404 {
    405   return TESTBIT(fFlags,kGausFitOK);
    406 }
    407 
    408 
    409 const Bool_t MHGausEvents::IsExpFitOK() const
    410 {
    411   return TESTBIT(fFlags,kExpFitOK);
    412 }
    413 
    414 const Bool_t MHGausEvents::IsExcluded() const
    415 {
    416   return TESTBIT(fFlags,kExcluded);
    417 }
    418 
    419 
    420 // -------------------------------------------------------------------
    421 //
    422 // The flag setters are to be used ONLY for Monte-Carlo!!
    423 //
    424 void  MHGausEvents::SetGausFitOK(const Bool_t b)
    425 {
    426   b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK);
    427 
    428 }
    429 // -------------------------------------------------------------------
    430 //
    431 // The flag setters are to be used ONLY for Monte-Carlo!!
    432 //
    433 void  MHGausEvents::SetExpFitOK(const Bool_t b)
    434 {
    435  
    436   b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK); 
    437 }
    438 
    439 // -------------------------------------------------------------------
    440 //
    441 // The flag setters are to be used ONLY for Monte-Carlo!!
    442 //
    443 void  MHGausEvents::SetFourierSpectrumOK(const Bool_t b)
    444 {
    445 
    446   b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK);   
     269
     270// -----------------------------------------------------------------------------
     271//
     272// Bypasses the Gauss fit by taking mean and RMS from the histogram
     273//
     274// Errors are determined in the following way:
     275// MeanErr  = RMS / Sqrt(entries)
     276// SigmaErr = RMS / (2.*Sqrt(entries) )
     277//
     278void MHGausEvents::BypassFit()
     279{
     280
     281  const Stat_t entries = fHGausHist.GetEntries();
     282 
     283  if (entries <= 0.)
     284    {
     285      *fLog << warn << GetDescriptor()
     286            << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
     287      return;
     288    }
     289 
     290  fMean     = fHGausHist.GetMean();
     291  fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
     292  fSigma    = fHGausHist.GetRMS() ;
     293  fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
     294}
     295
     296
     297
     298// --------------------------------------------------------------------------
     299//
     300// - Set fPixId to id
     301//
     302// Add id to names and titles of:
     303// - fHGausHist
     304//
     305void MHGausEvents::ChangeHistId(const Int_t id)
     306{
     307
     308  fPixId = id;
     309
     310  fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
     311  fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
     312
     313  fName  = Form("%s%d", fName.Data(),  id);
     314  fTitle = Form("%s%d", fTitle.Data(), id);
     315
     316}
     317
     318// -----------------------------------------------------------------------------
     319//
     320// Create the x-axis for the event graph
     321//
     322Float_t *MHGausEvents::CreateEventXaxis(Int_t n)
     323{
     324
     325  Float_t *xaxis = new Float_t[n];
     326
     327  if (fEventFrequency)
     328    for (Int_t i=0;i<n;i++)
     329      xaxis[i] = (Float_t)i/fEventFrequency;
     330  else
     331    for (Int_t i=0;i<n;i++)
     332      xaxis[i] = (Float_t)i;
     333
     334  return xaxis;
     335                 
    447336}
    448337
    … …  
    527416}
    528417
    529 // -------------------------------------------------------------------
    530 //
    531 // Fit fGausHist with a Gaussian after stripping zeros from both ends
    532 // and rebinned to the number of bins specified in fBinsAfterStripping
    533 //
    534 // The fit results are retrieved and stored in class-own variables. 
    535 //
    536 // A flag IsGausFitOK() is set according to whether the fit probability
    537 // is smaller or bigger than fProbLimit, whether the NDF is bigger than
    538 // fNDFLimit and whether results are NaNs.
    539 //
    540 Bool_t MHGausEvents::FitGaus(Option_t *option, const Double_t xmin, const Double_t xmax)
    541 {
    542 
    543   if (IsGausFitOK())
    544     return kTRUE;
    545 
    546   //
    547   // First, strip the zeros from the edges which contain only zeros and rebin
    548   // to about fBinsAfterStripping bins.
    549   //
    550   // (ATTENTION: The Chisquare method is more sensitive,
    551   // the _less_ bins, you have!)
    552   //
    553   StripZeros(&fHGausHist,fBinsAfterStripping);
    554  
    555   //
    556   // Get the fitting ranges
    557   //
    558   Axis_t rmin = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()) : xmin;
    559   Axis_t rmax = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast())  : xmax;
    560 
    561   //
    562   // First guesses for the fit (should be as close to reality as possible,
    563   //
    564   const Stat_t   entries     = fHGausHist.Integral("width");
    565   const Double_t mu_guess    = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
    566   const Double_t sigma_guess = fHGausHist.GetRMS();
    567   const Double_t area_guess  = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;
    568 
    569   fFGausFit = new TF1("GausFit","gaus",rmin,rmax);
    570 
    571   if (!fFGausFit)
    572     {
    573     *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
    574           << "in pixel: " << fPixId << endl;
    575     return kFALSE;
    576     }
    577  
    578   fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
    579   fFGausFit->SetParNames("Area","#mu","#sigma");
    580   fFGausFit->SetParLimits(0,0.,area_guess*1.5);
    581   fFGausFit->SetParLimits(1,rmin,rmax);
    582   fFGausFit->SetParLimits(2,0.,rmax-rmin);
    583   fFGausFit->SetRange(rmin,rmax);
    584 
    585   fHGausHist.Fit(fFGausFit,option);
    586 
    587 
    588   fMean     = fFGausFit->GetParameter(1);
    589   fSigma    = fFGausFit->GetParameter(2);
    590   fMeanErr  = fFGausFit->GetParError(1);
    591   fSigmaErr = fFGausFit->GetParError(2);
    592   fProb     = fFGausFit->GetProb();
    593   //
    594   // The fit result is accepted under condition:
    595   // 1) The results are not nan's
    596   // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
    597   // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
    598   //
    599   if (   TMath::IsNaN(fMean)
    600       || TMath::IsNaN(fMeanErr)
    601       || TMath::IsNaN(fProb)   
    602       || TMath::IsNaN(fSigma)
    603       || TMath::IsNaN(fSigmaErr)
    604       || fFGausFit->GetNDF() < fNDFLimit
    605       || fProb < fProbLimit )
    606     return kFALSE;
    607  
    608   SetGausFitOK(kTRUE);
    609   return kTRUE;
    610 }
    611 
    612 // -----------------------------------------------------------------------------
    613 //
    614 // If flag IsGausFitOK() is set (histogram already successfully fitted),
    615 // returns kTRUE
    616 //
    617 // If both fMean and fSigma are still zero, call FitGaus()
    618 //
    619 // Repeats the Gauss fit in a smaller range, defined by:
    620 //
    621 // min = GetMean() - fBlackoutLimit * GetSigma();
    622 // max = GetMean() + fPickupLimit   * GetSigma();
    623 //
    624 // The fit results are retrieved and stored in class-own variables. 
    625 //
    626 // A flag IsGausFitOK() is set according to whether the fit probability
    627 // is smaller or bigger than fProbLimit, whether the NDF is bigger than
    628 // fNDFLimit and whether results are NaNs.
    629 //
    630 Bool_t MHGausEvents::RepeatFit(const Option_t *option)
    631 {
    632 
    633   if (IsGausFitOK())
    634     return kTRUE;
    635 
    636   if ((fMean == 0.) && (fSigma == 0.))
    637     return FitGaus();
    638 
    639   //
    640   // Get new fitting ranges
    641   //
    642   Axis_t rmin = fMean - fBlackoutLimit * fSigma;
    643   Axis_t rmax = fMean + fPickupLimit   * fSigma;
    644 
    645   Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
    646   Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
    647 
    648   fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
    649 
    650   fHGausHist.Fit(fFGausFit,option);
    651 
    652   fMean     = fFGausFit->GetParameter(1);
    653   fMeanErr  = fFGausFit->GetParameter(2);
    654   fSigma    = fFGausFit->GetParError(1) ;
    655   fSigmaErr = fFGausFit->GetParError(2) ;
    656   fProb     = fFGausFit->GetProb()      ;     
    657 
    658   //
    659   // The fit result is accepted under condition:
    660   // 1) The results are not nan's
    661   // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
    662   // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
    663   //
    664   if (   TMath::IsNaN ( fMean     )
    665       || TMath::IsNaN ( fMeanErr  )
    666       || TMath::IsNaN ( fProb     )   
    667       || TMath::IsNaN ( fSigma    )
    668       || TMath::IsNaN ( fSigmaErr )
    669       || fFGausFit->GetNDF() < fNDFLimit
    670       || fProb < fProbLimit )
    671     return kFALSE;
    672  
    673   SetGausFitOK(kTRUE);
    674   return kTRUE;
    675 
    676 }
    677 
    678 // -----------------------------------------------------------------------------
    679 //
    680 // Bypasses the Gauss fit by taking mean and RMS from the histogram
    681 //
    682 // Errors are determined in the following way:
    683 // MeanErr  = RMS / Sqrt(entries)
    684 // SigmaErr = RMS / (2.*Sqrt(entries) )
    685 //
    686 void MHGausEvents::BypassFit()
    687 {
    688 
    689   const Stat_t entries = fHGausHist.GetEntries();
    690  
    691   if (entries <= 0.)
    692     {
    693       *fLog << warn << GetDescriptor()
    694             << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
    695       return;
    696     }
    697  
    698   fMean     = fHGausHist.GetMean();
    699   fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
    700   fSigma    = fHGausHist.GetRMS() ;
    701   fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
    702 }
    703 
    704 // -----------------------------------------------------------------------------------
    705 //
    706 // A default print
    707 //
    708 void MHGausEvents::Print(const Option_t *o) const
    709 {
    710  
    711   *fLog << all                                                        << endl;
    712   *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId     << endl;
    713   *fLog << all << "Mean: "             << GetMean()                   << endl;
    714   *fLog << all << "Sigma: "            << GetSigma()                  << endl;
    715   *fLog << all << "Chisquare: "        << GetChiSquare()              << endl;
    716   *fLog << all << "DoF: "              << GetNdf()                    << endl;
    717   *fLog << all << "Probability: "      << GetProb()                   << endl;
    718   *fLog << all                                                        << endl;
    719  
    720 }
    721 
    722 // --------------------------------------------------------------------------
    723 //
    724 // - Set fPixId to id
    725 //
    726 // Add id to names and titles of:
    727 // - fHGausHist
    728 //
    729 void MHGausEvents::ChangeHistId(const Int_t id)
    730 {
    731 
    732   fPixId = id;
    733 
    734   fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
    735   fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
    736 
    737   fName  = Form("%s%d", fName.Data(),  id);
    738   fTitle = Form("%s%d", fTitle.Data(), id);
    739 
    740 }
    741 
    742 // --------------------------------------------------------------------------
    743 //
    744 // Re-normalize the results, has to be overloaded
    745 //
    746 void  MHGausEvents::Renorm()
    747 {
    748 }
    749 
    750418// ----------------------------------------------------------------------------------
    751419//
    … …  
    783451}
    784452
    785 // -----------------------------------------------------------------------------
    786 //
    787 // Create the x-axis for the event graph
    788 //
    789 Float_t *MHGausEvents::CreateEventXaxis(Int_t n)
    790 {
    791 
    792   Float_t *xaxis = new Float_t[n];
    793 
    794   if (fEventFrequency)
    795     for (Int_t i=0;i<n;i++)
    796       xaxis[i] = (Float_t)i/fEventFrequency;
    797   else
    798     for (Int_t i=0;i<n;i++)
    799       xaxis[i] = (Float_t)i;
    800 
    801   return xaxis;
    802                  
    803 }
    804453
    805454// -----------------------------------------------------------------------------
    … …  
    942591}
    943592
    944 const Double_t MHGausEvents::GetPickup() const
     593
     594// --------------------------------------------------------------------------
     595//
     596// Fill fEvents with f
     597// If size of fEvents is 0, initializes it to 512
     598// If size of fEvents is smaller than fCurrentSize, double the size
     599// Increase fCurrentSize by 1
     600//
     601void MHGausEvents::FillArray(const Float_t f)
     602{
     603  if (fEvents.GetSize() == 0)
     604    fEvents.Set(512);
     605
     606  if (fCurrentSize >= fEvents.GetSize())
     607    fEvents.Set(fEvents.GetSize()*2);
     608 
     609  fEvents.AddAt(f,fCurrentSize++);
     610}
     611
     612
     613// --------------------------------------------------------------------------
     614//
     615// Fills fHGausHist with f
     616// Returns kFALSE, if overflow or underflow occurred, else kTRUE
     617//
     618Bool_t MHGausEvents::FillHist(const Float_t f)
     619{
     620  return fHGausHist.Fill(f) > -1;
     621}
     622
     623// --------------------------------------------------------------------------
     624//
     625// Executes:
     626// - FillArray()
     627// - FillHist()
     628//
     629Bool_t MHGausEvents::FillHistAndArray(const Float_t f)
     630{
     631
     632  FillArray(f);
     633  return FillHist(f);
     634}
     635
     636// -------------------------------------------------------------------
     637//
     638// Fit fGausHist with a Gaussian after stripping zeros from both ends
     639// and rebinned to the number of bins specified in fBinsAfterStripping
     640//
     641// The fit results are retrieved and stored in class-own variables. 
     642//
     643// A flag IsGausFitOK() is set according to whether the fit probability
     644// is smaller or bigger than fProbLimit, whether the NDF is bigger than
     645// fNDFLimit and whether results are NaNs.
     646//
     647Bool_t MHGausEvents::FitGaus(Option_t *option, const Double_t xmin, const Double_t xmax)
     648{
     649
     650  if (IsGausFitOK())
     651    return kTRUE;
     652
     653  //
     654  // First, strip the zeros from the edges which contain only zeros and rebin
     655  // to about fBinsAfterStripping bins.
     656  //
     657  // (ATTENTION: The Chisquare method is more sensitive,
     658  // the _less_ bins, you have!)
     659  //
     660  StripZeros(&fHGausHist,fBinsAfterStripping);
     661 
     662  //
     663  // Get the fitting ranges
     664  //
     665  Axis_t rmin = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()) : xmin;
     666  Axis_t rmax = (xmin==0.) && (xmax==0.) ? fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast())  : xmax;
     667
     668  //
     669  // First guesses for the fit (should be as close to reality as possible,
     670  //
     671  const Stat_t   entries     = fHGausHist.Integral("width");
     672  const Double_t mu_guess    = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin());
     673  const Double_t sigma_guess = fHGausHist.GetRMS();
     674  const Double_t area_guess  = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess;
     675
     676  fFGausFit = new TF1("GausFit","gaus",rmin,rmax);
     677
     678  if (!fFGausFit)
     679    {
     680    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
     681          << "in pixel: " << fPixId << endl;
     682    return kFALSE;
     683    }
     684 
     685  fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
     686  fFGausFit->SetParNames("Area","#mu","#sigma");
     687  fFGausFit->SetParLimits(0,0.,area_guess*1.5);
     688  fFGausFit->SetParLimits(1,rmin,rmax);
     689  fFGausFit->SetParLimits(2,0.,rmax-rmin);
     690  fFGausFit->SetRange(rmin,rmax);
     691
     692  fHGausHist.Fit(fFGausFit,option);
     693
     694
     695  fMean     = fFGausFit->GetParameter(1);
     696  fSigma    = fFGausFit->GetParameter(2);
     697  fMeanErr  = fFGausFit->GetParError(1);
     698  fSigmaErr = fFGausFit->GetParError(2);
     699  fProb     = fFGausFit->GetProb();
     700  //
     701  // The fit result is accepted under condition:
     702  // 1) The results are not nan's
     703  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
     704  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
     705  //
     706  if (   TMath::IsNaN(fMean)
     707      || TMath::IsNaN(fMeanErr)
     708      || TMath::IsNaN(fProb)   
     709      || TMath::IsNaN(fSigma)
     710      || TMath::IsNaN(fSigmaErr)
     711      || fFGausFit->GetNDF() < fNDFLimit
     712      || fProb < fProbLimit )
     713    return kFALSE;
     714 
     715  SetGausFitOK(kTRUE);
     716  return kTRUE;
     717}
     718
     719// -------------------------------------------------------------------------------
     720//
     721// Return the number of "blackout" events, which are events with values higher
     722// than fBlackoutLimit sigmas from the mean
     723//
     724//
     725const Double_t MHGausEvents::GetBlackout() const
    945726{
    946727 
    … …  
    948729    return -1.;
    949730
     731  const Int_t first = fHGausHist.GetXaxis()->GetFirst();
     732  const Int_t last  = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);
     733
     734  if (first >= last)
     735    return 0.;
     736 
     737  return fHGausHist.Integral(first, last, "width");
     738}
     739
     740const Double_t MHGausEvents::GetChiSquare()  const
     741{
     742  return ( fFGausFit ? fFGausFit->GetChisquare() : 0.);
     743}
     744
     745
     746const Double_t MHGausEvents::GetExpChiSquare()  const
     747{
     748  return ( fFExpFit ? fFExpFit->GetChisquare() : 0.);
     749}
     750
     751
     752const Int_t MHGausEvents::GetExpNdf()  const
     753{
     754  return ( fFExpFit ? fFExpFit->GetNDF() : 0);
     755}
     756
     757
     758const Double_t MHGausEvents::GetExpProb()  const
     759{
     760  return ( fFExpFit ? fFExpFit->GetProb() : 0.);
     761}
     762
     763
     764const Int_t MHGausEvents::GetNdf() const
     765{
     766  return ( fFGausFit ? fFGausFit->GetNDF() : 0);
     767}
     768
     769const Double_t MHGausEvents::GetOffset()  const
     770{
     771  return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.);
     772}
     773
     774
     775// -------------------------------------------------------------------------------
     776//
     777// Return the number of "pickup" events, which are events with values higher
     778// than fPickupLimit sigmas from the mean
     779//
     780//
     781const Double_t MHGausEvents::GetPickup() const
     782{
     783 
     784  if ((fMean == 0.) && (fSigma == 0.))
     785    return -1.;
     786
    950787  const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
    951788  const Int_t last  = fHGausHist.GetXaxis()->GetLast();
    … …  
    957794}
    958795
    959 const Double_t MHGausEvents::GetBlackout() const
    960 {
    961  
     796
     797const Double_t MHGausEvents::GetSlope()  const
     798{
     799  return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.);
     800}
     801
     802// --------------------------------------------------------------------------
     803//
     804// Default InitBins, can be overloaded.
     805//
     806// Executes:
     807// - fHGausHist.SetBins(fNbins,fFirst,fLast)
     808//
     809void MHGausEvents::InitBins()
     810{
     811  fHGausHist.SetBins(fNbins,fFirst,fLast);
     812}
     813
     814const Bool_t MHGausEvents::IsEmpty() const
     815{
     816    return !(fHGausHist.GetEntries());
     817}
     818
     819
     820const Bool_t MHGausEvents::IsExcluded() const
     821{
     822  return TESTBIT(fFlags,kExcluded);
     823}
     824
     825
     826const Bool_t MHGausEvents::IsExpFitOK() const
     827{
     828  return TESTBIT(fFlags,kExpFitOK);
     829}
     830
     831const Bool_t MHGausEvents::IsFourierSpectrumOK() const
     832{
     833  return TESTBIT(fFlags,kFourierSpectrumOK);
     834}
     835
     836
     837const Bool_t MHGausEvents::IsGausFitOK() const
     838{
     839  return TESTBIT(fFlags,kGausFitOK);
     840}
     841
     842
     843// -----------------------------------------------------------------------------------
     844//
     845// A default print
     846//
     847void MHGausEvents::Print(const Option_t *o) const
     848{
     849 
     850  *fLog << all                                                        << endl;
     851  *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId     << endl;
     852  *fLog << all << "Mean: "             << GetMean()                   << endl;
     853  *fLog << all << "Sigma: "            << GetSigma()                  << endl;
     854  *fLog << all << "Chisquare: "        << GetChiSquare()              << endl;
     855  *fLog << all << "DoF: "              << GetNdf()                    << endl;
     856  *fLog << all << "Probability: "      << GetProb()                   << endl;
     857  *fLog << all                                                        << endl;
     858 
     859}
     860
     861// --------------------------------------------------------------------------
     862//
     863// Re-normalize the results, has to be overloaded
     864//
     865void  MHGausEvents::Renorm()
     866{
     867}
     868
     869// -----------------------------------------------------------------------------
     870//
     871// If flag IsGausFitOK() is set (histogram already successfully fitted),
     872// returns kTRUE
     873//
     874// If both fMean and fSigma are still zero, call FitGaus()
     875//
     876// Repeats the Gauss fit in a smaller range, defined by:
     877//
     878// min = GetMean() - fBlackoutLimit * GetSigma();
     879// max = GetMean() + fPickupLimit   * GetSigma();
     880//
     881// The fit results are retrieved and stored in class-own variables. 
     882//
     883// A flag IsGausFitOK() is set according to whether the fit probability
     884// is smaller or bigger than fProbLimit, whether the NDF is bigger than
     885// fNDFLimit and whether results are NaNs.
     886//
     887Bool_t MHGausEvents::RepeatFit(const Option_t *option)
     888{
     889
     890  if (IsGausFitOK())
     891    return kTRUE;
     892
    962893  if ((fMean == 0.) && (fSigma == 0.))
    963     return -1.;
    964 
    965   const Int_t first = fHGausHist.GetXaxis()->GetFirst();
    966   const Int_t last  = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);
    967 
    968   if (first >= last)
    969     return 0.;
    970  
    971   return fHGausHist.Integral(first, last, "width");
    972 }
    973 
    974 
     894    return FitGaus();
     895
     896  //
     897  // Get new fitting ranges
     898  //
     899  Axis_t rmin = fMean - fBlackoutLimit * fSigma;
     900  Axis_t rmax = fMean + fPickupLimit   * fSigma;
     901
     902  Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
     903  Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
     904
     905  fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
     906
     907  fHGausHist.Fit(fFGausFit,option);
     908
     909  fMean     = fFGausFit->GetParameter(1);
     910  fMeanErr  = fFGausFit->GetParameter(2);
     911  fSigma    = fFGausFit->GetParError(1) ;
     912  fSigmaErr = fFGausFit->GetParError(2) ;
     913  fProb     = fFGausFit->GetProb()      ;     
     914
     915  //
     916  // The fit result is accepted under condition:
     917  // 1) The results are not nan's
     918  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
     919  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
     920  //
     921  if (   TMath::IsNaN ( fMean     )
     922      || TMath::IsNaN ( fMeanErr  )
     923      || TMath::IsNaN ( fProb     )   
     924      || TMath::IsNaN ( fSigma    )
     925      || TMath::IsNaN ( fSigmaErr )
     926      || fFGausFit->GetNDF() < fNDFLimit
     927      || fProb < fProbLimit )
     928    return kFALSE;
     929 
     930  SetGausFitOK(kTRUE);
     931  return kTRUE;
     932
     933}
     934
     935// --------------------------------------------------------------------------
     936//
     937// Default Reset(), can be overloaded.
     938//
     939// Executes:
     940// - Clear()
     941// - fHGausHist.Reset()
     942// - fEvents.Set(0)
     943//
     944void MHGausEvents::Reset()
     945{
     946
     947  Clear();
     948  fHGausHist.Reset();
     949  fEvents.Set(0);
     950
     951}
     952
     953// --------------------------------------------------------------------------
     954//
     955// Set Excluded bit from outside
     956//
     957void MHGausEvents::SetExcluded(const Bool_t b)
     958{
     959    b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
     960}
     961
     962
     963// -------------------------------------------------------------------
     964//
     965// The flag setters are to be used ONLY for Monte-Carlo!!
     966//
     967void  MHGausEvents::SetExpFitOK(const Bool_t b)
     968{
     969 
     970  b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK); 
     971}
     972
     973// -------------------------------------------------------------------
     974//
     975// The flag setters are to be used ONLY for Monte-Carlo!!
     976//
     977void  MHGausEvents::SetFourierSpectrumOK(const Bool_t b)
     978{
     979
     980  b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK);   
     981}
     982
     983
     984// -------------------------------------------------------------------
     985//
     986// The flag setters are to be used ONLY for Monte-Carlo!!
     987//
     988void  MHGausEvents::SetGausFitOK(const Bool_t b)
     989{
     990  b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK);
     991
     992}
Note: See TracChangeset for help on using the changeset viewer.