laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepGaussPdf.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2009 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include <iostream>
30 #include <vector>
31 using std::cout;
32 using std::cerr;
33 using std::endl;
34 
35 #include "TMath.h"
36 #include "TSystem.h"
37 
38 #include "LauConstants.hh"
39 #include "LauDaughters.hh"
40 #include "LauDPDepGaussPdf.hh"
41 #include "LauKinematics.hh"
42 
44 
45 
46 LauDPDepGaussPdf::LauDPDepGaussPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params,
47  Double_t minAbscissa, Double_t maxAbscissa,
48  const LauDaughters* daughters,
49  const std::vector<Double_t>& meanCoeffs,
50  const std::vector<Double_t>& sigmaCoeffs,
51  DPAxis dpAxis) :
52  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
53  kinematics_(daughters ? daughters->getKinematics() : 0),
54  mean_(0),
55  sigma_(0),
56  meanVal_(0.0),
57  sigmaVal_(0.0),
58  meanCoeffs_(meanCoeffs),
59  sigmaCoeffs_(sigmaCoeffs),
60  dpAxis_(dpAxis)
61 {
62  // Constructor for the Gaussian PDF.
63 
64  // Check we have a valid kinematics object
65  if ( ! kinematics_ ) {
66  cerr<<"ERROR in LauDPDepGaussPdf::LauDPDepGaussPdf : Have not been provided with a valid DP kinematics object."<<endl;
67  gSystem->Exit(EXIT_FAILURE);
68  }
69 
70  // The parameters in params are the mean and the sigma (half the width) of the gaussian.
71  // The last two arguments specify the range in which the PDF is defined, and the PDF
72  // will be normalised w.r.t. these limits.
73 
74  mean_ = this->findParameter("mean");
75  sigma_ = this->findParameter("sigma");
76 
77  if ((this->nParameters() != 2) || (mean_ == 0) || (sigma_ == 0)) {
78  cerr<<"ERROR in LauDPDepGaussPdf constructor: LauDPDepGaussPdf requires 2 parameters: \"mean\" and \"sigma\"."<<endl;
79  gSystem->Exit(EXIT_FAILURE);
80  }
81 
82  // Cache the normalisation factor
83  this->calcNorm();
84 }
85 
87 {
88  // Destructor
89 }
90 
92 {
93  this->setNorm(1.0);
94 }
95 
97 {
98  // Check that the given abscissa is within the allowed range
99  if (!this->checkRange(abscissas)) {
100  gSystem->Exit(EXIT_FAILURE);
101  }
102 
103  // Get our abscissa
104  Double_t abscissa = abscissas[0];
105 
106  // Get the up to date parameter values
109 
110  // Find out the DP position
111  Double_t dpPos(0.0);
112  UInt_t nVars = this->nInputVars();
113  if ( abscissas.size() == nVars+2 ) {
114  Double_t m13Sq = abscissas[nVars];
115  Double_t m23Sq = abscissas[nVars+1];
116 
117  if ( dpAxis_ == M12 ) {
118  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
119  } else if ( dpAxis_ == M13 ) {
120  dpPos = m13Sq;
121  } else if ( dpAxis_ == M23 ) {
122  dpPos = m23Sq;
123  } else if ( dpAxis_ == MMIN ) {
124  dpPos = TMath::Min( m13Sq, m23Sq );
125  } else if ( dpAxis_ == MMAX ) {
126  dpPos = TMath::Max( m13Sq, m23Sq );
127  } else {
128  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
129  }
130  }
131 
132  // Scale the parameters according to the DP position
133  this->scalePars( dpPos );
134 
135  // Calculate the value of the Gaussian for the given value of the abscissa.
136  Double_t arg = abscissa - meanVal_;
137 
138  Double_t exponent(0.0);
139  if (TMath::Abs(sigmaVal_) > 1e-10) {
140  exponent = -0.5*arg*arg/(sigmaVal_*sigmaVal_);
141  }
142 
143  Double_t value = TMath::Exp(exponent);
144 
145  Double_t norm(0.0);
146  Double_t scale = LauConstants::root2*sigmaVal_;
147  if (TMath::Abs(sigmaVal_) > 1e-10) {
148  norm = LauConstants::rootPiBy2*sigmaVal_*(TMath::Erf((this->getMaxAbscissa() - meanVal_)/scale) - TMath::Erf((this->getMinAbscissa() - meanVal_)/scale));
149  }
150 
151  value /= norm;
152 
153  this->setUnNormPDFVal(value);
154 }
155 
156 void LauDPDepGaussPdf::scalePars(Double_t dpPos)
157 {
158  Int_t power = 1;
159  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
160  Double_t coeff = (*iter);
161  meanVal_ += coeff * TMath::Power(dpPos,power);
162  ++power;
163  }
164 
165  power = 1;
166  for (std::vector<Double_t>::const_iterator iter = sigmaCoeffs_.begin(); iter != sigmaCoeffs_.end(); ++iter) {
167  Double_t coeff = (*iter);
168  sigmaVal_ += coeff * TMath::Power(dpPos,power);
169  ++power;
170  }
171 }
172 
174 {
175  // Get the up to date parameter values
177 
178  // Find out the DP position
179  Double_t dpPos(0.0);
180  if ( dpAxis_ == M12 ) {
181  dpPos = kinematics->getm12Sq();
182  } else if ( dpAxis_ == M13 ) {
183  dpPos = kinematics->getm13Sq();
184  } else if ( dpAxis_ == M23 ) {
185  dpPos = kinematics->getm23Sq();
186  } else if ( dpAxis_ == MMIN ) {
187  Double_t m13Sq = kinematics->getm13Sq();
188  Double_t m23Sq = kinematics->getm23Sq();
189  dpPos = TMath::Min( m13Sq, m23Sq );
190  } else if ( dpAxis_ == MMAX ) {
191  Double_t m13Sq = kinematics->getm13Sq();
192  Double_t m23Sq = kinematics->getm23Sq();
193  dpPos = TMath::Max( m13Sq, m23Sq );
194  } else {
195  dpPos = kinematics->distanceFromDPCentre();
196  }
197 
198  // Scale the parameters according to the DP position
199  this->scalePars( dpPos );
200 
201  LauAbscissas maxPoint(3);
202  maxPoint[0] = meanVal_;
203  maxPoint[1] = kinematics->getm13Sq();
204  maxPoint[2] = kinematics->getm23Sq();
205 
206  // Calculate the PDF height for the Gaussian function.
207  if (meanVal_>this->getMaxAbscissa()) {
208  maxPoint[0] = this->getMaxAbscissa();
209  } else if (meanVal_<this->getMinAbscissa()) {
210  maxPoint[0] = this->getMinAbscissa();
211  }
212  this->calcLikelihoodInfo(maxPoint);
213 
214  Double_t height = this->getUnNormLikelihood();
215  this->setMaxHeight(height);
216 }
217 
Double_t sigmaVal_
Gaussian sigma.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:383
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:131
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:47
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:210
File containing declaration of LauDaughters class.
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:339
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:227
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 Double_t unblindValue() const =0
The unblinded value of the parameter.
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.
Double_t calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const
Calculate the third invariant mass square from the two provided (e.g. mjkSq from mijSq and mikSq) ...
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:137
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:345
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.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:55
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:118
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:43
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:59