laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepBifurGaussPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2007 - 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 using std::vector;
21 
22 #include "TMath.h"
23 #include "TRandom.h"
24 #include "TSystem.h"
25 
26 #include "Lau1DHistPdf.hh"
27 #include "LauConstants.hh"
28 #include "LauComplex.hh"
29 #include "LauDaughters.hh"
30 #include "LauFitDataTree.hh"
31 #include "LauKinematics.hh"
32 #include "LauDPDepBifurGaussPdf.hh"
33 
35 
36 
37 LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf(const TString& theVarName, const vector<LauAbsRValue*>& params,
38  Double_t minAbscissa, Double_t maxAbscissa,
39  const LauDaughters* daughters,
40  const vector<Double_t>& meanCoeffs,
41  const vector<Double_t>& sigmaLCoeffs,
42  const vector<Double_t>& sigmaRCoeffs,
43  DPAxis dpAxis) :
44  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
45  kinematics_(daughters ? daughters->getKinematics() : 0),
46  mean_(0),
47  sigmaL_(0),
48  sigmaR_(0),
49  meanVal_(0.0),
50  sigmaLVal_(0.0),
51  sigmaRVal_(0.0),
52  meanCoeffs_(meanCoeffs),
53  sigmaLCoeffs_(sigmaLCoeffs),
54  sigmaRCoeffs_(sigmaRCoeffs),
55  dpAxis_(dpAxis),
56  scaleMethod_(poly)
57 {
58  // Check we have a valid kinematics object
59  if ( ! kinematics_ ) {
60  cerr<<"ERROR in LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf : Have not been provided with a valid DP kinematics object."<<endl;
61  gSystem->Exit(EXIT_FAILURE);
62  }
63 
64  // The parameters are:
65  // - the mean, the sigmaL and the sigmaR of the Bifurcated Gaussian
66  //
67  // The next two arguments specify the range in which the PDF is defined,
68  // and the PDF will be normalised w.r.t. these limits.
69  //
70  // The final three argument define whether the BF Gaussian parameters should be scaled or not
71 
72  mean_ = this->findParameter("mean");
73  sigmaL_ = this->findParameter("sigmaL");
74  sigmaR_ = this->findParameter("sigmaR");
75 
76  if ((this->nParameters() != 3) || (mean_ == 0) || (sigmaL_ == 0) || (sigmaR_ == 0) ) {
77  cerr<<"ERROR in LauDPDepBifurGaussPdf constructor: LauDPDepBifurGaussPdf requires 3 parameters:"
78  <<" \"mean\", \"sigmaL\" and \"sigmaR\"."<<endl;
79  gSystem->Exit(EXIT_FAILURE);
80  }
81 
82  // Cache the normalisation factor
83  this->calcNorm();
84 }
85 
87 {
88  // Destructor
89 }
90 
91 LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf(const LauDPDepBifurGaussPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa()),
92  kinematics_(other.kinematics_),
93  mean_(other.mean_),
94  sigmaL_(other.sigmaL_),
95  sigmaR_(other.sigmaR_),
96  meanCoeffs_(other.meanCoeffs_),
97  sigmaLCoeffs_(other.sigmaLCoeffs_),
98  sigmaRCoeffs_(other.sigmaRCoeffs_),
99  scaleMethod_(other.poly)
100 {
101  // Copy constructor
102  this->setRandomFun(other.getRandomFun());
103  this->calcNorm();
104 }
105 
107 {
108  // Check that the given abscissas are within the allowed range
109  if (!this->checkRange(abscissas)) {
110  gSystem->Exit(EXIT_FAILURE);
111  }
112 
113  // Get our abscissa
114  Double_t abscissa = abscissas[0];
115 
116  // Get the up to date parameter values
117  meanVal_ = mean_->value();
118  sigmaLVal_ = sigmaL_->value();
119  sigmaRVal_ = sigmaR_->value();
120 
121  // Find out the DP position
122  Double_t dpPos(0.0);
123  UInt_t nVars = this->nInputVars();
124  if ( abscissas.size() == nVars+2 ) {
125  Double_t m13Sq = abscissas[nVars];
126  Double_t m23Sq = abscissas[nVars+1];
127 
128  if ( dpAxis_ == M12 ) {
129  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
130  } else if ( dpAxis_ == M13 ) {
131  dpPos = m13Sq;
132  } else if ( dpAxis_ == M23 ) {
133  dpPos = m23Sq;
134  } else if ( dpAxis_ == MMIN ) {
135  dpPos = TMath::Min( m13Sq, m23Sq );
136  } else if ( dpAxis_ == MMAX ) {
137  dpPos = TMath::Max( m13Sq, m23Sq );
138  } else {
139  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
140  }
141  }
142 
143  // Scale the gaussian parameters by the dpPos (if appropriate)
144  ScaleMethod scale = this->scaleMethod();
145  if ( scale==poly ) {
146  this->scalePars_poly(dpPos);
147  } else if ( scale==polyNegPower ) {
148  this->scalePars_polyNegPower(dpPos);
149  } else {
150  cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
151  gSystem->Exit(EXIT_FAILURE);
152  }
153 
154  // Evaluate the Birfucated Gaussian PDF value
155  Double_t arg = abscissa - meanVal_;
156  Double_t coef(0.0);
157  Double_t value(0.0);
158 
159  if (arg < 0.0) {
160  if (TMath::Abs(sigmaLVal_) > 1e-30) {
161  coef = -0.5/(sigmaLVal_*sigmaLVal_);
162  }
163  } else {
164  if (TMath::Abs(sigmaRVal_) > 1e-30) {
165  coef = -0.5/(sigmaRVal_*sigmaRVal_);
166  }
167  }
168  value = TMath::Exp(coef*arg*arg);
169 
170  // Calculate the norm
171  Double_t xscaleL = LauConstants::root2*sigmaLVal_;
172  Double_t xscaleR = LauConstants::root2*sigmaRVal_;
173 
174  Double_t integral(0.0);
175  Double_t norm(0.0);
176  Double_t result(0.0);
177 
178 
179  if ( this->getMaxAbscissa() < meanVal_ ) {
180  integral = sigmaLVal_ * ( TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleL) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL));
181  } else if ( this->getMinAbscissa() > meanVal_ ) {
182  integral = sigmaRVal_ * (TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleR) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleR));
183  } else {
184  integral = sigmaRVal_*TMath::Erf((this->getMaxAbscissa() -meanVal_)/xscaleR) - sigmaLVal_*TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL);
185  }
186 
187  norm = LauConstants::rootPiBy2*integral;
188 
189  // the result
190  result = value/norm;
191  this->setUnNormPDFVal(result);
192 }
193 
195 {
196  Int_t power = 1;
197  for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
198  Double_t coeff = (*iter);
199  meanVal_ += coeff * TMath::Power(dpPos,power);
200  ++power;
201  }
202 
203  power = 1;
204  for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
205  Double_t coeff = (*iter);
206  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
207  ++power;
208  }
209 
210  power = 1;
211  for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
212  Double_t coeff = (*iter);
213  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
214  ++power;
215  }
216 }
217 
219 {
220  Int_t power = -1;
221  for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
222  Double_t coeff = (*iter);
223  meanVal_ += coeff * TMath::Power(dpPos,power);
224  --power;
225  }
226 
227  power = -1;
228  for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
229  Double_t coeff = (*iter);
230  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
231  --power;
232  }
233 
234  power = -1;
235  for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
236  Double_t coeff = (*iter);
237  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
238  --power;
239  }
240 }
241 
243 {
244  this->setNorm(1.0);
245 }
246 
248 {
249  // Get the up to date parameter values
250  meanVal_ = mean_->value();
251  sigmaLVal_ = sigmaL_->value();
252  sigmaRVal_ = sigmaR_->value();
253 
254  // Scale the gaussian parameters by the dpCentreDist (if appropriate)
255  Double_t dpCentreDist = kinematics->distanceFromDPCentre();
256  ScaleMethod scale = this->scaleMethod();
257  if ( scale==poly ) {
258  this->scalePars_poly(dpCentreDist);
259  } else if ( scale==polyNegPower ) {
260  this->scalePars_polyNegPower(dpCentreDist);
261  } else {
262  cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
263  gSystem->Exit(EXIT_FAILURE);
264  }
265 
266  // Calculate the PDF height for the Bifurcated Gaussian function.
267  LauAbscissas maxPoint(3);
268  maxPoint[0] = meanVal_;
269  maxPoint[1] = kinematics->getm13Sq();
270  maxPoint[2] = kinematics->getm23Sq();
271  if ( meanVal_ > this->getMaxAbscissa() ) {
272  maxPoint[0] = this->getMaxAbscissa();
273  } else if ( meanVal_ < this->getMinAbscissa() ) {
274  maxPoint[0] = this->getMinAbscissa();
275  }
276  this->calcLikelihoodInfo(maxPoint);
277  Double_t height = this->getUnNormLikelihood();
278 
279  // Multiply by a small factor to avoid problems from rounding errors
280  height *= 1.001;
281 
282  this->setMaxHeight(height);
283 }
284 
Double_t sigmaLVal_
Left 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 ~LauDPDepBifurGaussPdf()
Destructor.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:369
File containing declaration of LauFitDataTree class.
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:117
ClassImp(LauAbsCoeffSet)
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
const std::vector< Double_t > meanCoeffs_
Coefficients of Gaussian mean.
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
File containing declaration of LauDaughters class.
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.
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:325
void scalePars_poly(Double_t perEventDist)
Scale the Gaussian parameters with polynomial method.
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:213
const LauKinematics * kinematics_
The current DP kinematics.
Double_t getm23Sq() const
Get the m23 invariant mass square.
void scalePars_polyNegPower(Double_t perEventDist)
Scale the Gaussian parameters with negative power polynomial method.
DPAxis dpAxis_
The DP axis we depend on.
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:387
const Double_t rootPiBy2
Square root of Pi divided by two.
LauDPDepBifurGaussPdf(const TString &theVarName, const std::vector< LauAbsRValue * > &params, Double_t minAbscissa, Double_t maxAbscissa, const LauDaughters *daughters, const std::vector< Double_t > &meanCoeffs, const std::vector< Double_t > &sigmaLCoeffs, const std::vector< Double_t > &sigmaRCoeffs, DPAxis dpAxis)
Constructor.
const Double_t root2
Square root of two.
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:123
File containing declaration of LauComplex class.
File containing declaration of LauDPDepBifurGaussPdf class.
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
ScaleMethod scaleMethod() const
Retrieve scaling method information.
ScaleMethod
Define possibilties for the scaling method.
const std::vector< Double_t > sigmaLCoeffs_
Coefficients of left Gaussian sigma.
Double_t sigmaRVal_
Right Gaussian sigma.
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:41
const std::vector< Double_t > sigmaRCoeffs_
Coefficients of right Gaussian sigma.
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
virtual UInt_t nInputVars() const
Retrieve the number of abscissas.
Definition: LauAbsPdf.hh:104
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:233
virtual void calcNorm()
Calculate the normalisation.
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
Double_t meanVal_
Gaussian mean.
Class for defining a Bifurcated Gaussian PDF (DP dependent).
File containing declaration of Lau1DHistPdf class.