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