laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauEffModel.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 <cstdlib>
16 #include <iostream>
17 
18 #include "TSystem.h"
19 #include "Lau2DHistDP.hh"
20 #include "Lau2DSplineDP.hh"
21 #include "LauDaughters.hh"
22 #include "LauEffModel.hh"
23 #include "LauKinematics.hh"
24 #include "LauVetoes.hh"
25 
27 
28 
29 LauEffModel::LauEffModel(const LauDaughters* daughters, const LauVetoes* vetoes) :
30  daughters_( daughters ),
31  vetoes_( vetoes ),
32  effHisto_( 0 ),
33  squareDP_( kFALSE ),
34  fluctuateEffHisto_( kFALSE ),
35  lowBinWarningIssued_( kFALSE ),
36  highBinWarningIssued_( kFALSE )
37 {
38  if ( daughters_ == 0 ) {
39  std::cerr << "ERROR in LauEffModel Constructor : invalid pointer to daughters object supplied." << std::endl;
40  gSystem->Exit(EXIT_FAILURE);
41  }
42 }
43 
45 {
46  std::vector<Lau2DAbsDP*>::iterator it = effHisto_.begin();
47  std::vector<Lau2DAbsDP*>::iterator end = effHisto_.end();
48  for( ; it!=end; ++it) {
49  delete *it;
50  *it=0;
51  }
52  effHisto_.clear();
53 }
54 
55 void LauEffModel::setEffHisto(const TH2* effHisto, Bool_t useInterpolation,
56  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
57  Bool_t useUpperHalfOnly, Bool_t squareDP)
58 {
59  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
60  // with x = m_13^2, y = m_23^2.
61  Bool_t upperHalf( kFALSE );
62  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
63  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
64 
65  squareDP_ = squareDP;
66 
67  std::vector<Lau2DAbsDP*>::iterator it = effHisto_.begin();
68  std::vector<Lau2DAbsDP*>::iterator end = effHisto_.end();
69  for( ; it!=end; ++it) {
70  delete *it;
71  *it=0;
72  }
73  effHisto_.clear();
74 
75  // Copy the histogram.
76  effHisto_.push_back(new Lau2DHistDP(effHisto, daughters_,
77  useInterpolation, fluctuateBins,
78  avEff, absError, upperHalf, squareDP));
79  fluctuateEffHisto_ = fluctuateBins;
80 
81  if (avEff > 0.0 && absError > 0.0) {
82  fluctuateEffHisto_ = kTRUE;
83  }
84 }
85 
86 void LauEffModel::setEffHisto(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo, Bool_t useInterpolation,
87  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
88  Bool_t useUpperHalfOnly, Bool_t squareDP)
89 {
90  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
91  // with x = m_13^2, y = m_23^2.
92  Bool_t upperHalf( kFALSE );
93  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
94  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
95 
96  squareDP_ = squareDP;
97 
98  std::vector<Lau2DAbsDP*>::iterator it = effHisto_.begin();
99  std::vector<Lau2DAbsDP*>::iterator end = effHisto_.end();
100  for( ; it!=end; ++it) {
101  delete *it;
102  *it=0;
103  }
104  effHisto_.clear();
105 
106  // Copy the histogram.
107  effHisto_.push_back(new Lau2DHistDP(effHisto, errorHi, errorLo, daughters_,
108  useInterpolation, fluctuateBins,
109  avEff, absError, upperHalf, squareDP));
110  fluctuateEffHisto_ = fluctuateBins;
111 
112  if (avEff > 0.0 && absError > 0.0) {
113  fluctuateEffHisto_ = kTRUE;
114  }
115 }
116 
117 void LauEffModel::setEffSpline(const TH2* effHisto,
118  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
119  Bool_t useUpperHalfOnly, Bool_t squareDP)
120 {
121  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
122  // with x = m_13^2, y = m_23^2.
123  Bool_t upperHalf( kFALSE );
124  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
125  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
126 
127  squareDP_ = squareDP;
128 
129  std::vector<Lau2DAbsDP*>::iterator it = effHisto_.begin();
130  std::vector<Lau2DAbsDP*>::iterator end = effHisto_.end();
131  for( ; it!=end; ++it) {
132  delete *it;
133  *it=0;
134  }
135  effHisto_.clear();
136 
137  // Copy the histogram.
138  effHisto_.push_back(new Lau2DSplineDP(effHisto, daughters_,
139  fluctuateBins, avEff, absError, upperHalf, squareDP));
140  fluctuateEffHisto_ = fluctuateBins;
141 
142  if (avEff > 0.0 && absError > 0.0) {
143  fluctuateEffHisto_ = kTRUE;
144  }
145 }
146 
147 void LauEffModel::setEffSpline(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo,
148  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
149  Bool_t useUpperHalfOnly, Bool_t squareDP)
150 {
151  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
152  // with x = m_13^2, y = m_23^2.
153  Bool_t upperHalf( kFALSE );
154  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
155  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
156 
157  squareDP_ = squareDP;
158 
159  std::vector<Lau2DAbsDP*>::iterator it = effHisto_.begin();
160  std::vector<Lau2DAbsDP*>::iterator end = effHisto_.end();
161  for( ; it!=end; ++it) {
162  delete *it;
163  *it=0;
164  }
165  effHisto_.clear();
166 
167  // Copy the histogram.
168  effHisto_.push_back(new Lau2DSplineDP(effHisto, errorHi, errorLo, daughters_,
169  fluctuateBins, avEff, absError, upperHalf, squareDP));
170  fluctuateEffHisto_ = fluctuateBins;
171 
172  if (avEff > 0.0 && absError > 0.0) {
173  fluctuateEffHisto_ = kTRUE;
174  }
175 }
176 
177 void LauEffModel::addEffHisto(const TH2* effHisto, Bool_t useInterpolation,
178  Double_t avEff, Double_t absError,
179  Bool_t useUpperHalfOnly, Bool_t squareDP)
180 {
181  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
182  // with x = m_13^2, y = m_23^2.
183  Bool_t upperHalf( kFALSE );
184  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
185  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
186 
187  squareDP_ = squareDP;
188 
189  // Copy the histogram.
190  effHisto_.push_back(new Lau2DHistDP(effHisto, daughters_,
191  useInterpolation, fluctuateEffHisto_,
192  avEff, absError, upperHalf, squareDP));
193 }
194 
195 void LauEffModel::addEffHisto(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo, Bool_t useInterpolation,
196  Double_t avEff, Double_t absError,
197  Bool_t useUpperHalfOnly, Bool_t squareDP)
198 {
199  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
200  // with x = m_13^2, y = m_23^2.
201  Bool_t upperHalf( kFALSE );
202  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
203  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
204 
205  squareDP_ = squareDP;
206 
207  // Copy the histogram.
208  effHisto_.push_back(new Lau2DHistDP(effHisto, errorHi, errorLo, daughters_,
209  useInterpolation, fluctuateEffHisto_,
210  avEff, absError, upperHalf, squareDP));
211 }
212 
213 void LauEffModel::addEffSpline(const TH2* effHisto,
214  Double_t avEff, Double_t absError,
215  Bool_t useUpperHalfOnly, Bool_t squareDP)
216 {
217  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
218  // with x = m_13^2, y = m_23^2.
219  Bool_t upperHalf( kFALSE );
220  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
221  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
222 
223  squareDP_ = squareDP;
224 
225  // Copy the histogram.
226  effHisto_.push_back(new Lau2DSplineDP(effHisto, daughters_,
227  fluctuateEffHisto_, avEff, absError, upperHalf, squareDP));
228 }
229 
230 void LauEffModel::addEffSpline(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo,
231  Double_t avEff, Double_t absError,
232  Bool_t useUpperHalfOnly, Bool_t squareDP)
233 {
234  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
235  // with x = m_13^2, y = m_23^2.
236  Bool_t upperHalf( kFALSE );
237  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
238  std::cout<<"INFO in LauEffModel::setEffSpline : Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<std::endl;
239 
240  squareDP_ = squareDP;
241 
242  // Copy the histogram.
243  effHisto_.push_back(new Lau2DSplineDP(effHisto, errorHi, errorLo, daughters_,
244  fluctuateEffHisto_, avEff, absError, upperHalf, squareDP));
245 }
246 
247 Double_t LauEffModel::getEffHistValue(Double_t xVal, Double_t yVal) const
248 {
249  // Get the efficiency from the 2D histo.
250  Double_t eff(1.0);
251 
252  std::vector<Lau2DAbsDP*>::const_iterator it = effHisto_.begin();
253  std::vector<Lau2DAbsDP*>::const_iterator end = effHisto_.end();
254  for( ; it!=end; ++it) {
255  eff *= (*it)->interpolateXY(xVal, yVal);
256  }
257 
258  return eff;
259 }
260 
261 Double_t LauEffModel::calcEfficiency( const LauKinematics* kinematics ) const
262 {
263  // Routine to calculate the efficiency for the given event/point in
264  // the Dalitz plot. This routine uses the 2D histogram set by the
265  // setEffHisto() function.
266  Double_t eff(1.0);
267 
268  // Check for vetoes
269  Bool_t vetoOK(kTRUE);
270  if ( vetoes_ != 0 ) {
271  vetoOK = vetoes_->passVeto(kinematics);
272  }
273 
274  if (vetoOK == kFALSE) {
275  // We failed the veto.
276  eff = 0.0;
277  } else {
278  // We are using a 2D histogram representation of the efficiency across the Dalitz plot.
279  // First find out which bin we are in given the x and y Dalitz plot mass values
280  // Here, x = m13^2, y = m23^2 or if using square DP x = m', y = theta'.
281  if (squareDP_ == kTRUE) {
282  eff = this->getEffHistValue(kinematics->getmPrime(), kinematics->getThetaPrime());
283  } else {
284  eff = this->getEffHistValue(kinematics->getm13Sq(), kinematics->getm23Sq());
285  }
286  }
287 
288  // Check that the efficiency is in the allowed range (0-1)
289  // If we're using a spline then out-of-range efficiencies can be caused by adjacent bins that all contain a value of either zero or one.
290  // The spline requires the efficiency, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
291  // at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
292  // If two bins are zero but the third is not then the second bin will have a positive first derivative causing the spline to dip below zero
293  // between the two zero bins to remain smooth. The analogous case with adjacent maximised bins will cause peaks above one. Such dips are
294  // unavoidable but are correctly removed here.
295  if ( eff < 0.0 ) {
296  if(!lowBinWarningIssued_) {
297  std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
298  << " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
299  lowBinWarningIssued_=kTRUE;
300  }
301  eff = 0.0;
302  } else if ( eff > 1.0 ) {
303  if(!highBinWarningIssued_) {
304  std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is greater than 1 - setting to 1. You may want to check your histogram!" << std::endl
305  << " : If you are using a spline then this could be caused by adjacent full bins. Further warnings will be suppressed." << std::endl;
306  highBinWarningIssued_=kTRUE;
307  }
308  eff = 1.0;
309  }
310 
311  return eff;
312 }
313 
314 Bool_t LauEffModel::passVeto( const LauKinematics* kinematics ) const
315 {
316  Bool_t pass = kTRUE;
317  if ( vetoes_ != 0 ) {
318  pass = vetoes_->passVeto( kinematics );
319  }
320  return pass;
321 }
322 
File containing declaration of Lau2DSplineDP class.
void setEffHisto(const TH2 *effHisto, Bool_t useInterpolation=kTRUE, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Set the efficiency variation across the phase space using a predetermined 2D histogram.
Definition: LauEffModel.cc:55
ClassImp(LauAbsCoeffSet)
const LauDaughters * daughters_
The daughters object.
Definition: LauEffModel.hh:242
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
void addEffHisto(const TH2 *effHisto, Bool_t useInterpolation=kTRUE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Add a multiplicative efficiency variation across the phase space using a predetermined 2D histogram...
Definition: LauEffModel.cc:177
File containing declaration of LauDaughters class.
Double_t getm23Sq() const
Get the m23 invariant mass square.
Double_t getEffHistValue(Double_t xVal, Double_t yVal) const
Get the efficiency from a two-dimensional histogram by interpolating in x and y.
Definition: LauEffModel.cc:247
File containing declaration of LauKinematics class.
Bool_t gotSymmetricalDP() const
Is Dalitz plot symmetric, i.e. 2 identical particles.
Definition: LauDaughters.hh:66
Bool_t passVeto(const LauKinematics *kinematics) const
Determine whether the given DP position is outside the vetoes.
Definition: LauEffModel.cc:314
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:261
Bool_t highBinWarningIssued_
Flag to track whether a warning has been issued for bin values greater than one.
Definition: LauEffModel.hh:259
void setEffSpline(const TH2 *effHisto, Bool_t fluctuateBins=kFALSE, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Set the efficiency variation across the phase space using a spline based on a predetermined 2D histog...
Definition: LauEffModel.cc:117
const LauVetoes * vetoes_
The vetoes object.
Definition: LauEffModel.hh:245
File containing declaration of Lau2DHistDP class.
Bool_t lowBinWarningIssued_
Flag to track whether a warning has been issued for bin values less than zero.
Definition: LauEffModel.hh:256
File containing declaration of LauEffModel class.
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
Bool_t squareDP_
Use of the square Dalitz plot.
Definition: LauEffModel.hh:251
Double_t getm13Sq() const
Get the m13 invariant mass square.
Bool_t passVeto(Double_t &m12Sq, Double_t &m23Sq, Double_t &m13Sq) const
Check whether the specified Dalitz plot point passes the vetoes.
Definition: LauVetoes.cc:109
Class for defining a 2D DP histogram.
Definition: Lau2DHistDP.hh:34
Double_t getThetaPrime() const
Get theta&#39; value.
Class for calculating 3-body kinematic quantities.
std::vector< Lau2DAbsDP * > effHisto_
The efficiency histogram objects.
Definition: LauEffModel.hh:248
File containing declaration of LauVetoes class.
Bool_t fluctuateEffHisto_
Fluctuate histogram within the error.
Definition: LauEffModel.hh:253
void addEffSpline(const TH2 *effHisto, Double_t avEff=-1.0, Double_t absError=-1.0, Bool_t useUpperHalfOnly=kFALSE, Bool_t squareDP=kFALSE)
Add a multiplicative efficiency variation across the phase space using a spline based on a predetermi...
Definition: LauEffModel.cc:213
Class for defining vetoes within the Dalitz plot.
Definition: LauVetoes.hh:33
virtual ~LauEffModel()
Destructor.
Definition: LauEffModel.cc:44
Class for defining variations across a 2D DP using a spline.