laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
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 
92 {
93  // Check that the given abscissas are within the allowed range
94  if (!this->checkRange(abscissas)) {
95  gSystem->Exit(EXIT_FAILURE);
96  }
97 
98  // Get our abscissa
99  Double_t abscissa = abscissas[0];
100 
101  // Get the up to date parameter values
102  meanVal_ = mean_->unblindValue();
103  sigmaLVal_ = sigmaL_->unblindValue();
104  sigmaRVal_ = sigmaR_->unblindValue();
105 
106  // Find out the DP position
107  Double_t dpPos(0.0);
108  UInt_t nVars = this->nInputVars();
109  if ( abscissas.size() == nVars+2 ) {
110  Double_t m13Sq = abscissas[nVars];
111  Double_t m23Sq = abscissas[nVars+1];
112 
113  if ( dpAxis_ == M12 ) {
114  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
115  } else if ( dpAxis_ == M13 ) {
116  dpPos = m13Sq;
117  } else if ( dpAxis_ == M23 ) {
118  dpPos = m23Sq;
119  } else if ( dpAxis_ == MMIN ) {
120  dpPos = TMath::Min( m13Sq, m23Sq );
121  } else if ( dpAxis_ == MMAX ) {
122  dpPos = TMath::Max( m13Sq, m23Sq );
123  } else {
124  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
125  }
126  }
127 
128  // Scale the gaussian parameters by the dpPos (if appropriate)
129  ScaleMethod scale = this->scaleMethod();
130  if ( scale==poly ) {
131  this->scalePars_poly(dpPos);
132  } else if ( scale==polyNegPower ) {
133  this->scalePars_polyNegPower(dpPos);
134  } else {
135  cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
136  gSystem->Exit(EXIT_FAILURE);
137  }
138 
139  // Evaluate the Birfucated Gaussian PDF value
140  Double_t arg = abscissa - meanVal_;
141  Double_t coef(0.0);
142  Double_t value(0.0);
143 
144  if (arg < 0.0) {
145  if (TMath::Abs(sigmaLVal_) > 1e-30) {
146  coef = -0.5/(sigmaLVal_*sigmaLVal_);
147  }
148  } else {
149  if (TMath::Abs(sigmaRVal_) > 1e-30) {
150  coef = -0.5/(sigmaRVal_*sigmaRVal_);
151  }
152  }
153  value = TMath::Exp(coef*arg*arg);
154 
155  // Calculate the norm
156  Double_t xscaleL = LauConstants::root2*sigmaLVal_;
157  Double_t xscaleR = LauConstants::root2*sigmaRVal_;
158 
159  Double_t integral(0.0);
160  Double_t norm(0.0);
161  Double_t result(0.0);
162 
163 
164  if ( this->getMaxAbscissa() < meanVal_ ) {
165  integral = sigmaLVal_ * ( TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleL) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL));
166  } else if ( this->getMinAbscissa() > meanVal_ ) {
167  integral = sigmaRVal_ * (TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleR) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleR));
168  } else {
169  integral = sigmaRVal_*TMath::Erf((this->getMaxAbscissa() -meanVal_)/xscaleR) - sigmaLVal_*TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL);
170  }
171 
172  norm = LauConstants::rootPiBy2*integral;
173 
174  // the result
175  result = value/norm;
176  this->setUnNormPDFVal(result);
177 }
178 
180 {
181  Int_t power = 1;
182  for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
183  Double_t coeff = (*iter);
184  meanVal_ += coeff * TMath::Power(dpPos,power);
185  ++power;
186  }
187 
188  power = 1;
189  for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
190  Double_t coeff = (*iter);
191  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
192  ++power;
193  }
194 
195  power = 1;
196  for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
197  Double_t coeff = (*iter);
198  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
199  ++power;
200  }
201 }
202 
204 {
205  Int_t power = -1;
206  for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
207  Double_t coeff = (*iter);
208  meanVal_ += coeff * TMath::Power(dpPos,power);
209  --power;
210  }
211 
212  power = -1;
213  for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
214  Double_t coeff = (*iter);
215  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
216  --power;
217  }
218 
219  power = -1;
220  for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
221  Double_t coeff = (*iter);
222  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
223  --power;
224  }
225 }
226 
228 {
229  this->setNorm(1.0);
230 }
231 
233 {
234  // Get the up to date parameter values
235  meanVal_ = mean_->unblindValue();
236  sigmaLVal_ = sigmaL_->unblindValue();
237  sigmaRVal_ = sigmaR_->unblindValue();
238 
239  // Scale the gaussian parameters by the dpCentreDist (if appropriate)
240  Double_t dpCentreDist = kinematics->distanceFromDPCentre();
241  ScaleMethod scale = this->scaleMethod();
242  if ( scale==poly ) {
243  this->scalePars_poly(dpCentreDist);
244  } else if ( scale==polyNegPower ) {
245  this->scalePars_polyNegPower(dpCentreDist);
246  } else {
247  cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
248  gSystem->Exit(EXIT_FAILURE);
249  }
250 
251  // Calculate the PDF height for the Bifurcated Gaussian function.
252  LauAbscissas maxPoint(3);
253  maxPoint[0] = meanVal_;
254  maxPoint[1] = kinematics->getm13Sq();
255  maxPoint[2] = kinematics->getm23Sq();
256  if ( meanVal_ > this->getMaxAbscissa() ) {
257  maxPoint[0] = this->getMaxAbscissa();
258  } else if ( meanVal_ < this->getMinAbscissa() ) {
259  maxPoint[0] = this->getMinAbscissa();
260  }
261  this->calcLikelihoodInfo(maxPoint);
262  Double_t height = this->getUnNormLikelihood();
263 
264  // Multiply by a small factor to avoid problems from rounding errors
265  height *= 1.001;
266 
267  this->setMaxHeight(height);
268 }
269 
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.
const Double_t rootPiBy2
Square root of Pi divided by two.
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 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.