laura is hosted by Hepforge, IPPP Durham
Laura++  v3r5
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DSplineDP.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2013 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include <iostream>
30 
31 #include "TAxis.h"
32 #include "TH2.h"
33 #include "TRandom.h"
34 #include "TSystem.h"
35 
36 #include "Lau2DSplineDP.hh"
37 #include "Lau2DCubicSpline.hh"
38 #include "LauDaughters.hh"
39 #include "LauKinematics.hh"
40 #include "LauRandom.hh"
41 
43 
44 
45 Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const LauDaughters* daughters,
46  Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
47  Bool_t useUpperHalfOnly, Bool_t squareDP) :
48  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
49  spline_(0)
50 {
51  //We may need to modify the histogram so clone it
52  TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
53 
54  if ( ! tempHist ) {
55  std::cerr << "ERROR in Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
56  gSystem->Exit(EXIT_FAILURE);
57  }
58 
59  if (fluctuateBins) {
60  this->doBinFluctuation(tempHist);
61  }
62  if (avEff > 0.0 && avEffError > 0.0) {
63  this->raiseOrLowerBins(tempHist,avEff,avEffError);
64  }
65 
66  spline_ = new Lau2DCubicSpline(*tempHist);
67 
68  delete tempHist;
69 }
70 
71 Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
72  Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
73  Bool_t useUpperHalfOnly, Bool_t squareDP) :
74  Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
75  spline_(0)
76 {
77  //We may need to modify the histogram so clone it
78  TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
79  TH2* tempErrorHi(errorHi ? dynamic_cast<TH2*>(errorHi->Clone()) : 0);
80  TH2* tempErrorLo(errorLo ? dynamic_cast<TH2*>(errorLo->Clone()) : 0);
81 
82  if ( ! tempHist ) {
83  std::cerr << "ERROR in Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
84  gSystem->Exit(EXIT_FAILURE);
85  }
86  if ( ! tempErrorHi ) {
87  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram pointer is null." << std::endl;
88  gSystem->Exit(EXIT_FAILURE);
89  }
90  if ( ! tempErrorLo ) {
91  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram pointer is null." << std::endl;
92  gSystem->Exit(EXIT_FAILURE);
93  }
94 
95  TAxis* xAxis = tempHist->GetXaxis();
96  Double_t minX = static_cast<Double_t>(xAxis->GetXmin());
97  Double_t maxX = static_cast<Double_t>(xAxis->GetXmax());
98 
99  TAxis* yAxis = tempHist->GetYaxis();
100  Double_t minY = static_cast<Double_t>(yAxis->GetXmin());
101  Double_t maxY = static_cast<Double_t>(yAxis->GetXmax());
102 
103  Int_t nBinsX = static_cast<Int_t>(tempHist->GetNbinsX());
104  Int_t nBinsY = static_cast<Int_t>(tempHist->GetNbinsY());
105 
106  if(static_cast<Int_t>(tempErrorLo->GetNbinsX()) != nBinsX ||
107  static_cast<Int_t>(tempErrorLo->GetNbinsY()) != nBinsY) {
108  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different number of bins to the main histogram." << std::endl;
109  gSystem->Exit(EXIT_FAILURE);
110  }
111 
112  if(static_cast<Int_t>(tempErrorHi->GetNbinsX()) != nBinsX ||
113  static_cast<Int_t>(tempErrorHi->GetNbinsY()) != nBinsY) {
114  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different number of bins to the main histogram." << std::endl;
115  gSystem->Exit(EXIT_FAILURE);
116  }
117 
118  xAxis = tempErrorLo->GetXaxis();
119  yAxis = tempErrorLo->GetYaxis();
120 
121  if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
122  static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
123  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different x range to the main histogram." << std::endl;
124  gSystem->Exit(EXIT_FAILURE);
125  }
126 
127  if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
128  static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
129  std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different y range to the main histogram." << std::endl;
130  gSystem->Exit(EXIT_FAILURE);
131  }
132 
133  xAxis = tempErrorHi->GetXaxis();
134  yAxis = tempErrorHi->GetYaxis();
135 
136  if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
137  static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
138  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different x range to the main histogram." << std::endl;
139  gSystem->Exit(EXIT_FAILURE);
140  }
141 
142  if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
143  static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
144  std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different y range to the main histogram." << std::endl;
145  gSystem->Exit(EXIT_FAILURE);
146  }
147 
148 
149  if (fluctuateBins) {
150  this->doBinFluctuation(tempHist,tempErrorHi,tempErrorLo);
151  }
152  if (avEff > 0.0 && avEffError > 0.0) {
153  this->raiseOrLowerBins(tempHist,avEff,avEffError);
154  }
155 
156  spline_ = new Lau2DCubicSpline(*tempHist);
157 
158  delete tempHist;
159  delete tempErrorHi;
160  delete tempErrorLo;
161 }
162 
164 {
165  delete spline_;
166  spline_ = 0;
167 }
168 
169 Double_t Lau2DSplineDP::interpolateXY(Double_t x, Double_t y) const
170 {
171  // This function returns the interpolated value of the histogram function
172  // for the given values of x and y by finding the adjacent bins and extrapolating
173  // using weights based on the inverse distance of the point from the adajcent
174  // bin centres.
175  // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
176 
177  // If we're only using one half then flip co-ordinates
178  // appropriately for conventional or square DP
179  this->getUpperHalf(x,y);
180 
181  // First ask whether the point is inside the kinematic region.
182  if (this->withinDPBoundaries(x,y) == kFALSE) {
183  std::cerr << "WARNING in Lau2DSplineDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
184  return 0.0;
185  }
186 
187  return spline_->evaluate(x,y);
188 
189 }
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:47
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.