laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DSplineDP.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 "Lau2DSplineDP.hh"
23 #include "Lau2DCubicSpline.hh"
24 #include "LauDaughters.hh"
25 #include "LauKinematics.hh"
26 #include "LauRandom.hh"
27 
29 
30 
31 Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const LauDaughters* daughters,
32  Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
33  Bool_t useUpperHalfOnly, Bool_t squareDP) :
34  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
35  spline_(0)
36 {
37  //We may need to modify the histogram so clone it
38  TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
39 
40  if ( ! tempHist ) {
41  std::cerr << "ERROR in Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
42  gSystem->Exit(EXIT_FAILURE);
43  }
44 
45  if (fluctuateBins) {
46  this->doBinFluctuation(tempHist);
47  }
48  if (avEff > 0.0 && avEffError > 0.0) {
49  this->raiseOrLowerBins(tempHist,avEff,avEffError);
50  }
51 
52  spline_ = new Lau2DCubicSpline(*tempHist);
53 
54  delete tempHist;
55 }
56 
57 Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
58  Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
59  Bool_t useUpperHalfOnly, Bool_t squareDP) :
60  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
61  spline_(0)
62 {
63  //We may need to modify the histogram so clone it
64  TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
65  TH2* tempErrorHi(errorHi ? dynamic_cast<TH2*>(errorHi->Clone()) : 0);
66  TH2* tempErrorLo(errorLo ? dynamic_cast<TH2*>(errorLo->Clone()) : 0);
67 
68  if ( ! tempHist ) {
69  std::cerr << "ERROR in Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
70  gSystem->Exit(EXIT_FAILURE);
71  }
72  if ( ! tempErrorHi ) {
73  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram pointer is null." << std::endl;
74  gSystem->Exit(EXIT_FAILURE);
75  }
76  if ( ! tempErrorLo ) {
77  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram pointer is null." << std::endl;
78  gSystem->Exit(EXIT_FAILURE);
79  }
80 
81  TAxis* xAxis = tempHist->GetXaxis();
82  Double_t minX = static_cast<Double_t>(xAxis->GetXmin());
83  Double_t maxX = static_cast<Double_t>(xAxis->GetXmax());
84 
85  TAxis* yAxis = tempHist->GetYaxis();
86  Double_t minY = static_cast<Double_t>(yAxis->GetXmin());
87  Double_t maxY = static_cast<Double_t>(yAxis->GetXmax());
88 
89  Int_t nBinsX = static_cast<Int_t>(tempHist->GetNbinsX());
90  Int_t nBinsY = static_cast<Int_t>(tempHist->GetNbinsY());
91 
92  if(static_cast<Int_t>(tempErrorLo->GetNbinsX()) != nBinsX ||
93  static_cast<Int_t>(tempErrorLo->GetNbinsY()) != nBinsY) {
94  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different number of bins to the main histogram." << std::endl;
95  gSystem->Exit(EXIT_FAILURE);
96  }
97 
98  if(static_cast<Int_t>(tempErrorHi->GetNbinsX()) != nBinsX ||
99  static_cast<Int_t>(tempErrorHi->GetNbinsY()) != nBinsY) {
100  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different number of bins to the main histogram." << std::endl;
101  gSystem->Exit(EXIT_FAILURE);
102  }
103 
104  xAxis = tempErrorLo->GetXaxis();
105  yAxis = tempErrorLo->GetYaxis();
106 
107  if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
108  static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
109  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different x range to the main histogram." << std::endl;
110  gSystem->Exit(EXIT_FAILURE);
111  }
112 
113  if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
114  static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
115  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different y range to the main histogram." << std::endl;
116  gSystem->Exit(EXIT_FAILURE);
117  }
118 
119  xAxis = tempErrorHi->GetXaxis();
120  yAxis = tempErrorHi->GetYaxis();
121 
122  if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
123  static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
124  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different x range to the main histogram." << std::endl;
125  gSystem->Exit(EXIT_FAILURE);
126  }
127 
128  if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
129  static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
130  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different y range to the main histogram." << std::endl;
131  gSystem->Exit(EXIT_FAILURE);
132  }
133 
134 
135  if (fluctuateBins) {
136  this->doBinFluctuation(tempHist,tempErrorHi,tempErrorLo);
137  }
138  if (avEff > 0.0 && avEffError > 0.0) {
139  this->raiseOrLowerBins(tempHist,avEff,avEffError);
140  }
141 
142  spline_ = new Lau2DCubicSpline(*tempHist);
143 
144  delete tempHist;
145  delete tempErrorHi;
146  delete tempErrorLo;
147 }
148 
150 {
151  delete spline_;
152  spline_ = 0;
153 }
154 
155 Double_t Lau2DSplineDP::interpolateXY(Double_t x, Double_t y) const
156 {
157  // This function returns the interpolated value of the histogram function
158  // for the given values of x and y by finding the adjacent bins and extrapolating
159  // using weights based on the inverse distance of the point from the adajcent
160  // bin centres.
161  // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
162 
163  // If we're only using one half then flip co-ordinates
164  // appropriately for conventional or square DP
165  getUpperHalf(x,y);
166 
167  // First ask whether the point is inside the kinematic region.
168  if (withinDPBoundaries(x,y) == kFALSE) {
169  std::cerr << "WARNING in Lau2DSplineDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
170  return 0.0;
171  }
172 
173  return spline_->evaluate(x,y);
174 
175 }
File containing declaration of Lau2DSplineDP class.
Double_t interpolateXY(Double_t x, Double_t y) const
Perform the interpolation.
ClassImp(LauAbsCoeffSet)
virtual ~Lau2DSplineDP()
Destructor.
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
Class for defining a 2D cubic spline based on an input histogram.
File containing declaration of LauDaughters class.
File containing declaration of LauKinematics class.
virtual Double_t evaluate(Double_t x, Double_t y) const
Evaluate the function at given point.
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 LauRandom namespace.
File containing declaration of Lau2DCubicSpline class.
Lau2DSplineDP(const TH2 *hist, const LauDaughters *daughters, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t avEffError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Constructor.
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...
Lau2DCubicSpline * spline_
A 2D cubic spline generated from the histogram.
Class for defining variations across a 2D DP using a spline.