GausStepFitter.cxx
Go to the documentation of this file.
1 // GausStepFitter.cxx
2 
3 #include "GausStepFitter.h"
6 #include "TF1.h"
7 #include "TH1.h"
8 #include "TList.h"
9 #include <iostream>
10 
11 using std::cout;
12 using std::endl;
13 
14 //**********************************************************************
15 
16 GausStepFitter::GausStepFitter(double pos, double sigma, double height, Name fnam, Name fopt)
17 : m_pos(pos), m_sigma(sigma), m_height(height), m_fnam(fnam), m_fopt(fopt),
18  m_LogLevel(0) { }
19 
20 //**********************************************************************
21 
22 int GausStepFitter::fit(TH1* ph) const {
23  Name myname = "GausStepFitter::fit: ";
24  double sigma = m_sigma;
25  double height = m_height;
26  double sigfac = 2.0;
27  for ( int ifit=0; ifit<5; ++ifit ) {
28  TF1* pffix = gausTF1(height, m_pos, sigma, m_fnam);
29  double sigmax = 1.1*sigfac*sigma;
30  double sigmin = 0.9*sigma/sigfac;
31  if ( m_LogLevel >= 1 ) cout << myname << " Doing constrained fit " << ifit
32  << " with pos=" << m_pos << ", sigma=" << sigma
33  << " (" << sigmin << ", " << sigmax << ")" << endl;
34  pffix->SetParameter(2, sigma);
35  pffix->SetParLimits(2, sigmin, sigmax);
36  if ( m_height == 0.0 ) {
37  pffix->SetParLimits(0, 0.0, 1.e10); // Don't let height go negative.
38  } else if ( height > 0.0 ) {
39  pffix->SetParLimits(0, 0.1*height, 2.0*height); // Don't let height go negative.
40  } else if ( height < 0.0 ) {
41  pffix->SetParLimits(0, 2.0*height, 0.1*height); // Don't let height go positive.
42  }
43  Name fopt = m_fopt;
44  if ( m_LogLevel >= 2 ) fopt += "Q";
45  int fstat = quietHistFit(ph, pffix, fopt.c_str());
46  double signew = pffix->GetParameter(2);
47  bool atHiLimit = signew > 0.999*sigmax;
48  bool atLoLimit = signew < 1.001*sigmin;
49  if ( m_LogLevel >= 1 ) cout << myname << " status " << fstat << ", fit sigma=" << signew << endl;
50  if ( fstat == 0 && !atHiLimit && !atLoLimit ) {
51  if ( m_LogLevel >= 1 ) cout << myname << "Fit suceeded." << endl;
52  return 0;
53  }
54  // Otherwise we may get a crash when we try to view saved copy of histo:
55  ph->GetListOfFunctions()->Clear();
56  delete pffix;
57  if ( atLoLimit ) sigma /= sigfac;
58  else sigma *= sigfac;
59  }
60  if ( m_LogLevel >= 1 ) cout << myname << "Fit failed." << endl;
61  return 1;
62 }
63 
64 //**********************************************************************
GausStepFitter(double pos, double sigma, double height, Name fnam, Name fopt="WWS")
int quietHistFit(TH1 *ph, std::string fname, std::string fopt)
TF1 * gausTF1(double heightIn, double meanIn, double sigmaIn, std::string fname)
Definition: gausTF1.cxx:28
int fit(TH1 *ph) const
std::string Name
QTextStream & endl(QTextStream &s)