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