laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepGaussPdf.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 "TMath.h"
22 #include "TSystem.h"
23 
24 #include "LauConstants.hh"
25 #include "LauDaughters.hh"
26 #include "LauDPDepGaussPdf.hh"
27 #include "LauKinematics.hh"
28 
30 
31 
32 LauDPDepGaussPdf::LauDPDepGaussPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params,
33  Double_t minAbscissa, Double_t maxAbscissa,
34  const LauDaughters* daughters,
35  const std::vector<Double_t>& meanCoeffs,
36  const std::vector<Double_t>& sigmaCoeffs,
37  DPAxis dpAxis) :
38  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
39  kinematics_(daughters ? daughters->getKinematics() : 0),
40  mean_(0),
41  sigma_(0),
42  meanVal_(0.0),
43  sigmaVal_(0.0),
44  meanCoeffs_(meanCoeffs),
45  sigmaCoeffs_(sigmaCoeffs),
46  dpAxis_(dpAxis)
47 {
48  // Constructor for the Gaussian PDF.
49 
50  // Check we have a valid kinematics object
51  if ( ! kinematics_ ) {
52  cerr<<"ERROR in LauDPDepGaussPdf::LauDPDepGaussPdf : Have not been provided with a valid DP kinematics object."<<endl;
53  gSystem->Exit(EXIT_FAILURE);
54  }
55 
56  // The parameters in params are the mean and the sigma (half the width) of the gaussian.
57  // The last two arguments specify the range in which the PDF is defined, and the PDF
58  // will be normalised w.r.t. these limits.
59 
60  mean_ = this->findParameter("mean");
61  sigma_ = this->findParameter("sigma");
62 
63  if ((this->nParameters() != 2) || (mean_ == 0) || (sigma_ == 0)) {
64  cerr<<"ERROR in LauDPDepGaussPdf constructor: LauDPDepGaussPdf requires 2 parameters: \"mean\" and \"sigma\"."<<endl;
65  gSystem->Exit(EXIT_FAILURE);
66  }
67 
68  // Cache the normalisation factor
69  this->calcNorm();
70 }
71 
73 {
74  // Destructor
75 }
76 
78 {
79  this->setNorm(1.0);
80 }
81 
83 {
84  // Check that the given abscissa is within the allowed range
85  if (!this->checkRange(abscissas)) {
86  gSystem->Exit(EXIT_FAILURE);
87  }
88 
89  // Get our abscissa
90  Double_t abscissa = abscissas[0];
91 
92  // Get the up to date parameter values
93  meanVal_ = mean_->value();
94  sigmaVal_ = sigma_->value();
95 
96  // Find out the DP position
97  Double_t dpPos(0.0);
98  UInt_t nVars = this->nInputVars();
99  if ( abscissas.size() == nVars+2 ) {
100  Double_t m13Sq = abscissas[nVars];
101  Double_t m23Sq = abscissas[nVars+1];
102 
103  if ( dpAxis_ == M12 ) {
104  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
105  } else if ( dpAxis_ == M13 ) {
106  dpPos = m13Sq;
107  } else if ( dpAxis_ == M23 ) {
108  dpPos = m23Sq;
109  } else if ( dpAxis_ == MMIN ) {
110  dpPos = TMath::Min( m13Sq, m23Sq );
111  } else if ( dpAxis_ == MMAX ) {
112  dpPos = TMath::Max( m13Sq, m23Sq );
113  } else {
114  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
115  }
116  }
117 
118  // Scale the parameters according to the DP position
119  this->scalePars( dpPos );
120 
121  // Calculate the value of the Gaussian for the given value of the abscissa.
122  Double_t arg = abscissa - meanVal_;
123 
124  Double_t exponent(0.0);
125  if (TMath::Abs(sigmaVal_) > 1e-10) {
126  exponent = -0.5*arg*arg/(sigmaVal_*sigmaVal_);
127  }
128 
129  Double_t value = TMath::Exp(exponent);
130 
131  Double_t norm(0.0);
132  Double_t scale = LauConstants::root2*sigmaVal_;
133  if (TMath::Abs(sigmaVal_) > 1e-10) {
134  norm = LauConstants::rootPiBy2*sigmaVal_*(TMath::Erf((this->getMaxAbscissa() - meanVal_)/scale) - TMath::Erf((this->getMinAbscissa() - meanVal_)/scale));
135  }
136 
137  value /= norm;
138 
139  this->setUnNormPDFVal(value);
140 }
141 
142 void LauDPDepGaussPdf::scalePars(Double_t dpPos)
143 {
144  Int_t power = 1;
145  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
146  Double_t coeff = (*iter);
147  meanVal_ += coeff * TMath::Power(dpPos,power);
148  ++power;
149  }
150 
151  power = 1;
152  for (std::vector<Double_t>::const_iterator iter = sigmaCoeffs_.begin(); iter != sigmaCoeffs_.end(); ++iter) {
153  Double_t coeff = (*iter);
154  sigmaVal_ += coeff * TMath::Power(dpPos,power);
155  ++power;
156  }
157 }
158 
160 {
161  // Get the up to date parameter values
162  meanVal_ = mean_->value();
163 
164  // Find out the DP position
165  Double_t dpPos(0.0);
166  if ( dpAxis_ == M12 ) {
167  dpPos = kinematics->getm12Sq();
168  } else if ( dpAxis_ == M13 ) {
169  dpPos = kinematics->getm13Sq();
170  } else if ( dpAxis_ == M23 ) {
171  dpPos = kinematics->getm23Sq();
172  } else if ( dpAxis_ == MMIN ) {
173  Double_t m13Sq = kinematics->getm13Sq();
174  Double_t m23Sq = kinematics->getm23Sq();
175  dpPos = TMath::Min( m13Sq, m23Sq );
176  } else if ( dpAxis_ == MMAX ) {
177  Double_t m13Sq = kinematics->getm13Sq();
178  Double_t m23Sq = kinematics->getm23Sq();
179  dpPos = TMath::Max( m13Sq, m23Sq );
180  } else {
181  dpPos = kinematics->distanceFromDPCentre();
182  }
183 
184  // Scale the parameters according to the DP position
185  this->scalePars( dpPos );
186 
187  LauAbscissas maxPoint(3);
188  maxPoint[0] = meanVal_;
189  maxPoint[1] = kinematics->getm13Sq();
190  maxPoint[2] = kinematics->getm23Sq();
191 
192  // Calculate the PDF height for the Gaussian function.
193  if (meanVal_>this->getMaxAbscissa()) {
194  maxPoint[0] = this->getMaxAbscissa();
195  } else if (meanVal_<this->getMinAbscissa()) {
196  maxPoint[0] = this->getMinAbscissa();
197  }
198  this->calcLikelihoodInfo(maxPoint);
199 
200  Double_t height = this->getUnNormLikelihood();
201  this->setMaxHeight(height);
202 }
203 
Double_t sigmaVal_
Gaussian sigma.
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
virtual ~LauDPDepGaussPdf()
Destructor.
ClassImp(LauAbsCoeffSet)
DPAxis dpAxis_
The DP axis we depend on.
DPAxis
Define possibilties for the DP axes.
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
File containing declaration of LauDaughters class.
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
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 (...
File containing declaration of LauKinematics class.
File containing declaration of LauDPDepGaussPdf class.
const Double_t rootPiBy2
Square root of Pi divided by two.
virtual void calcNorm()
Calculate the normalisation.
const std::vector< Double_t > sigmaCoeffs_
Coefficients of Gaussian sigma.
const std::vector< Double_t > meanCoeffs_
Coefficients of Gaussian mean.
const Double_t root2
Square root of two.
LauAbsRValue * mean_
Gaussian mean.
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:123
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
Double_t meanVal_
Gaussian mean.
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
Double_t getm12Sq() const
Get the m12 invariant mass square.
Double_t getm13Sq() const
Get the m13 invariant mass square.
File containing LauConstants namespace.
virtual Double_t value() const =0
Return the value of the parameter.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
const LauKinematics * kinematics_
The current DP kinematics.
void scalePars(Double_t dpPos)
Scale parameters by their dependence on the DP position.
virtual UInt_t nInputVars() const
Retrieve the number of abscissas.
Definition: LauAbsPdf.hh:104
Class for defining a Gaussian PDF (DP dependent).
LauAbsRValue * sigma_
Gaussian sigma.
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:29
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45