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