laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
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 
27 ClassImp(Lau2DHistDP)
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  minX_(0.0), maxX_(0.0),
36  minY_(0.0), maxY_(0.0),
37  rangeX_(0.0), rangeY_(0.0),
38  binXWidth_(0.0), binYWidth_(0.0),
39  nBinsX_(0), nBinsY_(0),
40  useInterpolation_(useInterpolation)
41 {
42  if ( ! hist_ ) {
43  std::cerr << "ERROR in Lau2DHistDP constructor : the histogram pointer is null." << std::endl;
44  gSystem->Exit(EXIT_FAILURE);
45  }
46 
47  // Save various attributes of the histogram
48  // (axis ranges, number of bins, areas)
49  TAxis* xAxis = hist_->GetXaxis();
50  minX_ = static_cast<Double_t>(xAxis->GetXmin());
51  maxX_ = static_cast<Double_t>(xAxis->GetXmax());
52  rangeX_ = maxX_ - minX_;
53 
54  TAxis* yAxis = hist_->GetYaxis();
55  minY_ = static_cast<Double_t>(yAxis->GetXmin());
56  maxY_ = static_cast<Double_t>(yAxis->GetXmax());
57  rangeY_ = maxY_ - minY_;
58 
59  nBinsX_ = static_cast<Int_t>(hist_->GetNbinsX());
60  nBinsY_ = static_cast<Int_t>(hist_->GetNbinsY());
61 
62  binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
63  binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
64 
65  if (fluctuateBins) {
66  this->doBinFluctuation(hist_);
67  }
68  if (avEff > 0.0 && avEffError > 0.0) {
69  this->raiseOrLowerBins(hist_,avEff,avEffError);
70  }
71 }
72 
74 {
75  delete hist_;
76  hist_ = 0;
77 }
78 
79 Double_t Lau2DHistDP::getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
80 {
81  if (xBinNo < 0) {
82  xBinNo = 0;
83  } else if (xBinNo >= nBinsX_) {
84  return 0.0;
85  }
86 
87  if (yBinNo < 0) {
88  yBinNo = 0;
89  } else if (yBinNo >= nBinsY_) {
90  return 0.0;
91  }
92 
93  Double_t value = hist_->GetBinContent(xBinNo+1, yBinNo+1);
94  return value;
95 }
96 
97 Double_t Lau2DHistDP::interpolateXY(Double_t x, Double_t y) const
98 {
99  // This function returns the interpolated value of the histogram function
100  // for the given values of x and y by finding the adjacent bins and extrapolating
101  // using weights based on the inverse distance of the point from the adajcent
102  // bin centres.
103  // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
104 
105  // If we're only using one half then flip co-ordinates
106  // appropriately for conventional or square DP
107  getUpperHalf(x,y);
108 
109  // First ask whether the point is inside the kinematic region.
110  if (withinDPBoundaries(x,y) == kFALSE) {
111  std::cerr << "WARNING in Lau2DHistDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
112  return 0.0;
113  }
114 
115  // Find the 2D histogram bin for x and y
116  Int_t i = Int_t((x - minX_)*nBinsX_/rangeX_);
117  Int_t j = Int_t((y - minY_)*nBinsY_/rangeY_);
118 
119  if (i >= nBinsX_) {i = nBinsX_ - 1;}
120  if (j >= nBinsY_) {j = nBinsY_ - 1;}
121 
122  // Ask whether we want to do the interpolation, since this method is
123  // not reliable for low statistics histograms.
124  if (useInterpolation_ == kFALSE) {return getBinHistValue(i,j);}
125 
126  // Find the bin centres (actual co-ordinate positions, not histogram indices)
127  Double_t cbinx = Double_t(i+0.5)*rangeX_/nBinsX_ + minX_;
128  Double_t cbiny = Double_t(j+0.5)*rangeY_/nBinsY_ + minY_;
129 
130  // If bin centres are outside kinematic region, do not extrapolate
131  if (withinDPBoundaries(cbinx,cbiny) == kFALSE) {return getBinHistValue(i,j);}
132 
133  // Find the adjacent bins
134  Double_t deltax = x - cbinx;
135  Double_t deltay = y - cbiny;
136 
137  Int_t i_adj(0), j_adj(0);
138  if (deltax > 0.0) {
139  i_adj = i + 1;
140  } else {
141  i_adj = i - 1;
142  }
143  if (deltay > 0.0) {
144  j_adj = j + 1;
145  } else {
146  j_adj = j - 1;
147  }
148 
149  Bool_t isXBoundary(kFALSE), isYBoundary(kFALSE);
150 
151  Double_t value(0.0);
152 
153  if (i_adj >= nBinsX_ || i_adj < 0) {isYBoundary = kTRUE;}
154  if (j_adj >= nBinsY_ || j_adj < 0) {isXBoundary = kTRUE;}
155 
156  // In the corners, do no interpolation. Use entry in bin (i,j)
157  if (isXBoundary == kTRUE && isYBoundary == kTRUE) {
158 
159  value = getBinHistValue(i,j);
160 
161  } else if (isXBoundary == kTRUE && isYBoundary == kFALSE) {
162 
163  // Find the adjacent x bin centre
164  Double_t cbinx_adj = Double_t(i_adj+0.5)*rangeX_/nBinsX_ + minX_;
165 
166  if (withinDPBoundaries(cbinx_adj, y) == kFALSE) {
167 
168  // The adjacent bin is outside the DP range. Don't extrapolate.
169  value = getBinHistValue(i,j);
170 
171  } else {
172 
173  Double_t dx0 = TMath::Abs(x - cbinx);
174  Double_t dx1 = TMath::Abs(cbinx_adj - x);
175  Double_t inter_denom = dx0 + dx1;
176 
177  Double_t value1 = getBinHistValue(i,j);
178  Double_t value2 = getBinHistValue(i_adj,j);
179 
180  value = (value1*dx1 + value2*dx0)/inter_denom;
181 
182  }
183 
184  } else if (isYBoundary == kTRUE && isXBoundary == kFALSE) {
185 
186  // Find the adjacent y bin centre
187  Double_t cbiny_adj = Double_t(j_adj+0.5)*rangeY_/nBinsY_ + minY_;
188 
189  if (withinDPBoundaries(x, cbiny_adj) == kFALSE) {
190 
191  // The adjacent bin is outside the DP range. Don't extrapolate.
192  value = getBinHistValue(i,j);
193 
194  } else {
195 
196  Double_t dy0 = TMath::Abs(y - cbiny);
197  Double_t dy1 = TMath::Abs(cbiny_adj - y);
198  Double_t inter_denom = dy0 + dy1;
199 
200  Double_t value1 = getBinHistValue(i,j);
201  Double_t value2 = getBinHistValue(i,j_adj);
202 
203  value = (value1*dy1 + value2*dy0)/inter_denom;
204 
205  }
206 
207  } else {
208 
209  // 2D linear interpolation using inverse distance as weights.
210  // Find the adjacent x and y bin centres
211  Double_t cbinx_adj = Double_t(i_adj+0.5)*rangeX_/nBinsX_ + minX_;
212  Double_t cbiny_adj = Double_t(j_adj+0.5)*rangeY_/nBinsY_ + minY_;
213 
214  if (withinDPBoundaries(cbinx_adj, cbiny_adj) == kFALSE) {
215 
216  // The adjacent bin is outside the DP range. Don't extrapolate.
217  value = getBinHistValue(i,j);
218 
219  } else {
220 
221  Double_t dx0 = TMath::Abs(x - cbinx);
222  Double_t dx1 = TMath::Abs(cbinx_adj - x);
223  Double_t dy0 = TMath::Abs(y - cbiny);
224  Double_t dy1 = TMath::Abs(cbiny_adj - y);
225 
226  Double_t inter_denom = (dx0 + dx1)*(dy0 + dy1);
227 
228  Double_t value1 = getBinHistValue(i,j);
229  Double_t value2 = getBinHistValue(i_adj,j);
230  Double_t value3 = getBinHistValue(i,j_adj);
231  Double_t value4 = getBinHistValue(i_adj,j_adj);
232 
233  value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
234  value /= inter_denom;
235  }
236 
237  }
238 
239  return value;
240 
241 }
242 
Double_t minY_
The histogram y-axis minimum.
Definition: Lau2DHistDP.hh:87
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:103
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:79
File containing declaration of LauKinematics class.
Double_t minX_
The histogram x-axis minimum.
Definition: Lau2DHistDP.hh:83
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.
File containing declaration of Lau2DHistDP class.
Double_t interpolateXY(Double_t x, Double_t y) const
Perform the interpolation.
Definition: Lau2DHistDP.cc:97
File containing LauRandom namespace.
Bool_t useInterpolation_
Control boolean for using the linear interpolation.
Definition: Lau2DHistDP.hh:106
Class for defining a 2D DP histogram.
Definition: Lau2DHistDP.hh:34
virtual ~Lau2DHistDP()
Destructor.
Definition: Lau2DHistDP.cc:73
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:101
TH2 * hist_
The underlying histogram.
Definition: Lau2DHistDP.hh:80
Double_t rangeX_
The histogram x-axis range.
Definition: Lau2DHistDP.hh:91
Double_t rangeY_
The histogram y-axis range.
Definition: Lau2DHistDP.hh:93