laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
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 
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<LauAbsRValue*>(), 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 void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
87 {
88  if (this->heightUpToDate()) {
89  return;
90  }
91 
92  // Get the maximum height of the histogram
93  Int_t maxBin = hist_->GetMaximumBin();
94  Double_t height = hist_->GetBinContent(maxBin);
95  this->setMaxHeight(height);
96 }
97 
99 {
100  // Calculate the histogram normalisation.
101 
102  // Loop over the range to get the total area.
103  // Just sum the contributions up using 1e-3 increments of the range.
104  // Multiply the end result by dx.
105 
106  Double_t dx(1e-3*axisRange_);
107  Double_t area(0.0);
108 
109  Double_t x(axisMin_ + dx/2.0);
110  while (x > axisMin_ && x < axisMax_) {
111  area += this->interpolate(x);
112  x += dx;
113  }
114 
115  Double_t norm = area*dx;
116  this->setNorm(norm);
117 }
118 
120 {
121  Double_t dx(1e-3*axisRange_);
122  Double_t area(0.0);
123  Double_t areaNoNorm(0.0);
124 
125  Double_t x(axisMin_ + dx/2.0);
126  while (x > axisMin_ && x < axisMax_) {
127  area += this->interpolateNorm(x);
128  useInterpolation_ = kFALSE;
129  areaNoNorm += this->interpolate(x);
130  useInterpolation_ = kTRUE;
131  x += dx;
132  }
133  Double_t norm = area*dx;
134 
135  std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl;
136  std::cout << " : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl;
137  std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl;
138 }
139 
140 Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const
141 {
142  // Check that bin is in range [1 , nBins_]
143  if ((bin < 1) || (bin > nBins_)) {
144  return 0.0;
145  }
146  Double_t value = static_cast<Double_t>(hist_->GetBinContent(bin));
147  // protect against negative values
148  if ( value < 0.0 ) {
149  std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
150  value = 0.0;
151  }
152  return value;
153 }
154 
155 Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const
156 {
157  // Get the normalised interpolated value.
158  Double_t value = this->interpolate(x);
159  Double_t norm = this->getNorm();
160  return value/norm;
161 }
162 
164 {
165  // Check that the given abscissa is within the allowed range
166  if (!this->checkRange(abscissas)) {
167  gSystem->Exit(EXIT_FAILURE);
168  }
169 
170  // Get our abscissa
171  Double_t abscissa = abscissas[0];
172 
173  // Calculate the interpolated value
174  Double_t value = this->interpolate(abscissa);
175  this->setUnNormPDFVal(value);
176 }
177 
178 Double_t Lau1DHistPdf::interpolate(Double_t x) const
179 {
180  // This function returns the interpolated value of the histogram function
181  // for the given value of x by finding the adjacent bins and extrapolating
182  // using weights based on the inverse distance of the point from the adajcent
183  // bin centres.
184 
185  // Find the histogram bin
186  Int_t bin = hist_->FindFixBin(x);
187 
188  // Ask whether we want to do the interpolation, since this method is
189  // not reliable for low statistics histograms.
190  if (useInterpolation_ == kFALSE) {
191  return this->getBinHistValue(bin);
192  }
193 
194  // Find the bin centres (actual co-ordinate positions, not histogram indices)
195  Double_t cbinx = hist_->GetBinCenter(bin);
196 
197  // Find the adjacent bins
198  Double_t deltax = x - cbinx;
199 
200  Int_t bin_adj(0);
201  if (deltax > 0.0) {
202  bin_adj = bin + 1;
203  } else {
204  bin_adj = bin - 1;
205  }
206 
207 
208  Bool_t isBoundary(kFALSE);
209  if ( bin_adj > nBins_ || bin_adj < 1 ) {
210  isBoundary = kTRUE;
211  }
212 
213  // At the edges, do no interpolation, use entry in bin.
214  if (isBoundary == kTRUE) {
215  return this->getBinHistValue(bin);
216  }
217 
218  // Linear interpolation using inverse distance as weights.
219  // Find the adjacent bin centre
220  Double_t cbinx_adj = hist_->GetBinCenter(bin_adj);
221  Double_t deltax_adj = cbinx_adj - x;
222 
223  Double_t dx0 = TMath::Abs(deltax);
224  Double_t dx1 = TMath::Abs(deltax_adj);
225 
226  Double_t denom = dx0 + dx1;
227 
228  Double_t value0 = this->getBinHistValue(bin);
229  Double_t value1 = this->getBinHistValue(bin_adj);
230 
231  Double_t value = value0*dx1 + value1*dx0;
232  value /= denom;
233 
234  return value;
235 }
236 
238 {
239  TRandom* random = LauRandom::randomFun();
240  for (Int_t bin(0); bin<nBins_; bin++) {
241  Double_t currentContent = hist_->GetBinContent(bin+1);
242  Double_t currentError = hist_->GetBinError(bin+1);
243  Double_t newContent = random->Gaus(currentContent,currentError);
244  if (newContent<0.0) {
245  hist_->SetBinContent(bin+1,0.0);
246  } else {
247  hist_->SetBinContent(bin+1,newContent);
248  }
249  }
250 }
251 
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
virtual void calcNorm()
Calculate the normalisation.
Definition: Lau1DHistPdf.cc:98
TH1 * hist_
The underlying histogram.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:369
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
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:325
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 ~Lau1DHistPdf()
Destructor.
Definition: Lau1DHistPdf.cc:80
Int_t nBins_
The number of bins in the histogram.
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
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:86
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:29
virtual Double_t getNorm() const
Retrieve the normalisation factor.
Definition: LauAbsPdf.hh:202
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45
File containing declaration of Lau1DHistPdf class.