laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepCruijffPdf.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2009 - 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 
21 #include "TMath.h"
22 #include "TSystem.h"
23 
24 #include "LauConstants.hh"
25 #include "LauDaughters.hh"
26 #include "LauDPDepCruijffPdf.hh"
27 #include "LauKinematics.hh"
28 
29 ClassImp(LauDPDepCruijffPdf)
30 
31 
32 LauDPDepCruijffPdf::LauDPDepCruijffPdf(const TString& theVarName, const std::vector<LauParameter*>& params,
33  Double_t minAbscissa, Double_t maxAbscissa,
34  const LauDaughters* daughters,
35  const std::vector<Double_t>& meanCoeffs,
36  const std::vector<Double_t>& sigmaLCoeffs,
37  const std::vector<Double_t>& sigmaRCoeffs,
38  const std::vector<Double_t>& alphaLCoeffs,
39  const std::vector<Double_t>& alphaRCoeffs,
40  DPAxis dpAxis) :
41  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
42  kinematics_(daughters ? daughters->getKinematics() : 0),
43  mean_(0),
44  sigmaL_(0),
45  sigmaR_(0),
46  alphaL_(0),
47  alphaR_(0),
48  meanVal_(0.0),
49  sigmaLVal_(0.0),
50  sigmaRVal_(0.0),
51  alphaLVal_(0.0),
52  alphaRVal_(0.0),
53  meanCoeffs_(meanCoeffs),
54  sigmaLCoeffs_(sigmaLCoeffs),
55  sigmaRCoeffs_(sigmaRCoeffs),
56  alphaLCoeffs_(alphaLCoeffs),
57  alphaRCoeffs_(alphaRCoeffs),
58  dpAxis_(dpAxis)
59 {
60  // Constructor for the Dalitz-plot dependent Cruijff PDF.
61 
62  // Check we have a valid kinematics object
63  if ( ! kinematics_ ) {
64  cerr<<"ERROR in LauDPDepCruijffPdf::LauDPDepCruijffPdf : Have not been provided with a valid DP kinematics object."<<endl;
65  gSystem->Exit(EXIT_FAILURE);
66  }
67 
68  // The parameters in params are the mean, sigmaR, sigmaL, alphaR and alphaL.
69  // The last two arguments specify the range in which the PDF is defined, and the PDF
70  // will be normalised w.r.t. these limits.
71 
72  mean_ = this->findParameter("mean");
73  sigmaR_ = this->findParameter("sigmaR");
74  sigmaL_ = this->findParameter("sigmaL");
75  alphaR_ = this->findParameter("alphaR");
76  alphaL_ = this->findParameter("alphaL");
77 
78  if ((this->nParameters() != 5) || (mean_ == 0) || (sigmaR_ == 0) || (sigmaL_ == 0) || (alphaL_ == 0) || (alphaR_ == 0)) {
79  cerr<<"ERROR in LauDPDepCruijffPdf constructor: LauDPDepCruijffPdf requires 5 parameters: \"mean\", \"sigmaL\", \"sigmaR\", \"alphaR\" and \"alphaL\"."<<endl;
80  gSystem->Exit(EXIT_FAILURE);
81  }
82 
83  // Cache the normalisation factor
84  this->calcNorm();
85 }
86 
88 {
89  // Destructor
90 }
91 
92 LauDPDepCruijffPdf::LauDPDepCruijffPdf(const LauDPDepCruijffPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa()),
93  kinematics_( other.kinematics_ ),
94  mean_(0),
95  sigmaL_(0),
96  sigmaR_(0),
97  alphaL_(0),
98  alphaR_(0),
99  meanVal_(0.0),
100  sigmaLVal_(0.0),
101  sigmaRVal_(0.0),
102  alphaLVal_(0.0),
103  alphaRVal_(0.0),
104  meanCoeffs_( other.meanCoeffs_ ),
105  sigmaLCoeffs_( other.sigmaLCoeffs_ ),
106  sigmaRCoeffs_( other.sigmaRCoeffs_ ),
107  alphaLCoeffs_( other.alphaLCoeffs_ ),
108  alphaRCoeffs_( other.alphaRCoeffs_ ),
109  dpAxis_( other.dpAxis_ )
110 {
111  // Copy constructor
112  mean_ = this->findParameter("mean");
113  sigmaR_ = this->findParameter("sigmaR");
114  sigmaL_ = this->findParameter("sigmaL");
115  alphaR_ = this->findParameter("alphaR");
116  alphaL_ = this->findParameter("alphaL");
117  this->setRandomFun(other.getRandomFun());
118  this->calcNorm();
119 }
120 
122 {
123  this->setNorm(1.0);
124 }
125 
127 {
128  // Check that the given abscissa is within the allowed range
129  if (!this->checkRange(abscissas)) {
130  gSystem->Exit(EXIT_FAILURE);
131  }
132 
133  // Get our abscissa
134  Double_t abscissa = abscissas[0];
135 
136  // Get the up to date parameter values
137  meanVal_ = mean_->value();
138  sigmaLVal_ = sigmaL_->value();
139  sigmaRVal_ = sigmaR_->value();
140  alphaLVal_ = alphaL_->value();
141  alphaRVal_ = alphaR_->value();
142 
143  // Find out the DP position
144  Double_t dpPos(0.0);
145  UInt_t nVars = this->nInputVars();
146  if ( abscissas.size() == nVars+2 ) {
147  Double_t m13Sq = abscissas[nVars];
148  Double_t m23Sq = abscissas[nVars+1];
149 
150  if ( dpAxis_ == M12 ) {
151  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
152  } else if ( dpAxis_ == M13 ) {
153  dpPos = m13Sq;
154  } else if ( dpAxis_ == M23 ) {
155  dpPos = m23Sq;
156  } else if ( dpAxis_ == MMIN ) {
157  dpPos = TMath::Min( m13Sq, m23Sq );
158  } else if ( dpAxis_ == MMAX ) {
159  dpPos = TMath::Max( m13Sq, m23Sq );
160  } else {
161  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
162  }
163  }
164 
165  // Scale the parameters according to the DP position
166  this->scalePars( dpPos );
167 
168  // Calculate the PDF value
169  Double_t value = this->currentPDFValue( abscissa );
170 
171  // Calculate the normalisation
172  IntMethod sumMethod = this->integMethod();
173  Double_t normFac = (sumMethod == GaussLegendre) ? this->integrGaussLegendre() : this->integTrapezoid();
174 
175  value /= normFac;
176 
177  this->setUnNormPDFVal(value);
178 }
179 
180 Double_t LauDPDepCruijffPdf::currentPDFValue(Double_t abscissa) const
181 {
182  Double_t arg = abscissa - meanVal_;
183  Double_t coef(0.0);
184  if (arg < 0.0){
185  if (TMath::Abs(sigmaLVal_) > 1e-30) {
186  coef = -1.0/(2.0*sigmaLVal_*sigmaLVal_ + alphaLVal_*arg*arg);
187  }
188  } else {
189  if (TMath::Abs(sigmaRVal_) > 1e-30) {
190  coef = -1.0/(2.0*sigmaRVal_*sigmaRVal_ + alphaRVal_*arg*arg);
191  }
192  }
193  return TMath::Exp(coef*arg*arg);
194 }
195 
197 {
198  if (!this->normWeightsDone()) {
199  this->getNormWeights();
200  }
201 
202  const std::vector<LauAbscissas>& norm_abscissas = this->normAbscissas();
203  const std::vector<Double_t>& norm_weights = this->normWeights();
204 
205  // Now compute the integral
206  Double_t norm(0.0);
207  for (UInt_t i = 0; i < norm_weights.size(); i++) {
208  Double_t fun = this->currentPDFValue( norm_abscissas[i][0] );
209  Double_t intFactor = 0.5 * this->getRange();
210  norm += norm_weights[i]*intFactor*fun;
211  }
212 
213  return norm;
214 }
215 
217 {
218  static Double_t norm(0.0);
219  Double_t range = this->getRange();
220 
221  if (this->nNormPoints()==1){
222 
223  Double_t abscissa = this->getMinAbscissa();
224  Double_t funAbsMin = this->currentPDFValue(abscissa);
225 
226  abscissa = this->getMaxAbscissa();
227  Double_t funAbsMax = this->currentPDFValue(abscissa);
228 
229  norm = 0.5*range*(funAbsMin+funAbsMax);
230  return norm;
231 
232  } else {
233 
234  Double_t abscVal(0.0), tnm(0.0), sum(0.0), del(0.0);
235  Int_t it(0), j(0);
236 
237  for (it=1, j=1; j< this->nNormPoints()-1; j++) {it<<=1;}
238  tnm=it;
239  del=range/tnm;
240  abscVal= this->getMinAbscissa()+ 0.5*del;
241 
242  for (sum = 0.0, j=1; j<it; j++, abscVal+=del) {
243 
244  Double_t funVal = this->currentPDFValue(abscVal);
245  sum+=funVal;
246  }
247 
248  norm = 0.5*(norm + sum*range/tnm);
249  return norm;
250  }
251 }
252 
253 void LauDPDepCruijffPdf::scalePars(Double_t dpPos)
254 {
255  Int_t power = 1;
256  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
257  Double_t coeff = (*iter);
258  meanVal_ += coeff * TMath::Power(dpPos,power);
259  ++power;
260  }
261 
262  power = 1;
263  for (std::vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
264  Double_t coeff = (*iter);
265  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
266  ++power;
267  }
268 
269  power = 1;
270  for (std::vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
271  Double_t coeff = (*iter);
272  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
273  ++power;
274  }
275 
276  power = 1;
277  for (std::vector<Double_t>::const_iterator iter = alphaLCoeffs_.begin(); iter != alphaLCoeffs_.end(); ++iter) {
278  Double_t coeff = (*iter);
279  alphaLVal_ += coeff * TMath::Power(dpPos,power);
280  ++power;
281  }
282 
283  power = 1;
284  for (std::vector<Double_t>::const_iterator iter = alphaRCoeffs_.begin(); iter != alphaRCoeffs_.end(); ++iter) {
285  Double_t coeff = (*iter);
286  alphaRVal_ += coeff * TMath::Power(dpPos,power);
287  ++power;
288  }
289 }
290 
292 {
293  // Get the up to date parameter values
294  meanVal_ = mean_->value();
295  sigmaLVal_ = sigmaL_->value();
296  sigmaRVal_ = sigmaR_->value();
297  alphaLVal_ = alphaL_->value();
298  alphaRVal_ = alphaR_->value();
299 
300  // Find out the DP position
301  Double_t dpPos(0.0);
302  if ( dpAxis_ == M12 ) {
303  dpPos = kinematics->getm12Sq();
304  } else if ( dpAxis_ == M13 ) {
305  dpPos = kinematics->getm13Sq();
306  } else if ( dpAxis_ == M23 ) {
307  dpPos = kinematics->getm23Sq();
308  } else {
309  dpPos = kinematics->distanceFromDPCentre();
310  }
311 
312  // Scale the parameters according to the DP position
313  this->scalePars( dpPos );
314 
315  // Calculate the PDF height
316 
317  LauAbscissas maxPoint(3);
318  maxPoint[0] = meanVal_;
319  maxPoint[1] = kinematics->getm13Sq();
320  maxPoint[2] = kinematics->getm23Sq();
321 
322  if (meanVal_ < this->getMinAbscissa()) {
323  maxPoint[0] = this->getMinAbscissa();
324  } else if (meanVal_ > this->getMaxAbscissa()) {
325  maxPoint[0] = this->getMaxAbscissa();
326  }
327 
328  this->calcLikelihoodInfo(maxPoint);
329  Double_t height = this->getUnNormLikelihood();
330 
331  // Multiply by a small factor to avoid problems from rounding errors
332  height *= 1.01;
333 
334  this->setMaxHeight(height);
335 }
336 
Double_t range() const
The range allowed for the parameter.
virtual Double_t integTrapezoid()
Integrate the PDF using the simple trapezoid method.
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 void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:454
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:158
virtual Bool_t normWeightsDone() const
Check whether the normalisation weights have been calculated.
Definition: LauAbsPdf.hh:532
const std::vector< Double_t > alphaLCoeffs_
Coefficients of alpha for the left Gaussian.
Class for defining a Cruijff PDF (with DP dependence).
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
virtual const std::vector< LauAbscissas > & normAbscissas() const
Retrieve the abscissa points used for normalisation.
Definition: LauAbsPdf.hh:547
Double_t currentPDFValue(Double_t abscissa) const
Current PDF value.
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:278
File containing declaration of LauDaughters class.
const LauKinematics * kinematics_
The current DP kinematics.
void scalePars(Double_t dpPos)
Scale parameters by their dependence on the DP position.
virtual void setNorm(Double_t norm)
Set the normalisation factor.
Definition: LauAbsPdf.hh:410
Double_t sigmaRVal_
Sigma of right Gaussian.
virtual Bool_t checkRange(const LauAbscissas &abscissas) const
Check that all abscissas are within their allowed ranges.
Definition: LauAbsPdf.cc:213
virtual Int_t nNormPoints() const
Retrieve the number of points to integrate over when normalising.
Definition: LauAbsPdf.hh:361
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 TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:472
LauParameter * mean_
Gaussian mean.
const std::vector< Double_t > sigmaLCoeffs_
Coefficients of sigma for the left Gaussian.
virtual ~LauDPDepCruijffPdf()
Destructor.
const std::vector< Double_t > alphaRCoeffs_
Coefficients of alpha for the right Gaussian.
DPAxis dpAxis_
The DP axis we depend on.
virtual void calcNorm()
Calculate the normalisation.
File containing declaration of LauDPDepCruijffPdf class.
Double_t meanVal_
Gaussian mean.
virtual Double_t getMaxAbscissa() const
Retrieve the maximum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:164
const std::vector< Double_t > meanCoeffs_
Coefficients of Gaussian mean.
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:416
Double_t alphaRVal_
Alpha of right Gaussian.
LauParameter * alphaR_
Alpha of right Gaussian.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:31
Double_t alphaLVal_
Alpha of left Gaussian.
virtual const std::vector< Double_t > & normWeights() const
Retrieve the weights used for normalisation.
Definition: LauAbsPdf.hh:553
Double_t getm12Sq() const
Get the m12 invariant mass square.
LauParameter * alphaL_
Alpha of left Gaussian.
LauParameter * sigmaR_
Sigma of right Gaussian.
Double_t getm13Sq() const
Get the m13 invariant mass square.
Double_t sigmaLVal_
Sigma of left Gaussian.
virtual Double_t integrGaussLegendre()
Integrate the PDF using the Gauss-Legendre method.
File containing LauConstants namespace.
virtual void getNormWeights()
Calculate the weights and abscissas used for normalisation.
Definition: LauAbsPdf.cc:569
DPAxis
Define possibilties for the DP axes.
virtual LauParameter * findParameter(const TString &parName)
Retrieve the specified parameter.
Definition: LauAbsPdf.cc:381
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:40
const std::vector< Double_t > sigmaRCoeffs_
Coefficients of sigma for the right Gaussian.
virtual IntMethod integMethod() const
Retrieve the integration method used to normalise the PDF.
Definition: LauAbsPdf.hh:373
Class for calculating 3-body kinematic quantities.
LauParameter * sigmaL_
Sigma of left Gaussian.
Double_t value() const
The value of the parameter.
virtual UInt_t nInputVars() const
Retrieve the number of abscissas.
Definition: LauAbsPdf.hh:103
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
virtual Double_t getRange() const
Retrieve the range of the (primary) abscissa.
Definition: LauAbsPdf.hh:170
IntMethod
The possible numerical intergration methods.
Definition: LauAbsPdf.hh:47
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:315
LauDPDepCruijffPdf(const TString &theVarName, const std::vector< LauParameter * > &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, const std::vector< Double_t > &alphaLCoeffs, const std::vector< Double_t > &alphaRCoeffs, DPAxis dpAxis)
Constructor.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:44