laura is hosted by Hepforge, IPPP Durham
Laura++  v1r2
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 
29 ClassImp(LauDPDepGaussPdf)
30 
31 
32 LauDPDepGaussPdf::LauDPDepGaussPdf(const TString& theVarName, const std::vector<LauParameter*>& 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 
77 LauDPDepGaussPdf::LauDPDepGaussPdf(const LauDPDepGaussPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa()),
78  kinematics_( other.kinematics_ ),
79  mean_(0),
80  sigma_(0),
81  meanVal_(0.0),
82  sigmaVal_(0.0),
83  meanCoeffs_( other.meanCoeffs_ ),
84  sigmaCoeffs_( other.sigmaCoeffs_ ),
85  dpAxis_( other.dpAxis_ )
86 {
87  // Copy constructor
88  mean_ = this->findParameter("mean");
89  sigma_ = this->findParameter("sigma");
90  this->setRandomFun(other.getRandomFun());
91  this->calcNorm();
92 }
93 
95 {
96  this->setNorm(1.0);
97 }
98 
100 {
101  // Check that the given abscissa is within the allowed range
102  if (!this->checkRange(abscissas)) {
103  gSystem->Exit(EXIT_FAILURE);
104  }
105 
106  // Get our abscissa
107  Double_t abscissa = abscissas[0];
108 
109  // Get the up to date parameter values
110  meanVal_ = mean_->value();
111  sigmaVal_ = sigma_->value();
112 
113  // Find out the DP position
114  Double_t dpPos(0.0);
115  UInt_t nVars = this->nInputVars();
116  if ( abscissas.size() == nVars+2 ) {
117  Double_t m13Sq = abscissas[nVars];
118  Double_t m23Sq = abscissas[nVars+1];
119 
120  if ( dpAxis_ == M12 ) {
121  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
122  } else if ( dpAxis_ == M13 ) {
123  dpPos = m13Sq;
124  } else if ( dpAxis_ == M23 ) {
125  dpPos = m23Sq;
126  } else if ( dpAxis_ == MMIN ) {
127  dpPos = TMath::Min( m13Sq, m23Sq );
128  } else if ( dpAxis_ == MMAX ) {
129  dpPos = TMath::Max( m13Sq, m23Sq );
130  } else {
131  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
132  }
133  }
134 
135  // Scale the parameters according to the DP position
136  this->scalePars( dpPos );
137 
138  // Calculate the value of the Gaussian for the given value of the abscissa.
139  Double_t arg = abscissa - meanVal_;
140 
141  Double_t exponent(0.0);
142  if (TMath::Abs(sigmaVal_) > 1e-10) {
143  exponent = -0.5*arg*arg/(sigmaVal_*sigmaVal_);
144  }
145 
146  Double_t value = TMath::Exp(exponent);
147 
148  Double_t norm(0.0);
149  Double_t scale = LauConstants::root2*sigmaVal_;
150  if (TMath::Abs(sigmaVal_) > 1e-10) {
151  norm = LauConstants::rootPiBy2*sigmaVal_*(TMath::Erf((this->getMaxAbscissa() - meanVal_)/scale) - TMath::Erf((this->getMinAbscissa() - meanVal_)/scale));
152  }
153 
154  value /= norm;
155 
156  this->setUnNormPDFVal(value);
157 }
158 
159 void LauDPDepGaussPdf::scalePars(Double_t dpPos)
160 {
161  Int_t power = 1;
162  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
163  Double_t coeff = (*iter);
164  meanVal_ += coeff * TMath::Power(dpPos,power);
165  ++power;
166  }
167 
168  power = 1;
169  for (std::vector<Double_t>::const_iterator iter = sigmaCoeffs_.begin(); iter != sigmaCoeffs_.end(); ++iter) {
170  Double_t coeff = (*iter);
171  sigmaVal_ += coeff * TMath::Power(dpPos,power);
172  ++power;
173  }
174 }
175 
177 {
178  // Get the up to date parameter values
179  meanVal_ = mean_->value();
180 
181  // Find out the DP position
182  Double_t dpPos(0.0);
183  if ( dpAxis_ == M12 ) {
184  dpPos = kinematics->getm12Sq();
185  } else if ( dpAxis_ == M13 ) {
186  dpPos = kinematics->getm13Sq();
187  } else if ( dpAxis_ == M23 ) {
188  dpPos = kinematics->getm23Sq();
189  } else if ( dpAxis_ == MMIN ) {
190  Double_t m13Sq = kinematics->getm13Sq();
191  Double_t m23Sq = kinematics->getm23Sq();
192  dpPos = TMath::Min( m13Sq, m23Sq );
193  } else if ( dpAxis_ == MMAX ) {
194  Double_t m13Sq = kinematics->getm13Sq();
195  Double_t m23Sq = kinematics->getm23Sq();
196  dpPos = TMath::Max( m13Sq, m23Sq );
197  } else {
198  dpPos = kinematics->distanceFromDPCentre();
199  }
200 
201  // Scale the parameters according to the DP position
202  this->scalePars( dpPos );
203 
204  LauAbscissas maxPoint(3);
205  maxPoint[0] = meanVal_;
206  maxPoint[1] = kinematics->getm13Sq();
207  maxPoint[2] = kinematics->getm23Sq();
208 
209  // Calculate the PDF height for the Gaussian function.
210  if (meanVal_>this->getMaxAbscissa()) {
211  maxPoint[0] = this->getMaxAbscissa();
212  } else if (meanVal_<this->getMinAbscissa()) {
213  maxPoint[0] = this->getMinAbscissa();
214  }
215  this->calcLikelihoodInfo(maxPoint);
216 
217  Double_t height = this->getUnNormLikelihood();
218  this->setMaxHeight(height);
219 }
220 
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:454
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:158
virtual ~LauDPDepGaussPdf()
Destructor.
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:278
File containing declaration of LauDaughters class.
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
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.
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:472
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.
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:164
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:416
Double_t meanVal_
Gaussian mean.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:31
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.
LauParameter * mean_
Gaussian mean.
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 LauParameter * findParameter(const TString &parName)
Retrieve the specified parameter.
Definition: LauAbsPdf.cc:381
LauDPDepGaussPdf(const TString &theVarName, const std::vector< LauParameter * > &params, Double_t minAbscissa, Double_t maxAbscissa, const LauDaughters *daughters, const std::vector< Double_t > &meanCoeffs, const std::vector< Double_t > &sigmaCoeffs, DPAxis dpAxis)
Constructor.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:40
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:103
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:315
Class for defining a Gaussian PDF (DP dependent).
LauParameter * sigma_
Gaussian sigma.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:44