laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
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 using std::cout;
18 using std::cerr;
19 using std::endl;
20 
21 #include "TSystem.h"
22 #include "Lau2DHistDP.hh"
23 #include "Lau2DSplineDP.hh"
24 #include "LauDaughters.hh"
25 #include "LauEffModel.hh"
26 #include "LauKinematics.hh"
27 #include "LauVetoes.hh"
28 
29 ClassImp(LauEffModel)
30 
31 
32 LauEffModel::LauEffModel(const LauDaughters* daughters, const LauVetoes* vetoes) :
33  daughters_( daughters ),
34  vetoes_( vetoes ),
35  effHisto_( 0 ),
36  squareDP_( kFALSE ),
37  fluctuateEffHisto_( kFALSE ),
38  lowBinWarningIssued_( kFALSE ),
39  highBinWarningIssued_( kFALSE )
40 {
41  if ( daughters_ == 0 ) {
42  cerr << "ERROR in LauEffModel Constructor : invalid pointer to daughters object supplied." << endl;
43  gSystem->Exit(EXIT_FAILURE);
44  }
45 }
46 
47 /*LauEffModel::LauEffModel( const LauEffModel& rhs ) :
48  daughters_( rhs.daughters_ ),
49  vetoes_( rhs.vetoes_ ),
50  effHisto_( rhs.effHisto_ ? new Lau2DHistDP( *rhs.effHisto_ ) : 0 ),
51  squareDP_( rhs.squareDP_ ),
52  fluctuateEffHisto_( rhs.fluctuateEffHisto_ )
53 {
54 }*/
55 
57 {
58  if (effHisto_ != 0) {
59  delete effHisto_; effHisto_ = 0;
60  }
61 }
62 
63 void LauEffModel::setEffHisto(const TH2* effHisto, Bool_t useInterpolation,
64  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
65  Bool_t useUpperHalfOnly, Bool_t squareDP)
66 {
67  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
68  // with x = m_13^2, y = m_23^2.
69  Bool_t upperHalf( kFALSE );
70  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
71  cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
72 
73  squareDP_ = squareDP;
74 
75  if(effHisto_) {
76  delete effHisto_;
77  effHisto_=0;
78  }
79 
80  // Copy the histogram.
81  effHisto_ = new Lau2DHistDP(effHisto, daughters_,
82  useInterpolation, fluctuateBins,
83  avEff, absError, upperHalf, squareDP);
84  fluctuateEffHisto_ = fluctuateBins;
85 
86  if (avEff > 0.0 && absError > 0.0) {
87  fluctuateEffHisto_ = kTRUE;
88  }
89 }
90 
91 void LauEffModel::setEffSpline(const TH2* effHisto,
92  Bool_t fluctuateBins, Double_t avEff, Double_t absError,
93  Bool_t useUpperHalfOnly, Bool_t squareDP)
94 {
95  // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
96  // with x = m_13^2, y = m_23^2.
97  Bool_t upperHalf( kFALSE );
98  if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
99  cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
100 
101  squareDP_ = squareDP;
102 
103  if(effHisto_) {
104  delete effHisto_;
105  effHisto_=0;
106  }
107 
108  // Copy the histogram.
109  effHisto_ = new Lau2DSplineDP(effHisto, daughters_,
110  fluctuateBins, avEff, absError, upperHalf, squareDP);
111  fluctuateEffHisto_ = fluctuateBins;
112 
113  if (avEff > 0.0 && absError > 0.0) {
114  fluctuateEffHisto_ = kTRUE;
115  }
116 }
117 
118 Double_t LauEffModel::getEffHistValue(Double_t xVal, Double_t yVal) const
119 {
120  // Get the efficiency from the 2D histo.
121  Double_t eff(1.0);
122 
123  if (effHisto_ != 0) {
124  eff = effHisto_->interpolateXY(xVal, yVal);
125  }
126 
127  return eff;
128 }
129 
130 Double_t LauEffModel::calcEfficiency( const LauKinematics* kinematics ) const
131 {
132  // Routine to calculate the efficiency for the given event/point in
133  // the Dalitz plot. This routine uses the 2D histogram set by the
134  // setEffHisto() funciton.
135  Double_t eff(1.0);
136 
137  // Check for vetoes
138  Bool_t vetoOK(kTRUE);
139  if ( vetoes_ != 0 ) {
140  vetoOK = vetoes_->passVeto(kinematics);
141  }
142 
143  if (vetoOK == kFALSE) {
144  // We failed the veto.
145  eff = 0.0;
146  } else {
147  // We are using a 2D histogram representation of the efficiency across the Dalitz plot.
148  // First find out which bin we are in given the x and y Dalitz plot mass values
149  // Here, x = m13^2, y = m23^2 or if using square DP x = m', y = theta'.
150  if (squareDP_ == kTRUE) {
151  eff = this->getEffHistValue(kinematics->getmPrime(), kinematics->getThetaPrime());
152  } else {
153  eff = this->getEffHistValue(kinematics->getm13Sq(), kinematics->getm23Sq());
154  }
155  }
156 
157  // Check that the efficiency is in the allowed range (0-1)
158  // 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.
159  // The spline requires the efficiency, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
160  // at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
161  // 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
162  // between the two zero bins to remain smooth. The analogous case with adjacent maximised bins will cause peaks above one. Such dips are
163  // unavoidable but are correctly removed here.
164  if ( eff < 0.0 ) {
165  if(!lowBinWarningIssued_) {
166  std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
167  << " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
168  lowBinWarningIssued_=kTRUE;
169  }
170  eff = 0.0;
171  } else if ( eff > 1.0 ) {
172  if(!highBinWarningIssued_) {
173  std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is greater than 1 - setting to 1. You may want to check your histogram!" << std::endl
174  << " : If you are using a spline then this could be caused by adjacent full bins. Further warnings will be suppressed." << std::endl;
175  highBinWarningIssued_=kTRUE;
176  }
177  eff = 1.0;
178  }
179 
180  return eff;
181 }
182 
183 Bool_t LauEffModel::passVeto( const LauKinematics* kinematics ) const
184 {
185  Bool_t pass = kTRUE;
186  if ( vetoes_ != 0 ) {
187  pass = vetoes_->passVeto( kinematics );
188  }
189  return pass;
190 }
191 
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:63
const LauDaughters * daughters_
The daughters object.
Definition: LauEffModel.hh:114
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
File containing declaration of LauDaughters class.
Lau2DAbsDP * effHisto_
The efficiency histogram object.
Definition: LauEffModel.hh:120
virtual Double_t interpolateXY(Double_t x, Double_t y) const =0
Perform the interpolation.
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:118
File containing declaration of LauKinematics class.
Bool_t gotSymmetricalDP() const
Is Dalitz plot symmetric.
Definition: LauDaughters.hh:66
Bool_t passVeto(const LauKinematics *kinematics) const
Determine whether the given DP position is outside the vetoes.
Definition: LauEffModel.cc:183
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:130
Bool_t highBinWarningIssued_
Flag to track whether a warning has been issued for bin values greater than one.
Definition: LauEffModel.hh:131
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:91
const LauVetoes * vetoes_
The vetoes object.
Definition: LauEffModel.hh:117
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:128
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:123
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.
File containing declaration of LauVetoes class.
Bool_t fluctuateEffHisto_
Fluctuate histogram within the error.
Definition: LauEffModel.hh:125
Class for defining vetoes within the Dalitz plot.
Definition: LauVetoes.hh:33
virtual ~LauEffModel()
Destructor.
Definition: LauEffModel.cc:56
Class for defining variations across a 2D DP using a spline.