laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauBifurcatedGaussPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 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 
18 #include "TMath.h"
19 #include "TSystem.h"
20 
21 #include "LauBifurcatedGaussPdf.hh"
22 #include "LauConstants.hh"
23 
25 
26 
27 LauBifurcatedGaussPdf::LauBifurcatedGaussPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
28  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
29  mean_(0),
30  sigmaL_(0),
31  sigmaR_(0)
32 {
33  // Constructor for the bifurcated Gaussian PDF.
34  // The bifurcated Gaussian combines the left half of a
35  // Gaussian with resolution sigmaL with the right half
36  // of a Gaussian with resolution sigmaR, both having
37  // a common mean (or more correctly, peak position).
38 
39  // NB the parameters in params are the mean, sigmaL and sigmaR
40  // The last two arguments specify the range in which the PDF is defined, and the PDF
41  // will be normalised w.r.t. these limits.
42 
43  mean_ = this->findParameter("mean");
44  sigmaL_ = this->findParameter("sigmaL");
45  sigmaR_ = this->findParameter("sigmaR");
46 
47  if ((this->nParameters() != 3) || (mean_ == 0) || (sigmaL_ == 0) || (sigmaR_ == 0) ) {
48  std::cerr << "ERROR in LauBifurcatedGaussPdf constructor: LauBifurcatedGaussPdf requires 3 parameters: \"mean\", \"sigmaL\" and \"sigmaR\"." << std::endl;
49  gSystem->Exit(EXIT_FAILURE);
50  }
51 
52  // Cache the normalisation factor
53  this->calcNorm();
54 }
55 
57 {
58  // Destructor
59 }
60 
61 LauBifurcatedGaussPdf::LauBifurcatedGaussPdf(const LauBifurcatedGaussPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
62 {
63  // Copy constructor
64  this->setRandomFun(other.getRandomFun());
65  this->calcNorm();
66 }
67 
69 {
70 
71  // Check that the given abscissa is within the allowed range
72  if (!this->checkRange(abscissas)) {
73  gSystem->Exit(EXIT_FAILURE);
74  }
75 
76  // Get our abscissa
77  Double_t abscissa = abscissas[0];
78 
79  // Get the up to date parameter values
80  Double_t mean = mean_->value();
81  Double_t sigmaL = sigmaL_->value();
82  Double_t sigmaR = sigmaR_->value();
83 
84  // Evaluate the Birfucated Gaussian PDF value
85  Double_t arg = abscissa - mean;
86  Double_t coef(0.0);
87  Double_t value(0.0);
88 
89  if (arg < 0.0){
90  if (TMath::Abs(sigmaL) > 1e-30) {
91  coef = -0.5/(sigmaL*sigmaL);
92  }
93  } else {
94  if (TMath::Abs(sigmaR) > 1e-30) {
95  coef = -0.5/(sigmaR*sigmaR);
96  }
97  }
98  value = TMath::Exp(coef*arg*arg);
99 
100  // Calculate the norm
101  Double_t xscaleL = LauConstants::root2*sigmaL;
102  Double_t xscaleR = LauConstants::root2*sigmaR;
103 
104  Double_t integral(0.0);
105  Double_t norm(0.0);
106  Double_t result(0.0);
107 
108 
109  if(this->getMaxAbscissa() < mean){
110  integral = sigmaL * ( TMath::Erf((this->getMaxAbscissa() - mean)/xscaleL) - TMath::Erf((this->getMinAbscissa() - mean)/xscaleL));
111  }else if (this->getMinAbscissa() > mean){
112  integral = sigmaR * (TMath::Erf((this->getMaxAbscissa() - mean)/xscaleR) - TMath::Erf((this->getMinAbscissa() - mean)/xscaleR));
113  }else{
114  integral = sigmaR*TMath::Erf((this->getMaxAbscissa() -mean)/xscaleR) - sigmaL*TMath::Erf((this->getMinAbscissa() - mean)/xscaleL);
115  }
116 
117  norm = LauConstants::rootPiBy2*integral;
118 
119  // the result
120  result = value/norm;
121  this->setUnNormPDFVal(result);
122 
123 }
124 
126 {
127  // Nothing to do here, since it already normalized
128  this->setNorm(1.0);
129 }
130 
131 
133 {
134  if (this->heightUpToDate()) {
135  return;
136  }
137 
138  // Get the up to date parameter values
139  Double_t mean = mean_->value();
140 
141  LauAbscissas maxPoint(1);
142  maxPoint[0] = mean;
143 
144  // Calculate the PDF height for the Bifurcated Gaussian function.
145  if (mean>this->getMaxAbscissa()) {
146  maxPoint[0] = this->getMaxAbscissa();
147  } else if (mean<this->getMinAbscissa()) {
148  maxPoint[0] = this->getMinAbscissa();
149  }
150  this->calcLikelihoodInfo(maxPoint);
151 
152  Double_t height = this->getUnNormLikelihood();
153  this->setMaxHeight(height);
154 }
155 
156 
157 
158 
LauAbsRValue * sigmaL_
Sigma of left Gaussian.
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 Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
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
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
virtual ~LauBifurcatedGaussPdf()
Destructor.
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:387
LauBifurcatedGaussPdf(const TString &theVarName, const std::vector< LauAbsRValue * > &params, Double_t minAbscissa, Double_t maxAbscissa)
Constructor.
const Double_t rootPiBy2
Square root of Pi divided by two.
File containing declaration of LauBifurcatedGaussPdf class.
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
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
Class for defining a bifurcated Gaussian PDF.
LauAbsRValue * sigmaR_
Sigma of right Gaussian.
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
LauAbsRValue * mean_
Gaussian mean.
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.
virtual void calcNorm()
Calculate the normalisation.
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:233
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