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