laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
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 
30 
31 
32 LauDPDepCruijffPdf::LauDPDepCruijffPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& 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 
93 {
94  this->setNorm(1.0);
95 }
96 
98 {
99  // Check that the given abscissa is within the allowed range
100  if (!this->checkRange(abscissas)) {
101  gSystem->Exit(EXIT_FAILURE);
102  }
103 
104  // Get our abscissa
105  Double_t abscissa = abscissas[0];
106 
107  // Get the up to date parameter values
108  meanVal_ = mean_->value();
109  sigmaLVal_ = sigmaL_->value();
110  sigmaRVal_ = sigmaR_->value();
111  alphaLVal_ = alphaL_->value();
112  alphaRVal_ = alphaR_->value();
113 
114  // Find out the DP position
115  Double_t dpPos(0.0);
116  UInt_t nVars = this->nInputVars();
117  if ( abscissas.size() == nVars+2 ) {
118  Double_t m13Sq = abscissas[nVars];
119  Double_t m23Sq = abscissas[nVars+1];
120 
121  if ( dpAxis_ == M12 ) {
122  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
123  } else if ( dpAxis_ == M13 ) {
124  dpPos = m13Sq;
125  } else if ( dpAxis_ == M23 ) {
126  dpPos = m23Sq;
127  } else if ( dpAxis_ == MMIN ) {
128  dpPos = TMath::Min( m13Sq, m23Sq );
129  } else if ( dpAxis_ == MMAX ) {
130  dpPos = TMath::Max( m13Sq, m23Sq );
131  } else {
132  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
133  }
134  }
135 
136  // Scale the parameters according to the DP position
137  this->scalePars( dpPos );
138 
139  // Calculate the PDF value
140  Double_t value = this->currentPDFValue( abscissa );
141 
142  // Calculate the normalisation
143  IntMethod sumMethod = this->integMethod();
144  Double_t normFac = (sumMethod == GaussLegendre) ? this->integrGaussLegendre() : this->integTrapezoid();
145 
146  value /= normFac;
147 
148  this->setUnNormPDFVal(value);
149 }
150 
151 Double_t LauDPDepCruijffPdf::currentPDFValue(Double_t abscissa) const
152 {
153  Double_t arg = abscissa - meanVal_;
154  Double_t coef(0.0);
155  if (arg < 0.0){
156  if (TMath::Abs(sigmaLVal_) > 1e-30) {
157  coef = -1.0/(2.0*sigmaLVal_*sigmaLVal_ + alphaLVal_*arg*arg);
158  }
159  } else {
160  if (TMath::Abs(sigmaRVal_) > 1e-30) {
161  coef = -1.0/(2.0*sigmaRVal_*sigmaRVal_ + alphaRVal_*arg*arg);
162  }
163  }
164  return TMath::Exp(coef*arg*arg);
165 }
166 
168 {
169  if (!this->normWeightsDone()) {
170  this->getNormWeights();
171  }
172 
173  const std::vector<LauAbscissas>& norm_abscissas = this->normAbscissas();
174  const std::vector<Double_t>& norm_weights = this->normWeights();
175 
176  // Now compute the integral
177  Double_t norm(0.0);
178  for (UInt_t i = 0; i < norm_weights.size(); i++) {
179  Double_t fun = this->currentPDFValue( norm_abscissas[i][0] );
180  Double_t intFactor = 0.5 * this->getRange();
181  norm += norm_weights[i]*intFactor*fun;
182  }
183 
184  return norm;
185 }
186 
188 {
189  static Double_t norm(0.0);
190  Double_t range = this->getRange();
191 
192  if (this->nNormPoints()==1){
193 
194  Double_t abscissa = this->getMinAbscissa();
195  Double_t funAbsMin = this->currentPDFValue(abscissa);
196 
197  abscissa = this->getMaxAbscissa();
198  Double_t funAbsMax = this->currentPDFValue(abscissa);
199 
200  norm = 0.5*range*(funAbsMin+funAbsMax);
201  return norm;
202 
203  } else {
204 
205  Double_t abscVal(0.0), tnm(0.0), sum(0.0), del(0.0);
206  Int_t it(0), j(0);
207 
208  for (it=1, j=1; j< this->nNormPoints()-1; j++) {it<<=1;}
209  tnm=it;
210  del=range/tnm;
211  abscVal= this->getMinAbscissa()+ 0.5*del;
212 
213  for (sum = 0.0, j=1; j<it; j++, abscVal+=del) {
214 
215  Double_t funVal = this->currentPDFValue(abscVal);
216  sum+=funVal;
217  }
218 
219  norm = 0.5*(norm + sum*range/tnm);
220  return norm;
221  }
222 }
223 
224 void LauDPDepCruijffPdf::scalePars(Double_t dpPos)
225 {
226  Int_t power = 1;
227  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
228  Double_t coeff = (*iter);
229  meanVal_ += coeff * TMath::Power(dpPos,power);
230  ++power;
231  }
232 
233  power = 1;
234  for (std::vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
235  Double_t coeff = (*iter);
236  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
237  ++power;
238  }
239 
240  power = 1;
241  for (std::vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
242  Double_t coeff = (*iter);
243  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
244  ++power;
245  }
246 
247  power = 1;
248  for (std::vector<Double_t>::const_iterator iter = alphaLCoeffs_.begin(); iter != alphaLCoeffs_.end(); ++iter) {
249  Double_t coeff = (*iter);
250  alphaLVal_ += coeff * TMath::Power(dpPos,power);
251  ++power;
252  }
253 
254  power = 1;
255  for (std::vector<Double_t>::const_iterator iter = alphaRCoeffs_.begin(); iter != alphaRCoeffs_.end(); ++iter) {
256  Double_t coeff = (*iter);
257  alphaRVal_ += coeff * TMath::Power(dpPos,power);
258  ++power;
259  }
260 }
261 
263 {
264  // Get the up to date parameter values
265  meanVal_ = mean_->value();
266  sigmaLVal_ = sigmaL_->value();
267  sigmaRVal_ = sigmaR_->value();
268  alphaLVal_ = alphaL_->value();
269  alphaRVal_ = alphaR_->value();
270 
271  // Find out the DP position
272  Double_t dpPos(0.0);
273  if ( dpAxis_ == M12 ) {
274  dpPos = kinematics->getm12Sq();
275  } else if ( dpAxis_ == M13 ) {
276  dpPos = kinematics->getm13Sq();
277  } else if ( dpAxis_ == M23 ) {
278  dpPos = kinematics->getm23Sq();
279  } else {
280  dpPos = kinematics->distanceFromDPCentre();
281  }
282 
283  // Scale the parameters according to the DP position
284  this->scalePars( dpPos );
285 
286  // Calculate the PDF height
287 
288  LauAbscissas maxPoint(3);
289  maxPoint[0] = meanVal_;
290  maxPoint[1] = kinematics->getm13Sq();
291  maxPoint[2] = kinematics->getm23Sq();
292 
293  if (meanVal_ < this->getMinAbscissa()) {
294  maxPoint[0] = this->getMinAbscissa();
295  } else if (meanVal_ > this->getMaxAbscissa()) {
296  maxPoint[0] = this->getMaxAbscissa();
297  }
298 
299  this->calcLikelihoodInfo(maxPoint);
300  Double_t height = this->getUnNormLikelihood();
301 
302  // Multiply by a small factor to avoid problems from rounding errors
303  height *= 1.01;
304 
305  this->setMaxHeight(height);
306 }
307 
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) ...
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 normWeightsDone() const
Check whether the normalisation weights have been calculated.
Definition: LauAbsPdf.hh:447
ClassImp(LauAbsCoeffSet)
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
LauAbsRValue * alphaL_
Alpha of left Gaussian.
virtual const std::vector< LauAbscissas > & normAbscissas() const
Retrieve the abscissa points used for normalisation.
Definition: LauAbsPdf.hh:462
Double_t currentPDFValue(Double_t abscissa) const
Current PDF value.
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
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:325
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:282
Double_t getm23Sq() const
Get the m23 invariant mass square.
LauAbsRValue * sigmaR_
Sigma of right Gaussian.
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 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:123
const std::vector< Double_t > meanCoeffs_
Coefficients of Gaussian mean.
LauAbsRValue * mean_
Gaussian mean.
virtual void setMaxHeight(Double_t maxHeight)
Set the maximum height.
Definition: LauAbsPdf.hh:331
Double_t alphaRVal_
Alpha of right Gaussian.
Double_t alphaLVal_
Alpha of left Gaussian.
virtual const std::vector< Double_t > & normWeights() const
Retrieve the weights used for normalisation.
Definition: LauAbsPdf.hh:468
Double_t getm12Sq() const
Get the m12 invariant mass square.
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:470
DPAxis
Define possibilties for the DP axes.
virtual Double_t value() const =0
Return the value of the parameter.
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
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:294
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 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:129
IntMethod
The possible numerical intergration methods.
Definition: LauAbsPdf.hh:48
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:29
LauAbsRValue * alphaR_
Alpha of right Gaussian.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45