laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepMapPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2009 - 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 <iostream>
16 #include <vector>
17 using std::cout;
18 using std::cerr;
19 using std::endl;
20 
21 #include "TH1.h"
22 #include "TMath.h"
23 #include "TSystem.h"
24 
25 #include "LauConstants.hh"
26 #include "LauDaughters.hh"
27 #include "Lau2DHistDP.hh"
28 #include "LauParameter.hh"
29 #include "LauDPDepMapPdf.hh"
30 
31 ClassImp(LauDPDepMapPdf)
32 
33 
34 LauDPDepMapPdf::LauDPDepMapPdf(const std::vector<LauAbsPdf*>& pdfs,
35  const LauDaughters* daughters,
36  const TH2* dpHisto, Bool_t upperHalf) :
37  LauAbsPdf((!pdfs.empty() && pdfs[0]) ? pdfs[0]->varNames() : std::vector<TString>(), std::vector<LauParameter*>(), (!pdfs.empty() && pdfs[0]) ? pdfs[0]->getMinAbscissas() : LauFitData(), (!pdfs.empty() && pdfs[0]) ? pdfs[0]->getMaxAbscissas() : LauFitData()),
38  daughters_( new LauDaughters(*daughters) ),
39  pdfs_(pdfs),
40  dpDependence_( new Lau2DHistDP(dpHisto, daughters, kFALSE, kFALSE, 0, 0, upperHalf, daughters->squareDP()) ),
41  dpAxisDependence_(0),
42  dpAxis_( CentreDist ),
43  indices_(0)
44 {
45  // Constructor for the PDF map.
46  // The index into the PDF collection is read from the DP histogram.
47 
48 
49  // So the first thing we have to do is check the pointers are all valid.
50  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
51 
52  if (!(*iter)) {
53  cerr<<"ERROR in LauDPDepMapPdf constructor: one of the PDF pointers is null."<<endl;
54  gSystem->Exit(EXIT_FAILURE);
55  }
56 
57  // Next check that the abscissa ranges are the same for each PDF
58  if (pdfs_[0]->getMinAbscissa() != (*iter)->getMinAbscissa()) {
59  cerr<<"ERROR in LauDPDepMapPdf constructor: minimum abscissa values not the same for two PDFs."<<endl;
60  gSystem->Exit(EXIT_FAILURE);
61  }
62  if (pdfs_[0]->getMaxAbscissa() != (*iter)->getMaxAbscissa()) {
63  cerr<<"ERROR in LauDPDepMapPdf constructor: maximum abscissa values not the same for two PDFs."<<endl;
64  gSystem->Exit(EXIT_FAILURE);
65  }
66 
67  // Also check that both PDFs are expecting the same number of input variables
68  if (pdfs_[0]->nInputVars() != (*iter)->nInputVars()) {
69  cerr<<"ERROR in LauDPDepMapPdf constructor: number of input variables not the same for two PDFs."<<endl;
70  gSystem->Exit(EXIT_FAILURE);
71  }
72 
73  // Also check that both PDFs are expecting the same variable name(s)
74  if (pdfs_[0]->varNames() != (*iter)->varNames()) {
75  cerr<<"ERROR in LauDPDepMapPdf constructor: variable name(s) not the same for two PDFs."<<endl;
76  gSystem->Exit(EXIT_FAILURE);
77  }
78  }
79 
80  // Then we need to grab all the parameters and pass them to the base class.
81  // This is so that when we are asked for them they can be put into the fit.
82 
83 
84  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
85  std::vector<LauParameter*>& pdfpars = (*iter)->getParameters();
86  this->addParameters(pdfpars);
87 
88  }
89 
90  // Cache the normalisation factor
91  this->calcNorm();
92 }
93 
94 LauDPDepMapPdf::LauDPDepMapPdf(const std::vector<LauAbsPdf*>& pdfs,
95  const LauDaughters* daughters,
96  const TH1* dpAxisHisto,
97  DPAxis dpAxis) :
98  LauAbsPdf((!pdfs.empty() && pdfs[0]) ? pdfs[0]->varNames() : std::vector<TString>(), std::vector<LauParameter*>(), (!pdfs.empty() && pdfs[0]) ? pdfs[0]->getMinAbscissas() : LauFitData(), (!pdfs.empty() && pdfs[0]) ? pdfs[0]->getMaxAbscissas() : LauFitData()),
99  daughters_( new LauDaughters(*daughters) ),
100  pdfs_(pdfs),
101  dpDependence_( 0 ),
102  dpAxisDependence_(dpAxisHisto ? dynamic_cast<TH1*>(dpAxisHisto->Clone()) : 0),
103  dpAxis_( dpAxis ),
104  indices_( 0 )
105 {
106  // Constructor for the PDFs map.
107  // The index into the PDF collection is read from one of
108  // the DP axes.
109 
110  if ( dpAxisDependence_ ) {
111  dpAxisDependence_->SetDirectory(0);
112  }
113 
114  // So the first thing we have to do is check the pointers are all valid.
115  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
116 
117  if (!(*iter)) {
118  cerr<<"ERROR in LauDPDepMapPdf constructor: one of the PDF pointers is null."<<endl;
119  gSystem->Exit(EXIT_FAILURE);
120  }
121 
122  // Next check that the abscissa ranges are the same for each PDF
123  if (pdfs_[0]->getMinAbscissa() != (*iter)->getMinAbscissa()) {
124  cerr<<"ERROR in LauDPDepMapPdf constructor: minimum abscissa values not the same for two PDFs."<<endl;
125  gSystem->Exit(EXIT_FAILURE);
126  }
127  if (pdfs_[0]->getMaxAbscissa() != (*iter)->getMaxAbscissa()) {
128  cerr<<"ERROR in LauDPDepMapPdf constructor: maximum abscissa values not the same for two PDFs."<<endl;
129  gSystem->Exit(EXIT_FAILURE);
130  }
131 
132  // Also check that both PDFs are expecting the same number of input variables
133  if (pdfs_[0]->nInputVars() != (*iter)->nInputVars()) {
134  cerr<<"ERROR in LauDPDepMapPdf constructor: number of input variables not the same for two PDFs."<<endl;
135  gSystem->Exit(EXIT_FAILURE);
136  }
137 
138  // Also check that both PDFs are expecting the same variable name(s)
139  if (pdfs_[0]->varNames() != (*iter)->varNames()) {
140  cerr<<"ERROR in LauDPDepMapPdf constructor: variable name(s) not the same for two PDFs."<<endl;
141  gSystem->Exit(EXIT_FAILURE);
142  }
143  }
144 
145 
146  // Then we need to grab all the parameters and pass them to the base class.
147  // This is so that when we are asked for them they can be put into the fit.
148 
149  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
150  std::vector<LauParameter*>& pdfpars = (*iter)->getParameters();
151  this->addParameters(pdfpars);
152 
153  }
154 
155  // Cache the normalisation factor
156  this->calcNorm();
157 }
158 
160 {
161  // Destructor
162  delete daughters_; daughters_ = 0;
163  delete dpDependence_; dpDependence_ = 0;
165 }
166 
167 LauDPDepMapPdf::LauDPDepMapPdf(const LauDPDepMapPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa()),
168  daughters_( other.daughters_ ),
169  pdfs_( other.pdfs_ ),
170  dpDependence_( (other.dpDependence_) ? new Lau2DHistDP( *other.dpDependence_ ) : 0 ),
171  dpAxisDependence_( other.dpAxisDependence_ ? dynamic_cast<TH1*>( other.dpAxisDependence_->Clone()) : 0),
172  dpAxis_( other.dpAxis_ ),
173  indices_( other.indices_ )
174 {
175  // Copy constructor
176  this->setRandomFun(other.getRandomFun());
177  this->calcNorm();
178 }
179 
180 UInt_t LauDPDepMapPdf::determineDPRegion( Double_t m13Sq, Double_t m23Sq ) const
181 {
182  UInt_t regionIndex(0);
183 
184  LauKinematics* kinematics = daughters_->getKinematics();
185  kinematics->updateKinematics( m13Sq, m23Sq );
186 
187  if ( dpDependence_ ) {
188  if (daughters_->squareDP()){
189  Double_t mprime = kinematics->getmPrime();
190  Double_t thprime = kinematics->getThetaPrime();
191  regionIndex = static_cast<UInt_t>( dpDependence_->interpolateXY( mprime, thprime ) );
192  } else {
193  regionIndex = static_cast<UInt_t>( dpDependence_->interpolateXY( m13Sq, m23Sq ) );
194  }
195  } else if ( dpAxisDependence_ ) {
196  Double_t dpPos(0.0);
197  if ( dpAxis_ == M12 ) {
198  dpPos = kinematics->calcThirdMassSq(m13Sq,m23Sq);
199  } else if ( dpAxis_ == M13 ) {
200  dpPos = m13Sq;
201  } else if ( dpAxis_ == M23 ) {
202  dpPos = m23Sq;
203  } else if ( dpAxis_ == MMIN ) {
204  dpPos = TMath::Min( m13Sq, m23Sq );
205  } else if ( dpAxis_ == MMAX ) {
206  dpPos = TMath::Max( m13Sq, m23Sq );
207  } else {
208  dpPos = kinematics->distanceFromDPCentre(m13Sq,m23Sq);
209  }
210  Int_t bin = dpAxisDependence_->FindFixBin( dpPos );
211  regionIndex = static_cast<UInt_t>( dpAxisDependence_->GetBinContent( bin ) );
212  } else {
213  // This should never happen
214  cerr << "ERROR in LauDPDepMapPdf::determineDPRegion : No means of determining the region!" << endl;
215  gSystem->Exit(EXIT_FAILURE);
216  }
217 
218  return regionIndex;
219 }
220 
222 {
223  // Check that the given abscissa is within the allowed range
224  if (!this->checkRange(abscissas)) {
225  gSystem->Exit(EXIT_FAILURE);
226  }
227 
228  // Create a set of absissas that doesn't contain the DP variables
229  UInt_t nVars = this->nInputVars();
230  LauAbscissas noDPVars( nVars );
231  for ( UInt_t i(0); i < nVars; ++i ) {
232  noDPVars[i] = abscissas[i];
233  }
234 
235  // Find which region of the DP we're in
236  // The DP variables will be abscissas[nInputVars] and
237  // abscissas[nInputVars+1] (if present).
238  UInt_t regionIndex(0);
239  if ( abscissas.size() == nVars+2 ) {
240  Double_t m13Sq = abscissas[nVars];
241  Double_t m23Sq = abscissas[nVars+1];
242  regionIndex = this->determineDPRegion( m13Sq, m23Sq );
243  } else {
244  // This should never happen
245  cerr << "ERROR in LauDPDepMapPdf::calcLikelihoodInfo : DP vars not supplied, no means of determining the region!" << endl;
246  gSystem->Exit(EXIT_FAILURE);
247  }
248 
249  // Check that the region index is valid
250  if ( regionIndex >= pdfs_.size() ) {
251  cerr << "ERROR in LauDPDepMapPdf::calcLikelihoodInfo : No PDF supplied for region " << regionIndex << endl;
252  gSystem->Exit(EXIT_FAILURE);
253  }
254 
255  // Evaluate the normalised PDF values
256  LauAbsPdf* pdf = pdfs_[regionIndex];
257 
258  if ( pdf->isDPDependent() ) {
259  pdf->calcLikelihoodInfo(abscissas);
260  } else {
261  pdf->calcLikelihoodInfo(noDPVars);
262  }
263 
264  Double_t result = pdf->getUnNormLikelihood();
265  this->setUnNormPDFVal(result);
266 
267  Double_t norm = pdf->getNorm();
268  this->setNorm(norm);
269 }
270 
272 {
273  // Nothing to do here, since it is already normalized
274  this->setNorm(1.0);
275 }
276 
278 {
279  // This method gives you the maximum possible height of the PDF.
280  // It combines the maximum heights of the two individual PDFs.
281  // So it would give the true maximum if the two individual maxima coincided.
282  // It is guaranteed to always be >= the true maximum.
283 
284  Double_t m13Sq = kinematics->getm13Sq();
285  Double_t m23Sq = kinematics->getm23Sq();
286 
287  // Find which region of the DP we're in
288  UInt_t regionIndex = this->determineDPRegion( m13Sq, m23Sq );
289 
290  // Check that the region index is valid
291  if ( regionIndex >= pdfs_.size() ) {
292  cerr << "ERROR in LauDPDepMapPdf::calcPDFHeight : No PDF supplied for region " << regionIndex << endl;
293  gSystem->Exit(EXIT_FAILURE);
294  }
295 
296  // Evaluate the normalised PDF values
297  LauAbsPdf* pdf = pdfs_[regionIndex];
298 
299  // Update the heights of the individual PDFs
300  pdf->calcPDFHeight( kinematics );
301 
302  // Find the PDF maximum
303  Double_t height = pdf->getMaxHeight();
304  this->setMaxHeight(height);
305 }
306 
308 {
309  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
310  (*iter)->checkPositiveness();
311  }
312 }
313 
314 // Override the base class methods for cacheInfo and calcLikelihoodInfo(UInt_t iEvt).
315 // For both of these we delegate some functionality to the constituent PDFs.
316 
318 {
319  // delegate to the sub-PDFs to cache their information
320  for ( std::vector<LauAbsPdf*>::const_iterator iter = pdfs_.begin(); iter != pdfs_.end(); ++iter){
321  (*iter)->cacheInfo( inputData );
322  }
323 
324  // now check that the DP variables are included in the data
325  Bool_t hasBranch = inputData.haveBranch( "m13Sq" );
326  hasBranch &= inputData.haveBranch( "m23Sq" );
327  if (!hasBranch) {
328  cerr<<"ERROR in LauDPDepMapPdf::cacheInfo : Input data does not contain Dalitz plot variables m13Sq and/or m23Sq."<<endl;
329  return;
330  }
331 
332  // determine whether we are caching our PDF value
333  Bool_t doCaching( this->nFixedParameters() == this->nParameters() );
334  this->cachePDF( doCaching );
335 
336  if ( !doCaching ) {
337  // in this case we seem to be doing a fit where the parameters are floating
338  // so need to mark that the PDF height is no longer up to date
339  this->heightUpToDate(kFALSE);
340  }
341 
342  // clear the vector and reserve enough space
343  UInt_t nEvents = inputData.nEvents();
344  indices_.clear();
345  indices_.reserve(nEvents);
346 
347  // loop through the events, determine the fraction and store
348  for (UInt_t iEvt = 0; iEvt < nEvents; ++iEvt) {
349 
350  const LauFitData& dataValues = inputData.getData(iEvt);
351 
352  Double_t m13Sq = dataValues.find("m13Sq")->second;
353  Double_t m23Sq = dataValues.find("m23Sq")->second;
354 
355  UInt_t regionIndex = this->determineDPRegion( m13Sq, m23Sq );
356  indices_.push_back( regionIndex );
357  }
358 }
359 
361 {
362  // Get the fraction value for this event
363  UInt_t regionIndex = indices_[iEvt];
364 
365  // Evaluate the normalised PDF value
366  LauAbsPdf* pdf = pdfs_[regionIndex];
367  pdf->calcLikelihoodInfo(iEvt);
368 
369  Double_t result = pdf->getUnNormLikelihood();
370  this->setUnNormPDFVal(result);
371 
372  Double_t norm = pdf->getNorm();
373  this->setNorm(norm);
374 }
375 
Double_t calcThirdMassSq(Double_t firstMassSq, Double_t secondMassSq) const
Calculate the third invariant mass square from the two provided (e.g. mjkSq from mijSq and mikSq) ...
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:454
File containing declaration of LauDPDepMapPdf class.
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:158
LauDaughters * daughters_
Daughter particles.
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:349
virtual void addParameters(std::vector< LauParameter * > &params)
Add parameters to the PDF.
Definition: LauAbsPdf.cc:522
Double_t getmPrime() const
Get m&#39; value.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
virtual UInt_t nParameters() const
Retrieve the number of PDF parameters.
Definition: LauAbsPdf.hh:91
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:278
File containing declaration of LauDaughters class.
virtual void calcPDFHeight(const LauKinematics *kinematics)=0
Calculate the maximum height of the PDF.
UInt_t determineDPRegion(Double_t m13Sq, Double_t m23Sq) const
Determine the DP region.
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
virtual const std::vector< LauParameter * > & getParameters() const
Retrieve the parameters of the PDF, e.g. so that they can be loaded into a fit.
Definition: LauAbsPdf.hh:321
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:410
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:213
TH1 * dpAxisDependence_
1D histogram - DP axis
virtual ~LauDPDepMapPdf()
Destructor.
std::map< TString, Double_t > LauFitData
Type for holding event data.
Double_t getm23Sq() const
Get the m23 invariant mass square.
Bool_t squareDP() const
Determine to use or not the square Dalitz plot.
Definition: LauDaughters.hh:72
std::vector< LauAbsPdf * > pdfs_
The PDFs.
Double_t distanceFromDPCentre() const
Calculate the distance from the currently set (m13Sq, m23Sq) point to the centre of the Dalitz plot (...
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:472
LauDPDepMapPdf(const std::vector< LauAbsPdf * > &pdfs, const LauDaughters *daughters, const TH2 *dpHisto, Bool_t upperHalf=kFALSE)
Constructor - map described by 2D histogram.
virtual Double_t getMaxHeight() const
Retrieve the maximum height.
Definition: LauAbsPdf.hh:303
std::vector< UInt_t > indices_
Cached indices values.
virtual void cacheInfo(const LauFitDataTree &inputData)
Cache information from data.
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
virtual void checkPositiveness()
Check that PDF is positive.
File containing declaration of LauParameter class.
LauKinematics * getKinematics()
Retrieve the Dalitz plot kinematics.
virtual std::vector< TString > varNames() const
Retrieve the names of the abscissas.
Definition: LauAbsPdf.cc:103
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:164
virtual Bool_t isDPDependent() const
Specifies whether or not the PDF is DP dependent.
Definition: LauAbsPdf.hh:110
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:416
File containing declaration of Lau2DHistDP class.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:31
Double_t interpolateXY(Double_t x, Double_t y) const
Perform the interpolation.
Definition: Lau2DHistDP.cc:113
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:355
Double_t getm13Sq() const
Get the m13 invariant mass square.
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)=0
Calculate the likelihood (and all associated information) given value(s) of the abscissa(s) ...
Class for defining a 2D DP histogram.
Definition: Lau2DHistDP.hh:34
File containing LauConstants namespace.
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
Double_t getThetaPrime() const
Get theta&#39; value.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:40
Class for calculating 3-body kinematic quantities.
virtual UInt_t nFixedParameters() const
Retrieve the number of fixed PDF parameters.
Definition: LauAbsPdf.cc:113
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
UInt_t nEvents() const
Retrieve the number of events.
virtual void calcNorm()
Calculate the normalisation.
virtual UInt_t nInputVars() const
Retrieve the number of abscissas.
Definition: LauAbsPdf.hh:103
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:315
virtual Double_t getNorm() const
Retrieve the normalisation factor.
Definition: LauAbsPdf.hh:284
Class to allow having different PDFs in different regions of the DP.
Class to store the input fit variables.
DPAxis dpAxis_
The DP axis we depend on.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:44
Lau2DHistDP * dpDependence_
2D histogram - DP
DPAxis
Define possibilties for the DP axes.