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