laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauSumPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 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 
18 #include "TMath.h"
19 #include "TSystem.h"
20 
21 #include "LauConstants.hh"
22 #include "LauParameter.hh"
23 #include "LauSumPdf.hh"
24 
26 
27 
29  LauAbsPdf(pdf1 ? pdf1->varNames() : std::vector<TString>(), std::vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
30  pdf1_(pdf1),
31  pdf2_(pdf2),
32  frac_(frac)
33 {
34  // Constructor for the sum PDF.
35  // We are defining the sum as:
36  // f x (PDF1/S(PDF1)) + (1-f) x (PDF2/S(PDF2))
37  // where f is the fraction, x is multiplication, PDFi is the i'th PDF,
38  // and S(PDFi) is the integral of the i'th PDF.
39  //
40  // The last two arguments specify the range in which the PDF is defined, and the PDF
41  // will be normalised w.r.t. these limits.
42 
43  // So the first thing we have to do is check the pointers are all valid.
44  if (!frac || !pdf1 || !pdf2) {
45  std::cerr << "ERROR in LauSumPdf constructor : one of the 3 pointers is null." << std::endl;
46  gSystem->Exit(EXIT_FAILURE);
47  }
48 
49  // Next check that the abscissa ranges are the same for each PDF
50  if (pdf1->getMinAbscissa() != pdf2->getMinAbscissa()) {
51  std::cerr << "ERROR in LauSumPdf constructor : minimum abscissa values not the same for the two PDFs." << std::endl;
52  gSystem->Exit(EXIT_FAILURE);
53  }
54  if (pdf1->getMaxAbscissa() != pdf2->getMaxAbscissa()) {
55  std::cerr << "ERROR in LauSumPdf constructor : maximum abscissa values not the same for the two PDFs." << std::endl;
56  gSystem->Exit(EXIT_FAILURE);
57  }
58 
59  // Also check that both PDFs are expecting the same number of input variables
60  if (pdf1->nInputVars() != pdf2->nInputVars()) {
61  std::cerr << "ERROR in LauSumPdf constructor : number of input variables not the same for the two PDFs." << std::endl;
62  gSystem->Exit(EXIT_FAILURE);
63  }
64 
65  // Also check that both PDFs are expecting the same variable name
66  if (pdf1->varNames() != pdf2->varNames()) {
67  std::cerr << "ERROR in LauSumPdf constructor : variable name not the same for the two PDFs." << std::endl;
68  gSystem->Exit(EXIT_FAILURE);
69  }
70 
71  // Then we need to grab all the parameters and pass them to the base class.
72  // This is so that when we are asked for them they can be put into the fit.
73  // The number of parameters is the number in PDF1 + the number in PDF2 + 1 for the fraction.
74  UInt_t nPar(pdf1->nParameters()+pdf2->nParameters()+1);
75  std::vector<LauAbsRValue*> params; params.reserve(nPar);
76  std::vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
77  std::vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
78  params.push_back(frac);
79  for (std::vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
80  params.push_back(*iter);
81  }
82  for (std::vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
83  params.push_back(*iter);
84  }
85  this->addParameters(params);
86 
87  // Now check that we can find the fraction parameter ok
88  frac_ = this->findParameter("frac");
89  if (frac_ == 0) {
90  std::cerr << "ERROR in LauSumPdf constructor : parameter \"frac\" not found." << std::endl;
91  gSystem->Exit(EXIT_FAILURE);
92  }
93 
94  // Cache the normalisation factor
95  this->calcNorm();
96 }
97 
99 {
100  // Destructor
101 }
102 
104 {
105  // Check that the given abscissa is within the allowed range
106  if (!this->checkRange(abscissas)) {
107  gSystem->Exit(EXIT_FAILURE);
108  }
109 
110  // Get the up to date parameter values
111  Double_t frac = frac_->value();
112 
113  // Evaluate the normalised PDF values
114  pdf1_->calcLikelihoodInfo(abscissas);
115  pdf2_->calcLikelihoodInfo(abscissas);
116  Double_t result1 = pdf1_->getLikelihood();
117  Double_t result2 = pdf2_->getLikelihood();
118 
119  // Add them together
120  Double_t result = frac * result1 + (1.0-frac) * result2;
121  this->setUnNormPDFVal(result);
122 }
123 
125 {
126  // Nothing to do here, since it is already normalized
127  this->setNorm(1.0);
128 }
129 
130 void LauSumPdf::calcPDFHeight( const LauKinematics* kinematics )
131 {
132  if ( this->heightUpToDate() && ! this->isDPDependent() ) {
133  return;
134  }
135 
136  // This method gives you the maximum possible height of the PDF.
137  // It combines the maximum heights of the two individual PDFs.
138  // So it would give the true maximum if the two individual maxima coincided.
139  // It is guaranteed to always be >= the true maximum.
140 
141  // Update the heights of the individual PDFs
142  pdf1_->calcPDFHeight( kinematics );
143  pdf2_->calcPDFHeight( kinematics );
144 
145  // Get the up to date parameter values
146  Double_t frac = frac_->value();
147 
148  // Find the (un-normalised) individual PDF maxima
149  Double_t height1 = pdf1_->getMaxHeight();
150  Double_t height2 = pdf2_->getMaxHeight();
151 
152  // Get the individual PDF normalisation factors
153  Double_t norm1 = pdf1_->getNorm();
154  Double_t norm2 = pdf2_->getNorm();
155 
156  // Calculate the normalised individual PDF maxima
157  height1 /= norm1;
158  height2 /= norm2;
159 
160  // Combine these heights together
161  Double_t height = frac * height1 + (1-frac) * height2;
162  this->setMaxHeight(height);
163 }
164 
165 // Override the base class methods for cacheInfo and calcLikelihoodInfo(UInt_t iEvt).
166 // For both of these we delegate to the two constituent PDFs.
167 
168 void LauSumPdf::cacheInfo(const LauFitDataTree& inputData)
169 {
170  pdf1_->cacheInfo(inputData);
171  pdf2_->cacheInfo(inputData);
172 
173  // determine whether we are caching our PDF value
174  Bool_t doCaching ( this->nFixedParameters() == this->nParameters() );
175  this->cachePDF( doCaching );
176 
177  // Set whether the height is up-to-date only if it is true for both of the sub-PDFs
178  Bool_t hutd = pdf1_->heightUpToDate() && pdf2_->heightUpToDate();
179  this->heightUpToDate(hutd);
180 }
181 
183 {
184  // Get the up to date parameter values
185  Double_t frac = frac_->value();
186 
187  // Evaluate the normalised PDF values
188  pdf1_->calcLikelihoodInfo(iEvt);
189  pdf2_->calcLikelihoodInfo(iEvt);
190  Double_t result1 = pdf1_->getLikelihood();
191  Double_t result2 = pdf2_->getLikelihood();
192 
193  // Add them together
194  Double_t result = frac * result1 + (1.0-frac) * result2;
195  this->setUnNormPDFVal(result);
196 }
197 
File containing declaration of LauSumPdf class.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:369
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
virtual UInt_t nParameters() const
Retrieve the number of PDF parameters.
Definition: LauAbsPdf.hh:92
virtual ~LauSumPdf()
Destructor.
Definition: LauSumPdf.cc:98
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
std::map< TString, Double_t > LauFitData
Type for holding event data.
LauAbsPdf * pdf1_
First PDF.
Definition: LauSumPdf.hh:95
virtual Double_t getMaxHeight() const
Retrieve the maximum height.
Definition: LauAbsPdf.hh:221
File containing declaration of LauParameter class.
virtual void calcNorm()
Calculate the normalisation.
Definition: LauSumPdf.cc:124
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
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
Class for defining a PDF that is the sum of two other PDFs.
Definition: LauSumPdf.hh:32
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:276
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
Definition: LauSumPdf.cc:103
virtual Double_t getLikelihood() const
Retrieve the normalised likelihood value.
Definition: LauAbsPdf.cc:354
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)=0
Calculate the likelihood (and all associated information) given value(s) of the abscissa(s) ...
LauAbsPdf * pdf2_
Second PDF.
Definition: LauSumPdf.hh:98
File containing LauConstants namespace.
virtual Double_t value() const =0
Return the value of the parameter.
LauAbsRValue * frac_
Fractional contribution of first PDF to the final PDF.
Definition: LauSumPdf.hh:101
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Definition: LauSumPdf.cc:130
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
virtual Bool_t isDPDependent() const
Boolean for the DP dependence of PDFs in the sum.
Definition: LauSumPdf.hh:56
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
virtual void cacheInfo(const LauFitDataTree &inputData)
Cache information from data.
Definition: LauSumPdf.cc:168