laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepSumPdf.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 using std::vector;
21 
22 #include "TMath.h"
23 #include "TSystem.h"
24 
25 #include "LauConstants.hh"
26 #include "LauDaughters.hh"
27 #include "LauEffModel.hh"
28 #include "LauParameter.hh"
29 #include "LauDPDepSumPdf.hh"
30 
32 
33 
35  const LauDaughters* daughters,
36  const TH2* dpHisto, Bool_t upperHalf, Bool_t useSpline) :
37  LauAbsPdf(pdf1 ? pdf1->varNames() : vector<TString>(), vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
38  daughters_( new LauDaughters(*daughters) ),
39  pdf1_(pdf1),
40  pdf2_(pdf2),
41  frac_(0),
42  fracVal_(0.5),
43  dpDependence_( new LauEffModel(daughters, 0) ),
44  dpAxis_( CentreDist )
45 {
46  // Constructor for the sum PDF.
47  // We are defining the sum as:
48  // f x (PDF1/S(PDF1)) + (1-f) x (PDF2/S(PDF2))
49  // where f is the fraction, x is multiplication, PDFi is the i'th PDF,
50  // and S(PDFi) is the integral of the i'th PDF.
51  // The value of the fraction is read from the DP histogram.
52 
53  if(useSpline) {
54  dpDependence_->setEffSpline( dpHisto, kFALSE, 0.0, 0.0, upperHalf, daughters->squareDP());
55  } else {
56  dpDependence_->setEffHisto( dpHisto, kTRUE, kFALSE, 0.0, 0.0, upperHalf, daughters->squareDP() );
57  }
58 
59  // So the first thing we have to do is check the pointers are all valid.
60  if (!pdf1 || !pdf2) {
61  cerr<<"ERROR in LauDPDepSumPdf constructor: one of the 2 PDF pointers is null."<<endl;
62  gSystem->Exit(EXIT_FAILURE);
63  }
64 
65  // Next check that the abscissa ranges are the same for each PDF
66  if (pdf1->getMinAbscissa() != pdf2->getMinAbscissa()) {
67  cerr<<"ERROR in LauDPDepSumPdf constructor: minimum abscissa values not the same for the two PDFs."<<endl;
68  gSystem->Exit(EXIT_FAILURE);
69  }
70  if (pdf1->getMaxAbscissa() != pdf2->getMaxAbscissa()) {
71  cerr<<"ERROR in LauDPDepSumPdf constructor: maximum abscissa values not the same for the two PDFs."<<endl;
72  gSystem->Exit(EXIT_FAILURE);
73  }
74 
75  // Also check that both PDFs are expecting the same number of input variables
76  if (pdf1->nInputVars() != pdf2->nInputVars()) {
77  cerr<<"ERROR in LauDPDepSumPdf constructor: number of input variables not the same for the two PDFs."<<endl;
78  gSystem->Exit(EXIT_FAILURE);
79  }
80 
81  // Also check that both PDFs are expecting the same variable name(s)
82  if (pdf1->varNames() != pdf2->varNames()) {
83  cerr<<"ERROR in LauDPDepSumPdf constructor: variable name(s) not the same for the two PDFs."<<endl;
84  gSystem->Exit(EXIT_FAILURE);
85  }
86 
87  // Then we need to grab all the parameters and pass them to the base class.
88  // This is so that when we are asked for them they can be put into the fit.
89  // The number of parameters is the number in PDF1 + the number in PDF2.
90  UInt_t nPar(pdf1->nParameters()+pdf2->nParameters());
91  vector<LauAbsRValue*> params; params.reserve(nPar);
92  vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
93  vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
94  for (vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
95  params.push_back(*iter);
96  }
97  for (vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
98  params.push_back(*iter);
99  }
100  this->addParameters(params);
101 
102  // Cache the normalisation factor
103  this->calcNorm();
104 }
105 
107  LauParameter* frac,
108  const LauDaughters* daughters,
109  const std::vector<Double_t>& fracCoeffs,
110  DPAxis dpAxis) :
111  LauAbsPdf(pdf1 ? pdf1->varNames() : vector<TString>(), vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
112  daughters_( new LauDaughters(*daughters) ),
113  pdf1_(pdf1),
114  pdf2_(pdf2),
115  frac_(frac),
116  fracVal_( frac ? frac->value() : 0.0 ),
117  dpDependence_( 0 ),
118  fracCoeffs_( fracCoeffs ),
119  dpAxis_( dpAxis )
120 {
121  // Constructor for the sum PDF.
122  // We are defining the sum as:
123  // f x (PDF1/S(PDF1)) + (1-f) x (PDF2/S(PDF2))
124  // where f is the fraction, x is multiplication, PDFi is the i'th PDF,
125  // and S(PDFi) is the integral of the i'th PDF.
126  // The value of the fraction has a polynomial dependence on one of
127  // the DP axes.
128 
129  // So the first thing we have to do is check the pointers are all valid.
130  if (!pdf1 || !pdf2) {
131  cerr<<"ERROR in LauDPDepSumPdf constructor: one of the 2 PDF pointers is null."<<endl;
132  gSystem->Exit(EXIT_FAILURE);
133  }
134  if ( !frac ) {
135  cerr<<"ERROR in LauDPDepSumPdf constructor: the fraction parameter pointer is null."<<endl;
136  gSystem->Exit(EXIT_FAILURE);
137  }
138 
139  // Next check that the abscissa ranges are the same for each PDF
140  if (pdf1->getMinAbscissa() != pdf2->getMinAbscissa()) {
141  cerr<<"ERROR in LauDPDepSumPdf constructor: minimum abscissa values not the same for the two PDFs."<<endl;
142  gSystem->Exit(EXIT_FAILURE);
143  }
144  if (pdf1->getMaxAbscissa() != pdf2->getMaxAbscissa()) {
145  cerr<<"ERROR in LauDPDepSumPdf constructor: maximum abscissa values not the same for the two PDFs."<<endl;
146  gSystem->Exit(EXIT_FAILURE);
147  }
148 
149  // Also check that both PDFs are expecting the same number of input variables
150  if (pdf1->nInputVars() != pdf2->nInputVars()) {
151  cerr<<"ERROR in LauDPDepSumPdf constructor: number of input variables not the same for the two PDFs."<<endl;
152  gSystem->Exit(EXIT_FAILURE);
153  }
154 
155  // Also check that both PDFs are expecting the same variable name(s)
156  if (pdf1->varNames() != pdf2->varNames()) {
157  cerr<<"ERROR in LauDPDepSumPdf constructor: variable name(s) not the same for the two PDFs."<<endl;
158  gSystem->Exit(EXIT_FAILURE);
159  }
160 
161  // Then we need to grab all the parameters and pass them to the base class.
162  // This is so that when we are asked for them they can be put into the fit.
163  // The number of parameters is the number in PDF1 + the number in PDF2.
164  UInt_t nPar( pdf1->nParameters() + pdf2->nParameters() + 1 );
165  vector<LauAbsRValue*> params; params.reserve(nPar);
166  params.push_back(frac);
167  vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
168  vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
169  for (vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
170  params.push_back(*iter);
171  }
172  for (vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
173  params.push_back(*iter);
174  }
175  this->addParameters(params);
176 
177  // Now check that we can find the fraction parameter ok
178  frac_ = this->findParameter("frac");
179  if (frac_ == 0) {
180  cerr<<"ERROR in LauDPDepSumPdf constructor: parameter \"frac\" not found."<<endl;
181  gSystem->Exit(EXIT_FAILURE);
182  }
183 
184  // Cache the normalisation factor
185  this->calcNorm();
186 }
187 
189 {
190  // Destructor
191  delete daughters_; daughters_ = 0;
192  delete dpDependence_; dpDependence_ = 0;
193 }
194 
196 {
197  // Check that the given abscissa is within the allowed range
198  if (!this->checkRange(abscissas)) {
199  gSystem->Exit(EXIT_FAILURE);
200  }
201 
202  LauAbscissas noDPVars(1);
203  noDPVars[0] = abscissas[0];
204 
205  // Evaluate the normalised PDF values
206  if ( pdf1_->isDPDependent() ) {
207  pdf1_->calcLikelihoodInfo(abscissas);
208  } else {
209  pdf1_->calcLikelihoodInfo(noDPVars);
210  }
211  if ( pdf2_->isDPDependent() ) {
212  pdf2_->calcLikelihoodInfo(abscissas);
213  } else {
214  pdf2_->calcLikelihoodInfo(noDPVars);
215  }
216  Double_t result1 = pdf1_->getLikelihood();
217  Double_t result2 = pdf2_->getLikelihood();
218 
219  // Determine the fraction
220  // The DP variables will be abscissas[nInputVars] and
221  // abscissas[nInputVars+1] (if present).
222  UInt_t nVars = this->nInputVars();
223  if ( abscissas.size() == nVars+2 ) {
224  Double_t m13Sq = abscissas[nVars];
225  Double_t m23Sq = abscissas[nVars+1];
226  LauKinematics* kinematics = daughters_->getKinematics();
227  if ( dpDependence_ ) {
228  kinematics->updateKinematics( m13Sq, m23Sq );
229  fracVal_ = dpDependence_->calcEfficiency( kinematics );
230  } else {
231  fracVal_ = frac_->value();
232  Double_t dpPos(0.0);
233  if ( dpAxis_ == M12 ) {
234  dpPos = kinematics->calcThirdMassSq(m13Sq,m23Sq);
235  } else if ( dpAxis_ == M13 ) {
236  dpPos = m13Sq;
237  } else if ( dpAxis_ == M23 ) {
238  dpPos = m23Sq;
239  } else if ( dpAxis_ == MMIN ) {
240  dpPos = TMath::Min( m13Sq, m23Sq );
241  } else if ( dpAxis_ == MMAX ) {
242  dpPos = TMath::Max( m13Sq, m23Sq );
243  } else {
244  dpPos = kinematics->distanceFromDPCentre(m13Sq,m23Sq);
245  }
246  this->scaleFrac( dpPos );
247  }
248  }
249 
250  // Add them together
251  Double_t result = fracVal_ * result1 + (1.0-fracVal_) * result2;
252  this->setUnNormPDFVal(result);
253 }
254 
255 void LauDPDepSumPdf::scaleFrac( Double_t dpPos )
256 {
257  Int_t power = 1;
258  for (std::vector<Double_t>::const_iterator iter = fracCoeffs_.begin(); iter != fracCoeffs_.end(); ++iter) {
259  Double_t coeff = (*iter);
260  fracVal_ += coeff * TMath::Power(dpPos,power);
261  ++power;
262  }
263 }
264 
266 {
267  // Nothing to do here, since it is already normalized
268  this->setNorm(1.0);
269 }
270 
272 {
273  // This method gives you the maximum possible height of the PDF.
274  // It combines the maximum heights of the two individual PDFs.
275  // So it would give the true maximum if the two individual maxima coincided.
276  // It is guaranteed to always be >= the true maximum.
277 
278  // Update the heights of the individual PDFs
279  pdf1_->calcPDFHeight( kinematics );
280  pdf2_->calcPDFHeight( kinematics );
281 
282  // Get the up to date parameter values
283  if ( dpDependence_ ) {
284  fracVal_ = dpDependence_->calcEfficiency( kinematics );
285  } else {
286  fracVal_ = frac_->value();
287  Double_t dpPos(0.0);
288  if ( dpAxis_ == M12 ) {
289  dpPos = kinematics->getm12Sq();
290  } else if ( dpAxis_ == M13 ) {
291  dpPos = kinematics->getm13Sq();
292  } else if ( dpAxis_ == M23 ) {
293  dpPos = kinematics->getm23Sq();
294  } else if ( dpAxis_ == MMIN ) {
295  Double_t m13Sq = kinematics->getm13Sq();
296  Double_t m23Sq = kinematics->getm23Sq();
297  dpPos = TMath::Min( m13Sq, m23Sq );
298  } else if ( dpAxis_ == MMAX ) {
299  Double_t m13Sq = kinematics->getm13Sq();
300  Double_t m23Sq = kinematics->getm23Sq();
301  dpPos = TMath::Max( m13Sq, m23Sq );
302  } else {
303  dpPos = kinematics->distanceFromDPCentre();
304  }
305  this->scaleFrac( dpPos );
306  }
307 
308  // Find the (un-normalised) individual PDF maxima
309  Double_t height1 = pdf1_->getMaxHeight();
310  Double_t height2 = pdf2_->getMaxHeight();
311 
312  // Get the individual PDF normalisation factors
313  Double_t norm1 = pdf1_->getNorm();
314  Double_t norm2 = pdf2_->getNorm();
315 
316  // Calculate the normalised individual PDF maxima
317  height1 /= norm1;
318  height2 /= norm2;
319 
320  // Combine these heights together
321  Double_t height = fracVal_ * height1 + (1-fracVal_) * height2;
322  this->setMaxHeight(height);
323 }
324 
325 // Override the base class methods for cacheInfo and calcLikelihoodInfo(UInt_t iEvt).
326 // For both of these we delegate to the two constituent PDFs.
327 
329 {
330  // delegate to the two sub-PDFs to cache their information
331  pdf1_->cacheInfo(inputData);
332  pdf2_->cacheInfo(inputData);
333 
334  // now check that the DP variables are included in the data
335  Bool_t hasBranch = inputData.haveBranch( "m13Sq" );
336  hasBranch &= inputData.haveBranch( "m23Sq" );
337  if (!hasBranch) {
338  cerr<<"ERROR in LauDPDepSumPdf::cacheInfo : Input data does not contain Dalitz plot variables m13Sq and/or m23Sq."<<endl;
339  return;
340  }
341 
342  // determine whether we are caching our PDF value
343  Bool_t doCaching( this->nFixedParameters() == this->nParameters() );
344  this->cachePDF( doCaching );
345 
346  if ( !doCaching ) {
347  // in this case we seem to be doing a fit where the parameters are floating
348  // so need to mark that the PDF height is no longer up to date
349  this->heightUpToDate(kFALSE);
350  }
351 
352  // clear the vector and reserve enough space
353  UInt_t nEvents = inputData.nEvents();
354  fractions_.clear();
355  fractions_.reserve(nEvents);
356 
357  // loop through the events, determine the fraction and store
358  for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
359 
360  const LauFitData& dataValues = inputData.getData(iEvt);
361 
362  Double_t m13Sq = dataValues.find("m13Sq")->second;
363  Double_t m23Sq = dataValues.find("m23Sq")->second;
364 
365  LauKinematics* kinematics = daughters_->getKinematics();
366  if ( dpDependence_ ) {
367  // if we're using the histogram then just
368  // determine the fraction and store
369  kinematics->updateKinematics( m13Sq, m23Sq );
370  fracVal_ = dpDependence_->calcEfficiency( kinematics );
371  } else {
372  // if we're scaling the fraction parameter then we
373  // just store the scaling info since the parameter
374  // might be floating
375  fracVal_ = frac_->value();
376  Double_t dpPos(0.0);
377  if ( dpAxis_ == M12 ) {
378  dpPos = kinematics->calcThirdMassSq(m13Sq,m23Sq);
379  } else if ( dpAxis_ == M13 ) {
380  dpPos = m13Sq;
381  } else if ( dpAxis_ == M23 ) {
382  dpPos = m23Sq;
383  } else if ( dpAxis_ == MMIN ) {
384  dpPos = TMath::Min( m13Sq, m23Sq );
385  } else if ( dpAxis_ == MMAX ) {
386  dpPos = TMath::Max( m13Sq, m23Sq );
387  } else {
388  dpPos = kinematics->distanceFromDPCentre(m13Sq,m23Sq);
389  }
390  this->scaleFrac( dpPos );
391  fracVal_ -= frac_->value();
392  }
393 
394  fractions_.push_back( fracVal_ );
395  }
396 }
397 
399 {
400  // Get the fraction value for this event
401  fracVal_ = fractions_[iEvt];
402  if ( frac_ ) {
403  // if we're scaling the parameter then need to add the
404  // current value of the parameter
405  fracVal_ += frac_->value();
406  }
407 
408  // Evaluate the normalised PDF values
409  pdf1_->calcLikelihoodInfo(iEvt);
410  pdf2_->calcLikelihoodInfo(iEvt);
411  Double_t result1 = pdf1_->getLikelihood();
412  Double_t result2 = pdf2_->getLikelihood();
413 
414  // Add them together
415  Double_t result = fracVal_ * result1 + (1-fracVal_) * result2;
416  this->setUnNormPDFVal(result);
417 }
418 
virtual void cacheInfo(const LauFitDataTree &inputData)
Cache information from data.
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
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:369
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:117
DPAxis
Define possibilties for the DP axes.
LauDPDepSumPdf(LauAbsPdf *pdf1, LauAbsPdf *pdf2, const LauDaughters *daughters, const TH2 *dpHisto, Bool_t upperHalf=kFALSE, Bool_t useSpline=kFALSE)
Constructor - fraction determined by 2D histogram.
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
virtual LauAbsRValue * findParameter(const TString &parName)
Retrieve the specified parameter.
Definition: LauAbsPdf.cc:381
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:92
virtual void calcNorm()
Calculate the normalisation.
LauAbsPdf * pdf2_
Second PDF.
File containing declaration of LauDaughters class.
void scaleFrac(Double_t dpPos)
Scale fraction according to DP position.
virtual void calcPDFHeight(const LauKinematics *kinematics)=0
Calculate the maximum height of the PDF.
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:325
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:213
LauAbsRValue * frac_
Fractional contribution of first PDF to the final PDF.
Class for defining a PDF that is the DP-dependent sum of two other PDFs.
std::map< TString, Double_t > LauFitData
Type for holding event data.
Double_t getm23Sq() const
Get the m23 invariant mass square.
Double_t distanceFromDPCentre() const
Calculate the distance from the currently set (m13Sq, m23Sq) point to the centre of the Dalitz plot (...
std::vector< Double_t > fractions_
Cached values of fractions.
virtual Double_t getMaxHeight() const
Retrieve the maximum height.
Definition: LauAbsPdf.hh:221
virtual ~LauDPDepSumPdf()
Destructor.
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
File containing declaration of LauDPDepSumPdf class.
LauAbsPdf * pdf1_
First PDF.
Double_t calcEfficiency(const LauKinematics *kinematics) const
Determine the efficiency for a given point in the Dalitz plot.
Definition: LauEffModel.cc:261
File containing declaration of LauParameter class.
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
LauKinematics * getKinematics()
Retrieve the Dalitz plot kinematics.
virtual std::vector< TString > varNames() const
Retrieve the names of the abscissas.
Definition: LauAbsPdf.cc:103
const std::vector< Double_t > fracCoeffs_
Polynomial used to scale fraction.
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:123
Double_t fracVal_
Fractional contribution of first PDF to the final PDF.
virtual Bool_t isDPDependent() const
Specifies whether or not the PDF is DP dependent.
Definition: LauAbsPdf.hh:111
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
virtual void addParameters(std::vector< LauAbsRValue * > &params)
Add parameters to the PDF.
Definition: LauAbsPdf.cc:423
DPAxis dpAxis_
The DP axis we depend on.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
virtual const std::vector< LauAbsRValue * > & getParameters() const
Retrieve the parameters of the PDF, e.g. so that they can be loaded into a fit.
Definition: LauAbsPdf.hh:239
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
File containing declaration of LauEffModel class.
Double_t getm12Sq() const
Get the m12 invariant mass square.
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:276
Class that implements the efficiency description across the signal Dalitz plot.
Definition: LauEffModel.hh:37
virtual Double_t getLikelihood() const
Retrieve the normalised likelihood value.
Definition: LauAbsPdf.cc:354
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) ...
File containing LauConstants namespace.
virtual Double_t value() const =0
Return the value of the parameter.
Bool_t haveBranch(const TString &name) const
Check if the named branch is stored.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
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 cacheInfo(const LauFitDataTree &inputData)
Cache information from data.
Definition: LauAbsPdf.cc:241
Double_t value() const
The value of the parameter.
LauDaughters * daughters_
Daughter particles.
UInt_t nEvents() const
Retrieve the number of events.
LauEffModel * dpDependence_
DP dependence.
virtual UInt_t nInputVars() const
Retrieve the number of abscissas.
Definition: LauAbsPdf.hh:104
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:29
virtual Double_t getNorm() const
Retrieve the normalisation factor.
Definition: LauAbsPdf.hh:202
Class to store the input fit variables.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45