laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2p1
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 
103 LauSumPdf::LauSumPdf(const LauSumPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
104 {
105  // Copy constructor
106  pdf1_ = other.pdf1_;
107  pdf2_ = other.pdf2_;
108  this->setRandomFun(other.getRandomFun());
109  this->calcNorm();
110 }
111 
113 {
114  // Check that the given abscissa is within the allowed range
115  if (!this->checkRange(abscissas)) {
116  gSystem->Exit(EXIT_FAILURE);
117  }
118 
119  // Get the up to date parameter values
120  Double_t frac = frac_->value();
121 
122  // Evaluate the normalised PDF values
123  pdf1_->calcLikelihoodInfo(abscissas);
124  pdf2_->calcLikelihoodInfo(abscissas);
125  Double_t result1 = pdf1_->getLikelihood();
126  Double_t result2 = pdf2_->getLikelihood();
127 
128  // Add them together
129  Double_t result = frac * result1 + (1.0-frac) * result2;
130  this->setUnNormPDFVal(result);
131 }
132 
134 {
135  // Nothing to do here, since it is already normalized
136  this->setNorm(1.0);
137 }
138 
139 void LauSumPdf::calcPDFHeight( const LauKinematics* kinematics )
140 {
141  if ( this->heightUpToDate() && ! this->isDPDependent() ) {
142  return;
143  }
144 
145  // This method gives you the maximum possible height of the PDF.
146  // It combines the maximum heights of the two individual PDFs.
147  // So it would give the true maximum if the two individual maxima coincided.
148  // It is guaranteed to always be >= the true maximum.
149 
150  // Update the heights of the individual PDFs
151  pdf1_->calcPDFHeight( kinematics );
152  pdf2_->calcPDFHeight( kinematics );
153 
154  // Get the up to date parameter values
155  Double_t frac = frac_->value();
156 
157  // Find the (un-normalised) individual PDF maxima
158  Double_t height1 = pdf1_->getMaxHeight();
159  Double_t height2 = pdf2_->getMaxHeight();
160 
161  // Get the individual PDF normalisation factors
162  Double_t norm1 = pdf1_->getNorm();
163  Double_t norm2 = pdf2_->getNorm();
164 
165  // Calculate the normalised individual PDF maxima
166  height1 /= norm1;
167  height2 /= norm2;
168 
169  // Combine these heights together
170  Double_t height = frac * height1 + (1-frac) * height2;
171  this->setMaxHeight(height);
172 }
173 
174 // Override the base class methods for cacheInfo and calcLikelihoodInfo(UInt_t iEvt).
175 // For both of these we delegate to the two constituent PDFs.
176 
177 void LauSumPdf::cacheInfo(const LauFitDataTree& inputData)
178 {
179  pdf1_->cacheInfo(inputData);
180  pdf2_->cacheInfo(inputData);
181 
182  // determine whether we are caching our PDF value
183  Bool_t doCaching ( this->nFixedParameters() == this->nParameters() );
184  this->cachePDF( doCaching );
185 
186  // Set whether the height is up-to-date only if it is true for both of the sub-PDFs
187  Bool_t hutd = pdf1_->heightUpToDate() && pdf2_->heightUpToDate();
188  this->heightUpToDate(hutd);
189 }
190 
192 {
193  // Get the up to date parameter values
194  Double_t frac = frac_->value();
195 
196  // Evaluate the normalised PDF values
197  pdf1_->calcLikelihoodInfo(iEvt);
198  pdf2_->calcLikelihoodInfo(iEvt);
199  Double_t result1 = pdf1_->getLikelihood();
200  Double_t result2 = pdf2_->getLikelihood();
201 
202  // Add them together
203  Double_t result = frac * result1 + (1.0-frac) * result2;
204  this->setUnNormPDFVal(result);
205 }
206 
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:92
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:387
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:133
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
Class for defining the fit parameter objects.
Definition: LauParameter.hh:33
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:112
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:95
File containing LauConstants namespace.
virtual Double_t value() const =0
Return the value of the parameter.
LauSumPdf(LauAbsPdf *pdf1, LauAbsPdf *pdf2, LauParameter *frac1)
Constructor.
Definition: LauSumPdf.cc:28
LauAbsRValue * frac_
Fractional contribution of first PDF to the final PDF.
Definition: LauSumPdf.hh:98
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:139
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:59
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:233
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:177