laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau1DHistPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <cstdlib>
16 #include <iostream>
17 #include <vector>
18 
19 #include "TAxis.h"
20 #include "TH1.h"
21 #include "TRandom.h"
22 #include "TSystem.h"
23 
24 #include "Lau1DHistPdf.hh"
25 #include "LauRandom.hh"
26 
27 class LauParameter;
28 
29 ClassImp(Lau1DHistPdf)
30 
31 
32 Lau1DHistPdf::Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa,
33  Bool_t useInterpolation, Bool_t fluctuateBins) :
34  LauAbsPdf(theVarName, std::vector<LauParameter*>(), minAbscissa, maxAbscissa),
35  hist_(hist ? dynamic_cast<TH1*>(hist->Clone()) : 0),
36  useInterpolation_(useInterpolation),
37  fluctuateBins_(fluctuateBins),
38  nBins_(0),
39  axisMin_(0.0),
40  axisMax_(0.0),
41  axisRange_(0.0)
42 {
43  // Constructor
44 
45  // Set the directory for the histogram
46  hist_->SetDirectory(0);
47 
48  // Save various attributes of the histogram
49  nBins_ = hist_->GetNbinsX();
50  TAxis* xAxis = hist_->GetXaxis();
51  axisMin_ = xAxis->GetXmin();
52  axisMax_ = xAxis->GetXmax();
53  axisRange_ = axisMax_ - axisMin_;
54 
55  // Check that axis range corresponds to range of abscissa
56  if (TMath::Abs(this->getMinAbscissa() - axisMin_)>1e-6) {
57  std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis minimum: " << axisMin_ <<
58  " does not correspond to abscissa minimum: " << this->getMinAbscissa() << "." << std::endl;
59  gSystem->Exit(EXIT_FAILURE);
60  }
61  if (TMath::Abs(this->getMaxAbscissa() - axisMax_)>1e-6) {
62  std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis maximum: " << axisMax_ <<
63  " does not correspond to abscissa maximum: " << this->getMaxAbscissa() << "." << std::endl;
64  gSystem->Exit(EXIT_FAILURE);
65  }
66 
67  // If the bins are to be fluctuated then do so now before
68  // calculating anything that depends on the bin content.
69  if (fluctuateBins) {
70  this->doBinFluctuation();
71  }
72 
73  // Calculate the PDF normalisation.
74  this->calcNorm();
75 
76  // And check it is OK.
77  this->checkNormalisation();
78 }
79 
81 {
82  // Destructor
83  delete hist_; hist_ = 0;
84 }
85 
86 Lau1DHistPdf::Lau1DHistPdf(const Lau1DHistPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
87 {
88  hist_ = other.hist_ ? dynamic_cast<TH1*>(other.hist_->Clone()) : 0;
91  this->setRandomFun(other.getRandomFun());
92  this->calcNorm();
93 }
94 
95 void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
96 {
97  if (this->heightUpToDate()) {
98  return;
99  }
100 
101  // Get the maximum height of the histogram
102  Int_t maxBin = hist_->GetMaximumBin();
103  Double_t height = hist_->GetBinContent(maxBin);
104  this->setMaxHeight(height);
105 }
106 
108 {
109  // Calculate the histogram normalisation.
110 
111  // Loop over the range to get the total area.
112  // Just sum the contributions up using 1e-3 increments of the range.
113  // Multiply the end result by dx.
114 
115  Double_t dx(1e-3*axisRange_);
116  Double_t area(0.0);
117 
118  Double_t x(axisMin_ + dx/2.0);
119  while (x > axisMin_ && x < axisMax_) {
120  area += this->interpolate(x);
121  x += dx;
122  }
123 
124  Double_t norm = area*dx;
125  this->setNorm(norm);
126 }
127 
129 {
130  Double_t dx(1e-3*axisRange_);
131  Double_t area(0.0);
132  Double_t areaNoNorm(0.0);
133 
134  Double_t x(axisMin_ + dx/2.0);
135  while (x > axisMin_ && x < axisMax_) {
136  area += this->interpolateNorm(x);
137  useInterpolation_ = kFALSE;
138  areaNoNorm += this->interpolate(x);
139  useInterpolation_ = kTRUE;
140  x += dx;
141  }
142  Double_t norm = area*dx;
143 
144  std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl;
145  std::cout << " : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl;
146  std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl;
147 }
148 
149 Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const
150 {
151  // Check that bin is in range [1 , nBins_]
152  if ((bin < 1) || (bin > nBins_)) {
153  return 0.0;
154  }
155  Double_t value = static_cast<Double_t>(hist_->GetBinContent(bin));
156  // protect against negative values
157  if ( value < 0.0 ) {
158  std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
159  value = 0.0;
160  }
161  return value;
162 }
163 
164 Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const
165 {
166  // Get the normalised interpolated value.
167  Double_t value = this->interpolate(x);
168  Double_t norm = this->getNorm();
169  return value/norm;
170 }
171 
173 {
174  // Check that the given abscissa is within the allowed range
175  if (!this->checkRange(abscissas)) {
176  gSystem->Exit(EXIT_FAILURE);
177  }
178 
179  // Get our abscissa
180  Double_t abscissa = abscissas[0];
181 
182  // Calculate the interpolated value
183  Double_t value = this->interpolate(abscissa);
184  this->setUnNormPDFVal(value);
185 }
186 
187 Double_t Lau1DHistPdf::interpolate(Double_t x) const
188 {
189  // This function returns the interpolated value of the histogram function
190  // for the given value of x by finding the adjacent bins and extrapolating
191  // using weights based on the inverse distance of the point from the adajcent
192  // bin centres.
193 
194  // Find the histogram bin
195  Int_t bin = hist_->FindBin(x);
196 
197  // Ask whether we want to do the interpolation, since this method is
198  // not reliable for low statistics histograms.
199  if (useInterpolation_ == kFALSE) {
200  return this->getBinHistValue(bin);
201  }
202 
203  // Find the bin centres (actual co-ordinate positions, not histogram indices)
204  Double_t cbinx = hist_->GetBinCenter(bin);
205 
206  // Find the adjacent bins
207  Double_t deltax = x - cbinx;
208 
209  Int_t bin_adj(0);
210  if (deltax > 0.0) {
211  bin_adj = bin + 1;
212  } else {
213  bin_adj = bin - 1;
214  }
215 
216 
217  Bool_t isBoundary(kFALSE);
218  if ( bin_adj > nBins_ || bin_adj < 1 ) {
219  isBoundary = kTRUE;
220  }
221 
222  // At the edges, do no interpolation, use entry in bin.
223  if (isBoundary == kTRUE) {
224  return this->getBinHistValue(bin);
225  }
226 
227  // Linear interpolation using inverse distance as weights.
228  // Find the adjacent bin centre
229  Double_t cbinx_adj = hist_->GetBinCenter(bin_adj);
230  Double_t deltax_adj = cbinx_adj - x;
231 
232  Double_t dx0 = TMath::Abs(deltax);
233  Double_t dx1 = TMath::Abs(deltax_adj);
234 
235  Double_t denom = dx0 + dx1;
236 
237  Double_t value0 = this->getBinHistValue(bin);
238  Double_t value1 = this->getBinHistValue(bin_adj);
239 
240  Double_t value = value0*dx1 + value1*dx0;
241  value /= denom;
242 
243  return value;
244 }
245 
247 {
248  TRandom* random = LauRandom::zeroSeedRandom();
249  for (Int_t bin(0); bin<nBins_; bin++) {
250  Double_t currentContent = hist_->GetBinContent(bin+1);
251  Double_t currentError = hist_->GetBinError(bin+1);
252  Double_t newContent = random->Gaus(currentContent,currentError);
253  if (newContent<0.0) {
254  hist_->SetBinContent(bin+1,0.0);
255  } else {
256  hist_->SetBinContent(bin+1,newContent);
257  }
258  }
259 }
260 
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
virtual void calcNorm()
Calculate the normalisation.
TH1 * hist_
The underlying histogram.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:454
TRandom * zeroSeedRandom()
Access the singleton random number generator with seed set from machine clock time (within +-1 sec)...
Definition: LauRandom.cc:30
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:349
Bool_t fluctuateBins_
Control boolean for performing the fluctuation of the histogram bin contents.
void checkNormalisation()
Check the normalisation calculation.
Bool_t useInterpolation_
Control boolean for using the linear interpolation.
Double_t axisMax_
The histogram axis maximum.
Class for defining a 1D histogram PDF.
Definition: Lau1DHistPdf.hh:31
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:410
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:213
void doBinFluctuation()
Fluctuate the histogram bin contents in accorance with their errors.
Double_t getBinHistValue(Int_t bin) const
Get the bin content from the histogram.
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:472
virtual ~Lau1DHistPdf()
Destructor.
Definition: Lau1DHistPdf.cc:80
Int_t nBins_
The number of bins in the histogram.
Lau1DHistPdf(const TString &theVarName, const TH1 *hist, Double_t minAbscissa, Double_t maxAbscissa, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE)
Constructor.
Definition: Lau1DHistPdf.cc:32
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:416
Class for defining the fit parameter objects.
Definition: LauParameter.hh:31
Double_t interpolateNorm(Double_t x) const
Perform the interpolation and divide by the normalisation.
File containing LauRandom namespace.
Double_t axisMin_
The histogram axis minimum.
Double_t axisRange_
The histogram axis range.
Double_t interpolate(Double_t x) const
Perform the interpolation (unnormalised)
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Definition: Lau1DHistPdf.cc:95
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:40
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:315
virtual Double_t getNorm() const
Retrieve the normalisation factor.
Definition: LauAbsPdf.hh:284
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:44
File containing declaration of Lau1DHistPdf class.