laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauArgusPdf.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 "LauArgusPdf.hh"
22 #include "LauConstants.hh"
23 
25 
26 
27 LauArgusPdf::LauArgusPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
28  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
29  xi_(0),
30  m0_(0)
31 {
32  // Constructor for the ARGUS PDF.
33  //
34  // The parameters in params are the shape, xi, and the end-point of the curve.
35  // The last two arguments specify the range in which the PDF is defined, and the PDF
36  // will be normalised w.r.t. these limits.
37 
38  xi_ = this->findParameter("xi");
39  m0_ = this->findParameter("m0");
40 
41  if ((this->nParameters() != 2) || (xi_ == 0) || (m0_ == 0)) {
42  std::cerr << "ERROR in LauArgusPdf constructor: LauArgusPdf requires 2 parameters: argus shape parameter, \"xi\", and end-point, \"m0\"." << std::endl;
43  gSystem->Exit(EXIT_FAILURE);
44  }
45 
46  // Cache the normalisation factor.
47  this->calcNorm();
48 }
49 
51 {
52  // Destructor
53 }
54 
55 LauArgusPdf::LauArgusPdf(const LauArgusPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
56 {
57  // Copy constructor
58  this->setRandomFun(other.getRandomFun());
59  this->calcNorm();
60 }
61 
63 {
64  // Check that the given abscissa is within the allowed range
65  if (!this->checkRange(abscissas)) {
66  gSystem->Exit(EXIT_FAILURE);
67  }
68 
69  // Get our abscissa
70  Double_t abscissa = abscissas[0];
71 
72  // Get the up to date parameter values
73  Double_t xi = xi_->value();
74  Double_t m0 = m0_->value();
75 
76  // Calculate the value of the ARGUS function for the given value of the abscissa.
77  Double_t x = abscissa/m0;
78 
79  Double_t term = 1.0 - x*x;
80  if (term < 0.0) {term = 0.0;} // In case |x| > 1.0 (which should happen rarely).
81 
82  Double_t value = abscissa*TMath::Sqrt(term)*TMath::Exp(-xi*term);
83  this->setUnNormPDFVal(value);
84 
85  // if the parameters are floating then we
86  // need to recalculate the normalisation
87  if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
88  this->calcNorm();
89  }
90 }
91 
93 {
94  // Calculate the PDF normalisation and cache it
95 
96  // Get the up to date parameter values
97  Double_t xi = xi_->value();
98  Double_t m0 = m0_->value();
99 
100  // Since the PDF is 0 above m0 by definition need to check whether m0 is within the range, above it or below it
101  Double_t min = (this->getMinAbscissa() < m0) ? this->getMinAbscissa() : m0;
102  Double_t max = (this->getMaxAbscissa() < m0) ? this->getMaxAbscissa() : m0;
103 
104  // Define variables equivalent to "term" in calcLikelihoodInfo above but at the min and max points
105  Double_t termMin = 1.0 - (min/m0)*(min/m0);
106  Double_t termMax = 1.0 - (max/m0)*(max/m0);
107 
108  // Calculate the various terms in the integrals
109  Double_t norm1 = TMath::Sqrt(termMax)*TMath::Exp(-xi*termMax) - TMath::Sqrt(termMin)*TMath::Exp(-xi*termMin);
110  Double_t norm2 = LauConstants::rootPi/(2.0*TMath::Sqrt(xi)) * ( TMath::Erf(TMath::Sqrt(xi*termMax)) - TMath::Erf(TMath::Sqrt(xi*termMin)) );
111 
112  // Combine them and set the normalisation
113  Double_t norm = m0*m0*(norm1 - norm2)/(2.0*xi);
114 
115  this->setNorm(norm);
116 }
117 
118 void LauArgusPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
119 {
120  if (this->heightUpToDate()) {
121  return;
122  }
123 
124  // Calculate the PDF height of the ARGUS function.
125  // Get the up to date parameter values
126  Double_t xi = xi_->value();
127  Double_t m0 = m0_->value();
128 
129  // First make sure that the limits are not larger than the end-point.
130  // (Btw, use the logarithmic derivative to derive this formula)
131  Double_t term = xi*xi + 1.0;
132  Double_t x = TMath::Sqrt((TMath::Sqrt(term) - 1.0 + xi)/(2.0*xi));
133  x = (x*m0 >= this->getMinAbscissa()) ? x*m0 : this->getMinAbscissa();
134 
135  LauAbscissas abscissa(1);
136  abscissa[0] = x;
137  this->calcLikelihoodInfo(abscissa);
138 
139  Double_t height = this->getUnNormLikelihood();
140  this->setMaxHeight(height);
141 }
142 
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
LauArgusPdf(const TString &theVarName, const std::vector< LauAbsRValue * > &params, Double_t minAbscissa, Double_t maxAbscissa)
Constructor.
Definition: LauArgusPdf.cc:27
ClassImp(LauAbsCoeffSet)
virtual ~LauArgusPdf()
Destructor.
Definition: LauArgusPdf.cc:50
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
virtual void calcNorm()
Calculate the normalisation.
Definition: LauArgusPdf.cc:92
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 Bool_t withinNormCalc() const
Check whether the calcNorm method is running.
Definition: LauAbsPdf.hh:423
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:387
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Definition: LauArgusPdf.cc:118
File containing declaration of LauArgusPdf class.
LauAbsRValue * m0_
Endpoint of curve.
Definition: LauArgusPdf.hh:76
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:123
const Double_t rootPi
Square root of Pi.
Definition: LauConstants.hh:91
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
virtual Bool_t withinGeneration() const
Check whether the generate method is running.
Definition: LauAbsPdf.hh:435
LauAbsRValue * xi_
Shape of curve.
Definition: LauArgusPdf.hh:73
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:270
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.
Class for defining an ARGUS PDF.
Definition: LauArgusPdf.hh:34
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
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
Definition: LauArgusPdf.cc:62
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45