laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauBkgndDPModel.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 "TRandom.h"
22 #include "TSystem.h"
23 
24 #include "Lau2DHistDPPdf.hh"
25 #include "Lau2DSplineDPPdf.hh"
26 #include "LauBkgndDPModel.hh"
27 #include "LauDaughters.hh"
28 #include "LauFitDataTree.hh"
29 #include "LauKinematics.hh"
30 #include "LauRandom.hh"
31 #include "LauVetoes.hh"
32 
33 ClassImp(LauBkgndDPModel)
34 
35 
37  LauAbsBkgndDPModel(daughters, vetoes),
38  symmetricalDP_(kFALSE),
39  squareDP_(kFALSE),
40  bgHistDPPdf_(0),
41  curEvtHistVal_(0.0),
42  maxPdfHeight_(1.0),
43  pdfNorm_(1.0),
44  doneGenWarning_(kFALSE),
45  lowBinWarningIssued_(kFALSE)
46 {
47  if (daughters != 0) {
48  symmetricalDP_ = daughters->gotSymmetricalDP();
49  }
50 }
51 
53 {
54  if (bgHistDPPdf_ != 0) {
55  delete bgHistDPPdf_; bgHistDPPdf_ = 0;
56  }
57 }
58 
59 void LauBkgndDPModel::setBkgndHisto(const TH2* histo, Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
60 {
61  Bool_t upperHalf = kFALSE;
62  if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
63  cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
64 
65  squareDP_ = squareDP;
66 
67  LauKinematics* kinematics = this->getKinematics();
68  const LauVetoes* vetoes = this->getVetoes();
69  bgHistDPPdf_ = new Lau2DHistDPPdf(histo, kinematics, vetoes, useInterpolation, fluctuateBins, upperHalf, squareDP);
70 
73 }
74 
75 void LauBkgndDPModel::setBkgndSpline(const TH2* histo, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
76 {
77  Bool_t upperHalf = kFALSE;
78  if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
79  cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
80 
81  squareDP_ = squareDP;
82 
83  LauKinematics* kinematics = this->getKinematics();
84  const LauVetoes* vetoes = this->getVetoes();
85  bgHistDPPdf_ = new Lau2DSplineDPPdf(histo, kinematics, vetoes, fluctuateBins, upperHalf, squareDP);
86 
89 }
90 
91 Double_t LauBkgndDPModel::calcHistValue(Double_t xVal, Double_t yVal)
92 {
93  // Get the likelihood value of the background in the Dalitz plot.
94 
95  // Check that we have a valid histogram PDF
96  if (bgHistDPPdf_ == 0) {
97  cerr << "WARNING in LauBkgndDPModel::calcHistValue : We don't have a histogram so assuming the likelihood is flat in the Dalitz plot." << endl;
98  this->setBkgndHisto( 0, kFALSE, kFALSE, kFALSE, kFALSE );
99  }
100 
101  // Find out the un-normalised PDF value
102  Double_t value = bgHistDPPdf_->interpolateXY(xVal, yVal);
103 
104  // Check that the value is greater than zero
105  // If we're using a spline then negative values can be caused by adjacent bins that all contain a value of zero.
106  // The spline requires the value, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
107  // at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
108  // 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
109  // between the two zero bins to remain smooth. Such dips are unavoidable but are correctly removed here.
110  if ( value < 0.0 ) {
111  if(!lowBinWarningIssued_) {
112  std::cerr << "WARNING in LauBkgndDPModel::calcHistValue : Value " << value << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
113  << " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
114  lowBinWarningIssued_=kTRUE;
115  }
116  return 0.0;
117  }
118 
119  LauKinematics* kinematics = this->getKinematics();
120 
121  // For square DP co-ordinates, need to divide by Jacobian
122  if (squareDP_ == kTRUE) {
123  // Make sure that square DP kinematics are correct, then get Jacobian
124  kinematics->updateSqDPKinematics(xVal, yVal);
125  Double_t jacobian = kinematics->calcSqDPJacobian();
126  value /= jacobian;
127  }
128 
129  return value;
130 }
131 
133 {
134 }
135 
137 {
138  // Routine to generate the background, using data provided by an
139  // already defined histogram.
140 
141  LauKinematics* kinematics = this->getKinematics();
142 
143  Bool_t gotBG(kFALSE);
144 
145  while (gotBG == kFALSE) {
146 
147  if (squareDP_ == kTRUE) {
148  // Generate a point in m', theta' space. By construction, this point
149  // is already within the DP region.
150  Double_t mPrime(0.0), thetaPrime(0.0);
151  kinematics->genFlatSqDP(mPrime, thetaPrime);
152 
153  // If we're in a symmetrical DP then we should only generate events in one half
154  if ( symmetricalDP_ && thetaPrime > 0.5 ) {
155  thetaPrime = 1.0 - thetaPrime;
156  }
157 
158  // Calculate histogram height for DP point and
159  // compare with the maximum height
160  if ( bgHistDPPdf_ != 0 ) {
161  Double_t bgContDP = bgHistDPPdf_->interpolateXY(mPrime, thetaPrime)/maxPdfHeight_;
162  if (LauRandom::randomFun()->Rndm() < bgContDP) {
163  kinematics->updateSqDPKinematics(mPrime, thetaPrime);
164  gotBG = kTRUE;
165  }
166  } else {
167  if ( !doneGenWarning_ ) {
168  cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the square DP, which won't be flat in the conventional DP!" << endl;
169  cerr << "WARNING in LauBkgndDPModel::generate : This should never happen!! What have you done?!" << endl;
170  doneGenWarning_ = kTRUE;
171  }
172  kinematics->updateSqDPKinematics(mPrime, thetaPrime);
173  gotBG = kTRUE;
174  }
175  } else {
176  // Generate a point in the Dalitz plot (phase-space).
177  Double_t m13Sq(0.0), m23Sq(0.0);
178  kinematics->genFlatPhaseSpace(m13Sq, m23Sq);
179 
180  // If we're in a symmetrical DP then we should only generate events in one half
181  if ( symmetricalDP_ && m13Sq > m23Sq ) {
182  Double_t tmpSq = m13Sq;
183  m13Sq = m23Sq;
184  m23Sq = tmpSq;
185  }
186 
187  // Calculate histogram height for DP point and
188  // compare with the maximum height
189  if ( bgHistDPPdf_ != 0 ) {
190  Double_t bgContDP = bgHistDPPdf_->interpolateXY(m13Sq, m23Sq)/maxPdfHeight_;
191  if (LauRandom::randomFun()->Rndm() < bgContDP) {
192  kinematics->updateKinematics(m13Sq, m23Sq);
193  gotBG = kTRUE;
194  }
195  } else {
196  if ( !doneGenWarning_ ) {
197  cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the DP." << endl;
198  doneGenWarning_ = kTRUE;
199  }
200  kinematics->updateKinematics(m13Sq, m23Sq);
201  gotBG = kTRUE;
202  }
203  }
204  }
205 
206  // Implement veto
207  Bool_t vetoOK(kTRUE);
208  const LauVetoes* vetoes = this->getVetoes();
209  if (vetoes) {
210  vetoOK = vetoes->passVeto(kinematics);
211  }
212  // Call this function recusively until we pass the veto.
213  if (vetoOK == kFALSE) {
214  this->generate();
215  }
216  return kTRUE;
217 }
218 
220 {
221  // In LauFitDataTree, the first two variables should always be
222  // m13^2 and m23^2. Other variables follow thus: charge.
223 
224  Int_t nBranches = inputFitTree.nBranches();
225  if (nBranches < 2) {
226  cout<<"Error in LauBkgndDPModel::fillDataTree. Expecting at least 2 variables "
227  <<"in input data tree, but there are "<<nBranches<<"! Make sure you have "
228  <<"the right number of variables in your input data file!"<<endl;
229  return;
230  }
231 
232  Double_t m13Sq(0.0), m23Sq(0.0), mPrime(0.0), thetaPrime(0.0);
233 
234  UInt_t nEvents = inputFitTree.nEvents();
235  LauKinematics* kinematics = this->getKinematics();
236 
237  // clear and resize the data vector
238  bgData_.clear();
239  bgData_.resize(nEvents);
240 
241  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
242 
243  const LauFitData& dataValues = inputFitTree.getData(iEvt);
244  LauFitData::const_iterator iter = dataValues.find("m13Sq");
245  m13Sq = iter->second;
246  iter = dataValues.find("m23Sq");
247  m23Sq = iter->second;
248 
249  // Update the kinematics. This will also update m' and theta' if squareDP = true
250  kinematics->updateKinematics(m13Sq, m23Sq);
251 
252  if (squareDP_ == kTRUE) {
253  mPrime = kinematics->getmPrime();
254  thetaPrime = kinematics->getThetaPrime();
255  curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
256  } else {
257  curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
258  }
259 
260  bgData_[iEvt] = curEvtHistVal_;
261  }
262 }
263 
264 Double_t LauBkgndDPModel::getUnNormValue(UInt_t iEvt)
265 {
266  // Retrieve the likelihood for the given event
267  this->setDataEventNo(iEvt);
268  return curEvtHistVal_;
269 }
270 
271 Double_t LauBkgndDPModel::getLikelihood(UInt_t iEvt)
272 {
273  // Retrieve the likelihood for the given event
274  this->setDataEventNo(iEvt);
275  Double_t llhd = curEvtHistVal_ / this->getPdfNorm();
276  return llhd;
277 }
278 
280 {
281  // Retrieve the data for event iEvt
282  if (bgData_.size() > iEvt) {
283  curEvtHistVal_ = bgData_[iEvt];
284  } else {
285  cerr<<"ERROR in LauBkgndDPModel::setDataEventNo : Event index too large: "<<iEvt<<" >= "<<bgData_.size()<<"."<<endl;
286  }
287 }
288 
The abstract interface for a background Dalitz plot model.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
File containing declaration of LauFitDataTree class.
virtual void initialise()
Initialise the model.
Bool_t lowBinWarningIssued_
Flag to track whether a warning has been issued for bin values less than zero.
Double_t maxPdfHeight_
Maximum height of PDF.
Lau2DAbsDPPdf * bgHistDPPdf_
PDF of Dalitz plot background, from a 2D histogram.
void updateSqDPKinematics(Double_t mPrime, Double_t thetaPrime)
Update all kinematic quantities based on the square DP co-ordinates m&#39; and theta&#39;.
File containing declaration of LauBkgndDPModel class.
void setBkgndSpline(const TH2 *histo, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
Set the background histogram and generate a spline.
std::vector< Double_t > bgData_
Cached histogram values for each event.
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
virtual Double_t getHistNorm() const =0
Retrieve PDF normalisation.
File containing declaration of LauDaughters class.
virtual Bool_t generate()
Generate a toy MC event from the model.
Class for defining a 2D DP histogram PDF.
const LauKinematics * getKinematics() const
Get the Dalitz plot kinematics.
std::map< TString, Double_t > LauFitData
Type for holding event data.
UInt_t nBranches() const
Obtain the number of branches in the tree.
Bool_t symmetricalDP_
Flags whether the DP is symmetrical or not.
File containing declaration of Lau2DHistDPPdf class.
virtual Double_t getMaxHeight() const =0
Retrieve maximum height.
Class for defining a 2D DP spline PDF.
File containing declaration of LauKinematics class.
File containing declaration of Lau2DSplineDPPdf class.
virtual Double_t interpolateXY(Double_t x, Double_t y) const =0
Perform the interpolation (unnormalised)
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
Bool_t doneGenWarning_
Boolean to indicate if the warning that there is no histogram has already been issued.
void setBkgndHisto(const TH2 *histo, Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP=kFALSE)
Set background histogram.
virtual Double_t getLikelihood(UInt_t iEvt)
Get likelihood for a given event.
void genFlatSqDP(Double_t &mPrime, Double_t &thetaPrime) const
Routine to generate events flat in the square Dalitz plot.
const LauVetoes * getVetoes() const
Get vetoes in the Dalitz plot.
Double_t calcHistValue(Double_t xVal, Double_t yVal)
Calulate histogram value at a given point.
Double_t curEvtHistVal_
Histogram value for the current event.
Bool_t squareDP_
Flags whether or not to work in square DP coordinates.
Class for defining a histogram-based background Dalitz plot model.
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
File containing LauRandom namespace.
virtual Double_t getUnNormValue(UInt_t iEvt)
Get unnormalised likelihood for a given event.
virtual ~LauBkgndDPModel()
Destructor.
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
Double_t getThetaPrime() const
Get theta&#39; value.
virtual void fillDataTree(const LauFitDataTree &fitDataTree)
Cache the input data and (if appropriate) the per-event likelihood values.
virtual void setDataEventNo(UInt_t iEvt)
Set data event number.
void genFlatPhaseSpace(Double_t &m13Sq, Double_t &m23Sq) const
Routine to generate events flat in phase-space.
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
Double_t pdfNorm_
Normalisation of PDF.
UInt_t nEvents() const
Retrieve the number of events.
File containing declaration of LauVetoes class.
Double_t calcSqDPJacobian()
Calculate the Jacobian for the transformation m23^2, m13^2 -&gt; m&#39;, theta&#39; (square DP) ...
Class for defining vetoes within the Dalitz plot.
Definition: LauVetoes.hh:33
Class to store the input fit variables.
virtual Double_t getPdfNorm() const
Get PDF normalisation constant.