laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DHistDPPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 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 <iostream>
16 
17 #include "TAxis.h"
18 #include "TH2.h"
19 #include "TRandom.h"
20 
21 #include "Lau2DHistDPPdf.hh"
22 #include "LauKinematics.hh"
23 #include "LauRandom.hh"
24 #include "LauVetoes.hh"
25 
26 ClassImp(Lau2DHistDPPdf)
27 
28 
29 Lau2DHistDPPdf::Lau2DHistDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
30  Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP) :
31  Lau2DAbsHistDPPdf(kinematics,vetoes,useUpperHalfOnly,squareDP),
32  hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
33  minX_(0.0), maxX_(0.0), minY_(0.0), maxY_(0.0), rangeX_(0.0), rangeY_(0.0),
34  binXWidth_(0.0), binYWidth_(0.0),
35  invBinXWidth_(0.0), invBinYWidth_(0.0),
36  nBinsX_(0), nBinsY_(0),
37  norm_(0.0),
38  useInterpolation_(useInterpolation)
39 {
40  // For square Dalitz plots, the co-ordinates must be m', theta'
41  // The input data to the fit is still in m13^2, m23^2 co-ords, and a
42  // transformation is applied to go from one co-ordinate system to the
43  // other.
44 
45  // Save various attributes of the histogram
46  // (axis ranges, number of bins, areas)
47  if ( hist_ ) {
48  TAxis* xAxis = hist_->GetXaxis();
49  minX_ = xAxis->GetXmin();
50  maxX_ = xAxis->GetXmax();
51 
52  TAxis* yAxis = hist_->GetYaxis();
53  minY_ = yAxis->GetXmin();
54  maxY_ = yAxis->GetXmax();
55 
56  nBinsX_ = hist_->GetNbinsX();
57  nBinsY_ = hist_->GetNbinsY();
58  } else {
59  minX_ = getKinematics()->getm13SqMin();
60  maxX_ = getKinematics()->getm13SqMax();
61 
62  minY_ = getKinematics()->getm23SqMin();
63  maxY_ = getKinematics()->getm23SqMax();
64 
65  nBinsX_ = 1;
66  nBinsY_ = 1;
67  }
68 
69  rangeX_ = maxX_ - minX_;
70  rangeY_ = maxY_ - minY_;
71  if (nBinsX_ > 0) { binXWidth_ = TMath::Abs(rangeX_)/(nBinsX_*1.0); }
72  if (nBinsY_ > 0) { binYWidth_ = TMath::Abs(rangeY_)/(nBinsY_*1.0); }
73  if (binXWidth_ > 1e-10) { invBinXWidth_ = 1.0/binXWidth_; }
74  if (binYWidth_ > 1e-10) { invBinYWidth_ = 1.0/binYWidth_; }
75 
76  // If the bins are to be fluctuated then do so now before
77  // calculating anything that depends on the bin content.
78  if (fluctuateBins) {
79  this->doBinFluctuation(hist_);
80  }
81 
82  // Calculate the PDF normalisation.
83  this->calcHistNorm();
84 
85  // And check it is OK.
86  this->checkNormalisation();
87 
88  // Also obtain the maximum height
89  this->calcMaxHeight(hist_);
90 }
91 
93 {
94  // Destructor
95  delete hist_; hist_ = 0;
96 }
97 
99 {
100  // Calculate the histogram normalisation. We must integrate
101  // over the allowed kinematic DP region.
102 
103  // Loop over the total x and y range to get the total area
104  // under x and y. Just sum the contributions up using 1e-3 increments
105  // of the range in x and y. Multiply the end result by dx and dy.
106 
107  const Double_t axisMin = TMath::Min( minX_, minY_ );
108  const Double_t axisMax = TMath::Max( maxX_, maxY_ );
109  const Double_t axisRange = axisMax - axisMin;
110 
111  const Double_t dx(1e-3 * axisRange);
112  const Double_t dy(dx);
113 
114  Double_t area(0.0);
115 
116  Double_t x(axisMin + dx/2.0);
117  while (x < axisMax) {
118  Double_t y(axisMin + dy/2.0);
119  while (y < axisMax) {
120  area += this->interpolateXY(x,y);
121  y += dy;
122  } // y while loop
123  x += dx;
124  } // x while loop
125 
126  norm_ = area*dx*dy;
127 
128  std::cout << "INFO in Lau2DHistDPPdf::calcHistNorm : Norm = " << norm_ << ", bX*bY = " << binXWidth_ << "*" << binYWidth_ << " = " << binXWidth_*binYWidth_ << std::endl;
129 }
130 
131 Double_t Lau2DHistDPPdf::getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
132 {
133  if (xBinNo < 0) {
134  xBinNo = 0;
135  } else if (xBinNo >= nBinsX_) {
136  return 0.0;
137  }
138 
139  if (yBinNo < 0) {
140  yBinNo = 0;
141  } else if (yBinNo >= nBinsY_) {
142  return 0.0;
143  }
144 
145  if ( hist_ == 0 ) {
146  return 1.0;
147  }
148 
149  Double_t value = hist_->GetBinContent(xBinNo+1, yBinNo+1);
150  return value;
151 }
152 
153 Double_t Lau2DHistDPPdf::interpolateXYNorm(Double_t x, Double_t y) const
154 {
155  // Get the normalised interpolated value.
156  Double_t value = this->interpolateXY(x,y);
157  return value/norm_;
158 
159 }
160 
161 Double_t Lau2DHistDPPdf::interpolateXY(Double_t x, Double_t y) const
162 {
163  // This function returns the interpolated value of the histogram function
164  // for the given values of x and y by finding the adjacent bins and extrapolating
165  // using weights based on the inverse distance of the point from the adajcent
166  // bin centres.
167  // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
168 
169  // If we're only using one half then flip co-ordinates
170  // appropriately for conventional or square DP
171  getUpperHalf(x,y);
172 
173  // First ask whether the point is inside the kinematic region.
174  if (withinDPBoundaries(x,y) == kFALSE) {return 0.0;}
175 
176  // Update the kinematics to the position of interest.
177  updateKinematics(x,y);
178 
179  // Check that we're not inside a veto
180  Bool_t vetoOK(kTRUE);
181  if (getVetoes()) {
182  vetoOK = getVetoes()->passVeto(getKinematics());
183  }
184  if (vetoOK == kFALSE) {return 0.0;}
185 
186  // Find the 2D histogram bin for x and y
187  Int_t i = TMath::FloorNint((x - minX_)*invBinXWidth_);
188  Int_t j = TMath::FloorNint((y - minY_)*invBinYWidth_);
189 
190  if (i >= nBinsX_) {i = nBinsX_ - 1;}
191  if (j >= nBinsY_) {j = nBinsY_ - 1;}
192 
193  // Ask whether we want to do the interpolation, since this method is
194  // not reliable for low statistics histograms.
195  if (useInterpolation_ == kFALSE) {
196  return this->getBinHistValue(i,j);
197  }
198 
199  // Find the bin centres (actual co-ordinate positions, not histogram indices)
200  Double_t cbinx = (i+0.5)*binXWidth_ + minX_;
201  Double_t cbiny = (j+0.5)*binYWidth_ + minY_;
202 
203  // If bin centres are outside kinematic region, do not extrapolate
204  if (withinDPBoundaries(cbinx, cbiny) == kFALSE) {
205  return this->getBinHistValue(i,j);
206  }
207 
208  // Find the adjacent bins
209  Double_t deltax = x - cbinx;
210  Double_t deltay = y - cbiny;
211 
212  Int_t i_adj(0), j_adj(0);
213  if (deltax > 0.0) {
214  i_adj = i + 1;
215  } else {
216  i_adj = i - 1;
217  }
218  if (deltay > 0.0) {
219  j_adj = j + 1;
220  } else {
221  j_adj = j - 1;
222  }
223 
224  Bool_t isXBoundary(kFALSE), isYBoundary(kFALSE);
225 
226  Double_t value(0.0);
227 
228  if (i_adj >= nBinsX_ || i_adj < 0) {isYBoundary = kTRUE;}
229  if (j_adj >= nBinsY_ || j_adj < 0) {isXBoundary = kTRUE;}
230 
231  // In the corners, do no interpolation. Use entry in bin (i,j)
232  if (isXBoundary == kTRUE && isYBoundary == kTRUE) {
233 
234  value = this->getBinHistValue(i,j);
235 
236  } else if (isXBoundary == kTRUE && isYBoundary == kFALSE) {
237 
238  // Find the adjacent x bin centre
239  Double_t cbinx_adj = Double_t(i_adj+0.5)*binXWidth_ + minX_;
240 
241  if (withinDPBoundaries(cbinx_adj, y) == kFALSE) {
242 
243  // The adjacent bin is outside the DP range. Don't extrapolate.
244  value = this->getBinHistValue(i,j);
245 
246  } else {
247 
248  Double_t dx0 = TMath::Abs(x - cbinx);
249  Double_t dx1 = TMath::Abs(cbinx_adj - x);
250  Double_t inter_denom = dx0 + dx1;
251 
252  Double_t value1 = this->getBinHistValue(i,j);
253  Double_t value2 = this->getBinHistValue(i_adj,j);
254 
255  value = (value1*dx1 + value2*dx0)/inter_denom;
256 
257  }
258 
259  } else if (isYBoundary == kTRUE && isXBoundary == kFALSE) {
260 
261  // Find the adjacent y bin centre
262  Double_t cbiny_adj = Double_t(j_adj+0.5)*binYWidth_ + minY_;
263 
264  if (withinDPBoundaries(x, cbiny_adj) == kFALSE) {
265 
266  // The adjacent bin is outside the DP range. Don't extrapolate.
267  value = this->getBinHistValue(i,j);
268 
269  } else {
270 
271  Double_t dy0 = TMath::Abs(y - cbiny);
272  Double_t dy1 = TMath::Abs(cbiny_adj - y);
273  Double_t inter_denom = dy0 + dy1;
274 
275  Double_t value1 = this->getBinHistValue(i,j);
276  Double_t value2 = this->getBinHistValue(i,j_adj);
277 
278  value = (value1*dy1 + value2*dy0)/inter_denom;
279 
280  }
281 
282  } else {
283 
284  // 2D linear interpolation using inverse distance as weights.
285  // Find the adjacent x and y bin centres
286  Double_t cbinx_adj = Double_t(i_adj+0.5)*binXWidth_ + minX_;
287  Double_t cbiny_adj = Double_t(j_adj+0.5)*binYWidth_ + minY_;
288 
289  if (withinDPBoundaries(cbinx_adj, cbiny_adj) == kFALSE) {
290 
291  // The adjacent bin is outside the DP range. Don't extrapolate.
292  value = this->getBinHistValue(i,j);
293 
294  } else {
295 
296  Double_t dx0 = TMath::Abs(x - cbinx);
297  Double_t dx1 = TMath::Abs(cbinx_adj - x);
298  Double_t dy0 = TMath::Abs(y - cbiny);
299  Double_t dy1 = TMath::Abs(cbiny_adj - y);
300 
301  Double_t inter_denom = (dx0 + dx1)*(dy0 + dy1);
302 
303  Double_t value1 = this->getBinHistValue(i,j);
304  Double_t value2 = this->getBinHistValue(i_adj,j);
305  Double_t value3 = this->getBinHistValue(i,j_adj);
306  Double_t value4 = this->getBinHistValue(i_adj,j_adj);
307 
308  value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
309  value /= inter_denom;
310  }
311 
312  }
313 
314  return value;
315 
316 }
317 
319 {
320  // Loop over the total x and y range to get the total area
321  // under x and y. Just sum the contributions up using 1e-3 increments
322  // of the range in x and y. Multiply the end result by dx and dy.
323 
324  const Double_t axisMin = TMath::Min( minX_, minY_ );
325  const Double_t axisMax = TMath::Max( maxX_, maxY_ );
326  const Double_t axisRange = axisMax - axisMin;
327 
328  const Double_t dx(1e-3 * axisRange);
329  const Double_t dy(dx);
330 
331  Double_t area(0.0);
332  Double_t areaNoNorm(0.0);
333 
334  // Preserve the value of a variable we change temporarily
335  Bool_t interpolate = useInterpolation_;
336 
337  Double_t x(axisMin + dx/2.0);
338  while (x < axisMax) {
339  Double_t y(axisMin + dy/2.0);
340  while (y < axisMax) {
341  area += this->interpolateXYNorm(x,y);
342  useInterpolation_ = kFALSE;
343  areaNoNorm += this->interpolateXY(x,y);
344  useInterpolation_ = kTRUE;
345  y += dy;
346  } // y while loop
347  x += dx;
348  } // x while loop
349  Double_t norm = area*dx*dy;
350 
351  useInterpolation_ = interpolate; //Restore old value
352 
353  std::cout << "INFO in Lau2DHistDPPdf::checkNormalisation : Area = " << area << ", dx = " << dx << ", dy = " << dy << ", dx*dy = " << dx*dy << std::endl;
354  std::cout << " : Area with no norm = " << areaNoNorm << "*dx*dy = " << areaNoNorm*dx*dy << std::endl;
355  std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl;
356 }
357 
Double_t minY_
The histogram y-axis minimum.
Double_t norm_
The histogram normalisation.
Double_t invBinYWidth_
The histogram y-axis inverse bin width.
Double_t maxX_
The histogram x-axis maximum.
Double_t interpolateXYNorm(Double_t x, Double_t y) const
Perform the interpolation and divide by the normalisation.
Abstract base class for defining a variation across a 2D DP based on a histogram. ...
Bool_t withinDPBoundaries(Double_t x, Double_t y) const
Check whether the given co-ordinates are within the kinematic boundary.
Class for defining a 2D DP histogram PDF.
const LauVetoes * getVetoes() const
Get the vetoes object.
Int_t nBinsX_
The number of bins on the x-axis of the histogram.
File containing declaration of Lau2DHistDPPdf class.
Double_t binXWidth_
The histogram x-axis bin width.
File containing declaration of LauKinematics class.
TH2 * hist_
The underlying histogram.
void checkNormalisation()
Check the normalisation calculation.
Double_t minX_
The histogram x-axis minimum.
void getUpperHalf(Double_t &x, Double_t &y) const
If only using the upper half of the (symmetric) DP then transform into the correct half...
virtual ~Lau2DHistDPPdf()
Destructor.
Bool_t useInterpolation_
Control boolean for using the linear interpolation.
File containing LauRandom namespace.
Double_t invBinXWidth_
The histogram x-axis inverse bin width.
Bool_t passVeto(Double_t &m12Sq, Double_t &m23Sq, Double_t &m13Sq) const
Check whether the specified Dalitz plot point passes the vetoes.
Definition: LauVetoes.cc:109
Double_t getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
Get the bin content from the histogram.
const LauKinematics * getKinematics() const
Get the kinematics object.
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
void updateKinematics(Double_t x, Double_t y) const
Update the current co-ordinates in the kinematic space.
File containing declaration of LauVetoes class.
Double_t maxY_
The histogram y-axis maximum.
void calcHistNorm()
Calculate the PDF normalisation.
Double_t binYWidth_
The histogram y-axis bin width.
Int_t nBinsY_
The number of bins on the y-axis of the histogram.
Double_t interpolateXY(Double_t x, Double_t y) const
Perform the interpolation (unnormalised)
Class for defining vetoes within the Dalitz plot.
Definition: LauVetoes.hh:33