laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DHistDP.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 #include "TSystem.h"
21 
22 #include "Lau2DHistDP.hh"
23 #include "LauDaughters.hh"
24 #include "LauKinematics.hh"
25 #include "LauRandom.hh"
26 
28 
29 
30 Lau2DHistDP::Lau2DHistDP(const TH2* hist, const LauDaughters* daughters,
31  Bool_t useInterpolation, Bool_t fluctuateBins,
32  Double_t avEff, Double_t avEffError, Bool_t useUpperHalfOnly, Bool_t squareDP) :
33  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
34  hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
35  errorHi_(0), errorLo_(0),
36  minX_(0.0), maxX_(0.0),
37  minY_(0.0), maxY_(0.0),
38  rangeX_(0.0), rangeY_(0.0),
39  binXWidth_(0.0), binYWidth_(0.0),
40  nBinsX_(0), nBinsY_(0),
41  useInterpolation_(useInterpolation)
42 {
43  if ( ! hist_ ) {
44  std::cerr << "ERROR in Lau2DHistDP constructor : the histogram pointer is null." << std::endl;
45  gSystem->Exit(EXIT_FAILURE);
46  }
47 
48  // Save various attributes of the histogram
49  // (axis ranges, number of bins, areas)
50  TAxis* xAxis = hist_->GetXaxis();
51  minX_ = static_cast<Double_t>(xAxis->GetXmin());
52  maxX_ = static_cast<Double_t>(xAxis->GetXmax());
53  rangeX_ = maxX_ - minX_;
54 
55  TAxis* yAxis = hist_->GetYaxis();
56  minY_ = static_cast<Double_t>(yAxis->GetXmin());
57  maxY_ = static_cast<Double_t>(yAxis->GetXmax());
58  rangeY_ = maxY_ - minY_;
59 
60  nBinsX_ = static_cast<Int_t>(hist_->GetNbinsX());
61  nBinsY_ = static_cast<Int_t>(hist_->GetNbinsY());
62 
63  binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
64  binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
65 
66  if (fluctuateBins) {
67  this->doBinFluctuation(hist_);
68  }
69  if (avEff > 0.0 && avEffError > 0.0) {
70  this->raiseOrLowerBins(hist_,avEff,avEffError);
71  }
72 }
73 
74 Lau2DHistDP::Lau2DHistDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
75  Bool_t useInterpolation, Bool_t fluctuateBins,
76  Double_t avEff, Double_t avEffError, Bool_t useUpperHalfOnly, Bool_t squareDP) :
77  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
78  hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
79  errorHi_(errorHi ? dynamic_cast<TH2*>(errorHi->Clone()) : 0),
80  errorLo_(errorLo ? dynamic_cast<TH2*>(errorLo->Clone()) : 0),
81  minX_(0.0), maxX_(0.0),
82  minY_(0.0), maxY_(0.0),
83  rangeX_(0.0), rangeY_(0.0),
84  binXWidth_(0.0), binYWidth_(0.0),
85  nBinsX_(0), nBinsY_(0),
86  useInterpolation_(useInterpolation)
87 {
88  if ( ! hist_ ) {
89  std::cerr << "ERROR in Lau2DHistDP constructor : the histogram pointer is null." << std::endl;
90  gSystem->Exit(EXIT_FAILURE);
91  }
92  if ( ! errorHi_ ) {
93  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram pointer is null." << std::endl;
94  gSystem->Exit(EXIT_FAILURE);
95  }
96  if ( ! errorLo_ ) {
97  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram pointer is null." << std::endl;
98  gSystem->Exit(EXIT_FAILURE);
99  }
100 
101  // Save various attributes of the histogram
102  // (axis ranges, number of bins, areas)
103  TAxis* xAxis = hist_->GetXaxis();
104  minX_ = static_cast<Double_t>(xAxis->GetXmin());
105  maxX_ = static_cast<Double_t>(xAxis->GetXmax());
106  rangeX_ = maxX_ - minX_;
107 
108  TAxis* yAxis = hist_->GetYaxis();
109  minY_ = static_cast<Double_t>(yAxis->GetXmin());
110  maxY_ = static_cast<Double_t>(yAxis->GetXmax());
111  rangeY_ = maxY_ - minY_;
112 
113  nBinsX_ = static_cast<Int_t>(hist_->GetNbinsX());
114  nBinsY_ = static_cast<Int_t>(hist_->GetNbinsY());
115 
116  binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
117  binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
118 
119  if(static_cast<Int_t>(errorLo_->GetNbinsX()) != nBinsX_ ||
120  static_cast<Int_t>(errorLo_->GetNbinsY()) != nBinsY_) {
121  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different number of bins to the main histogram." << std::endl;
122  gSystem->Exit(EXIT_FAILURE);
123  }
124 
125  if(static_cast<Int_t>(errorHi_->GetNbinsX()) != nBinsX_ ||
126  static_cast<Int_t>(errorHi_->GetNbinsY()) != nBinsY_) {
127  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different number of bins to the main histogram." << std::endl;
128  gSystem->Exit(EXIT_FAILURE);
129  }
130 
131  xAxis = errorLo_->GetXaxis();
132  yAxis = errorLo_->GetYaxis();
133 
134  if(static_cast<Double_t>(xAxis->GetXmin()) != minX_ ||
135  static_cast<Double_t>(xAxis->GetXmax()) != maxX_) {
136  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different x range to the main histogram." << std::endl;
137  gSystem->Exit(EXIT_FAILURE);
138  }
139 
140  if(static_cast<Double_t>(yAxis->GetXmin()) != minY_ ||
141  static_cast<Double_t>(yAxis->GetXmax()) != maxY_) {
142  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different y range to the main histogram." << std::endl;
143  gSystem->Exit(EXIT_FAILURE);
144  }
145 
146  xAxis = errorHi_->GetXaxis();
147  yAxis = errorHi_->GetYaxis();
148 
149  if(static_cast<Double_t>(xAxis->GetXmin()) != minX_ ||
150  static_cast<Double_t>(xAxis->GetXmax()) != maxX_) {
151  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different x range to the main histogram." << std::endl;
152  gSystem->Exit(EXIT_FAILURE);
153  }
154 
155  if(static_cast<Double_t>(yAxis->GetXmin()) != minY_ ||
156  static_cast<Double_t>(yAxis->GetXmax()) != maxY_) {
157  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different y range to the main histogram." << std::endl;
158  gSystem->Exit(EXIT_FAILURE);
159  }
160 
161  if (fluctuateBins) {
163  }
164  if (avEff > 0.0 && avEffError > 0.0) {
165  this->raiseOrLowerBins(hist_,avEff,avEffError);
166  }
167 }
168 
170 {
171  delete hist_;
172  hist_ = 0;
173  if(errorHi_) {
174  delete errorHi_;
175  errorHi_=0;
176  }
177  if(errorLo_) {
178  delete errorLo_;
179  errorLo_=0;
180  }
181 }
182 
183 Double_t Lau2DHistDP::getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
184 {
185  if (xBinNo < 0) {
186  xBinNo = 0;
187  } else if (xBinNo >= nBinsX_) {
188  return 0.0;
189  }
190 
191  if (yBinNo < 0) {
192  yBinNo = 0;
193  } else if (yBinNo >= nBinsY_) {
194  return 0.0;
195  }
196 
197  Double_t value = hist_->GetBinContent(xBinNo+1, yBinNo+1);
198  return value;
199 }
200 
201 Double_t Lau2DHistDP::interpolateXY(Double_t x, Double_t y) const
202 {
203  // This function returns the interpolated value of the histogram function
204  // for the given values of x and y by finding the adjacent bins and extrapolating
205  // using weights based on the inverse distance of the point from the adajcent
206  // bin centres.
207  // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
208 
209  // If we're only using one half then flip co-ordinates
210  // appropriately for conventional or square DP
211  getUpperHalf(x,y);
212 
213  // First ask whether the point is inside the kinematic region.
214  if (withinDPBoundaries(x,y) == kFALSE) {
215  std::cerr << "WARNING in Lau2DHistDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
216  return 0.0;
217  }
218 
219  // Find the 2D histogram bin for x and y
220  Int_t i = Int_t((x - minX_)*nBinsX_/rangeX_);
221  Int_t j = Int_t((y - minY_)*nBinsY_/rangeY_);
222 
223  if (i >= nBinsX_) {i = nBinsX_ - 1;}
224  if (j >= nBinsY_) {j = nBinsY_ - 1;}
225 
226  // Ask whether we want to do the interpolation, since this method is
227  // not reliable for low statistics histograms.
228  if (useInterpolation_ == kFALSE) {return getBinHistValue(i,j);}
229 
230  // Find the bin centres (actual co-ordinate positions, not histogram indices)
231  Double_t cbinx = Double_t(i+0.5)*rangeX_/nBinsX_ + minX_;
232  Double_t cbiny = Double_t(j+0.5)*rangeY_/nBinsY_ + minY_;
233 
234  // If bin centres are outside kinematic region, do not extrapolate
235  if (withinDPBoundaries(cbinx,cbiny) == kFALSE) {return getBinHistValue(i,j);}
236 
237  // Find the adjacent bins
238  Double_t deltax = x - cbinx;
239  Double_t deltay = y - cbiny;
240 
241  Int_t i_adj(0), j_adj(0);
242  if (deltax > 0.0) {
243  i_adj = i + 1;
244  } else {
245  i_adj = i - 1;
246  }
247  if (deltay > 0.0) {
248  j_adj = j + 1;
249  } else {
250  j_adj = j - 1;
251  }
252 
253  Bool_t isXBoundary(kFALSE), isYBoundary(kFALSE);
254 
255  Double_t value(0.0);
256 
257  if (i_adj >= nBinsX_ || i_adj < 0) {isYBoundary = kTRUE;}
258  if (j_adj >= nBinsY_ || j_adj < 0) {isXBoundary = kTRUE;}
259 
260  // In the corners, do no interpolation. Use entry in bin (i,j)
261  if (isXBoundary == kTRUE && isYBoundary == kTRUE) {
262 
263  value = getBinHistValue(i,j);
264 
265  } else if (isXBoundary == kTRUE && isYBoundary == kFALSE) {
266 
267  // Find the adjacent x bin centre
268  Double_t cbinx_adj = Double_t(i_adj+0.5)*rangeX_/nBinsX_ + minX_;
269 
270  if (withinDPBoundaries(cbinx_adj, y) == kFALSE) {
271 
272  // The adjacent bin is outside the DP range. Don't extrapolate.
273  value = getBinHistValue(i,j);
274 
275  } else {
276 
277  Double_t dx0 = TMath::Abs(x - cbinx);
278  Double_t dx1 = TMath::Abs(cbinx_adj - x);
279  Double_t inter_denom = dx0 + dx1;
280 
281  Double_t value1 = getBinHistValue(i,j);
282  Double_t value2 = getBinHistValue(i_adj,j);
283 
284  value = (value1*dx1 + value2*dx0)/inter_denom;
285 
286  }
287 
288  } else if (isYBoundary == kTRUE && isXBoundary == kFALSE) {
289 
290  // Find the adjacent y bin centre
291  Double_t cbiny_adj = Double_t(j_adj+0.5)*rangeY_/nBinsY_ + minY_;
292 
293  if (withinDPBoundaries(x, cbiny_adj) == kFALSE) {
294 
295  // The adjacent bin is outside the DP range. Don't extrapolate.
296  value = getBinHistValue(i,j);
297 
298  } else {
299 
300  Double_t dy0 = TMath::Abs(y - cbiny);
301  Double_t dy1 = TMath::Abs(cbiny_adj - y);
302  Double_t inter_denom = dy0 + dy1;
303 
304  Double_t value1 = getBinHistValue(i,j);
305  Double_t value2 = getBinHistValue(i,j_adj);
306 
307  value = (value1*dy1 + value2*dy0)/inter_denom;
308 
309  }
310 
311  } else {
312 
313  // 2D linear interpolation using inverse distance as weights.
314  // Find the adjacent x and y bin centres
315  Double_t cbinx_adj = Double_t(i_adj+0.5)*rangeX_/nBinsX_ + minX_;
316  Double_t cbiny_adj = Double_t(j_adj+0.5)*rangeY_/nBinsY_ + minY_;
317 
318  if (withinDPBoundaries(cbinx_adj, cbiny_adj) == kFALSE) {
319 
320  // The adjacent bin is outside the DP range. Don't extrapolate.
321  value = getBinHistValue(i,j);
322 
323  } else {
324 
325  Double_t dx0 = TMath::Abs(x - cbinx);
326  Double_t dx1 = TMath::Abs(cbinx_adj - x);
327  Double_t dy0 = TMath::Abs(y - cbiny);
328  Double_t dy1 = TMath::Abs(cbiny_adj - y);
329 
330  Double_t inter_denom = (dx0 + dx1)*(dy0 + dy1);
331 
332  Double_t value1 = getBinHistValue(i,j);
333  Double_t value2 = getBinHistValue(i_adj,j);
334  Double_t value3 = getBinHistValue(i,j_adj);
335  Double_t value4 = getBinHistValue(i_adj,j_adj);
336 
337  value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
338  value /= inter_denom;
339  }
340 
341  }
342 
343  return value;
344 
345 }
346 
Double_t minY_
The histogram y-axis minimum.
Definition: Lau2DHistDP.hh:112
ClassImp(LauAbsCoeffSet)
void raiseOrLowerBins(TH2 *hist, const Double_t avEff, const Double_t avEffError)
Rescale the histogram bin contents based on the desired average efficiency and its uncertainty...
void doBinFluctuation(TH2 *hist)
Fluctuate the contents of each histogram bin independently, in accordance with their errors...
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
Int_t nBinsY_
The number of bins on the y-axis of the histogram.
Definition: Lau2DHistDP.hh:128
Double_t binYWidth_
The histogram y-axis bin width.
Definition: Lau2DHistDP.hh:123
File containing declaration of LauDaughters class.
Double_t getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
Get the raw bin content from the histogram.
Definition: Lau2DHistDP.cc:183
File containing declaration of LauKinematics class.
Double_t minX_
The histogram x-axis minimum.
Definition: Lau2DHistDP.hh:108
Double_t binXWidth_
The histogram x-axis bin width.
Definition: Lau2DHistDP.hh:121
Abstract base class for defining a variation across a 2D DP based on a histogram. ...
Lau2DHistDP(const TH2 *hist, const LauDaughters *daughters, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t avEffError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Constructor.
Definition: Lau2DHistDP.cc:30
Bool_t withinDPBoundaries(Double_t x, Double_t y) const
Check whether the given co-ordinates are within the kinematic boundary.
File containing declaration of Lau2DHistDP class.
Double_t interpolateXY(Double_t x, Double_t y) const
Perform the interpolation.
Definition: Lau2DHistDP.cc:201
File containing LauRandom namespace.
Bool_t useInterpolation_
Control boolean for using the linear interpolation.
Definition: Lau2DHistDP.hh:131
Double_t maxX_
The histogram x-axis maximum.
Definition: Lau2DHistDP.hh:110
Double_t maxY_
The histogram y-axis maximum.
Definition: Lau2DHistDP.hh:114
Class for defining a 2D DP histogram.
Definition: Lau2DHistDP.hh:34
virtual ~Lau2DHistDP()
Destructor.
Definition: Lau2DHistDP.cc:169
TH2 * errorLo_
The histogram containing the lower errors.
Definition: Lau2DHistDP.hh:105
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...
Double_t value() const
The value of the parameter.
Int_t nBinsX_
The number of bins on the x-axis of the histogram.
Definition: Lau2DHistDP.hh:126
TH2 * errorHi_
The histogram containing the upper errors.
Definition: Lau2DHistDP.hh:103
TH2 * hist_
The underlying histogram.
Definition: Lau2DHistDP.hh:101
Double_t rangeX_
The histogram x-axis range.
Definition: Lau2DHistDP.hh:116
Double_t rangeY_
The histogram y-axis range.
Definition: Lau2DHistDP.hh:118