laura is hosted by Hepforge, IPPP Durham
Laura++  v3r5
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPDepCruijffPdf.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2009 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include <iostream>
30 #include <vector>
31 using std::cout;
32 using std::cerr;
33 using std::endl;
34 
35 #include "TMath.h"
36 #include "TSystem.h"
37 
38 #include "LauConstants.hh"
39 #include "LauDaughters.hh"
40 #include "LauDPDepCruijffPdf.hh"
41 #include "LauKinematics.hh"
42 
44 
45 
46 LauDPDepCruijffPdf::LauDPDepCruijffPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params,
47  Double_t minAbscissa, Double_t maxAbscissa,
48  const LauDaughters* daughters,
49  const std::vector<Double_t>& meanCoeffs,
50  const std::vector<Double_t>& sigmaLCoeffs,
51  const std::vector<Double_t>& sigmaRCoeffs,
52  const std::vector<Double_t>& alphaLCoeffs,
53  const std::vector<Double_t>& alphaRCoeffs,
54  DPAxis dpAxis) :
55  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
56  kinematics_(daughters ? daughters->getKinematics() : 0),
57  mean_(0),
58  sigmaL_(0),
59  sigmaR_(0),
60  alphaL_(0),
61  alphaR_(0),
62  meanVal_(0.0),
63  sigmaLVal_(0.0),
64  sigmaRVal_(0.0),
65  alphaLVal_(0.0),
66  alphaRVal_(0.0),
67  meanCoeffs_(meanCoeffs),
68  sigmaLCoeffs_(sigmaLCoeffs),
69  sigmaRCoeffs_(sigmaRCoeffs),
70  alphaLCoeffs_(alphaLCoeffs),
71  alphaRCoeffs_(alphaRCoeffs),
72  dpAxis_(dpAxis)
73 {
74  // Constructor for the Dalitz-plot dependent Cruijff PDF.
75 
76  // Check we have a valid kinematics object
77  if ( ! kinematics_ ) {
78  cerr<<"ERROR in LauDPDepCruijffPdf::LauDPDepCruijffPdf : Have not been provided with a valid DP kinematics object."<<endl;
79  gSystem->Exit(EXIT_FAILURE);
80  }
81 
82  // The parameters in params are the mean, sigmaR, sigmaL, alphaR and alphaL.
83  // The last two arguments specify the range in which the PDF is defined, and the PDF
84  // will be normalised w.r.t. these limits.
85 
86  mean_ = this->findParameter("mean");
87  sigmaR_ = this->findParameter("sigmaR");
88  sigmaL_ = this->findParameter("sigmaL");
89  alphaR_ = this->findParameter("alphaR");
90  alphaL_ = this->findParameter("alphaL");
91 
92  if ((this->nParameters() != 5) || (mean_ == 0) || (sigmaR_ == 0) || (sigmaL_ == 0) || (alphaL_ == 0) || (alphaR_ == 0)) {
93  cerr<<"ERROR in LauDPDepCruijffPdf constructor: LauDPDepCruijffPdf requires 5 parameters: \"mean\", \"sigmaL\", \"sigmaR\", \"alphaR\" and \"alphaL\"."<<endl;
94  gSystem->Exit(EXIT_FAILURE);
95  }
96 
97  // Cache the normalisation factor
98  this->calcNorm();
99 }
100 
102 {
103  // Destructor
104 }
105 
107 {
108  this->setNorm(1.0);
109 }
110 
112 {
113  // Check that the given abscissa is within the allowed range
114  if (!this->checkRange(abscissas)) {
115  gSystem->Exit(EXIT_FAILURE);
116  }
117 
118  // Get our abscissa
119  Double_t abscissa = abscissas[0];
120 
121  // Get the up to date parameter values
127 
128  // Find out the DP position
129  Double_t dpPos(0.0);
130  UInt_t nVars = this->nInputVars();
131  if ( abscissas.size() == nVars+2 ) {
132  Double_t m13Sq = abscissas[nVars];
133  Double_t m23Sq = abscissas[nVars+1];
134 
135  if ( dpAxis_ == M12 ) {
136  dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
137  } else if ( dpAxis_ == M13 ) {
138  dpPos = m13Sq;
139  } else if ( dpAxis_ == M23 ) {
140  dpPos = m23Sq;
141  } else if ( dpAxis_ == MMIN ) {
142  dpPos = TMath::Min( m13Sq, m23Sq );
143  } else if ( dpAxis_ == MMAX ) {
144  dpPos = TMath::Max( m13Sq, m23Sq );
145  } else {
146  dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
147  }
148  }
149 
150  // Scale the parameters according to the DP position
151  this->scalePars( dpPos );
152 
153  // Calculate the PDF value
154  Double_t value = this->currentPDFValue( abscissa );
155 
156  // Calculate the normalisation
157  IntMethod sumMethod = this->integMethod();
158  Double_t normFac = (sumMethod == GaussLegendre) ? this->integrGaussLegendre() : this->integTrapezoid();
159 
160  value /= normFac;
161 
162  this->setUnNormPDFVal(value);
163 }
164 
165 Double_t LauDPDepCruijffPdf::currentPDFValue(Double_t abscissa) const
166 {
167  Double_t arg = abscissa - meanVal_;
168  Double_t coef(0.0);
169  if (arg < 0.0){
170  if (TMath::Abs(sigmaLVal_) > 1e-30) {
171  coef = -1.0/(2.0*sigmaLVal_*sigmaLVal_ + alphaLVal_*arg*arg);
172  }
173  } else {
174  if (TMath::Abs(sigmaRVal_) > 1e-30) {
175  coef = -1.0/(2.0*sigmaRVal_*sigmaRVal_ + alphaRVal_*arg*arg);
176  }
177  }
178  return TMath::Exp(coef*arg*arg);
179 }
180 
182 {
183  if (!this->normWeightsDone()) {
184  this->getNormWeights();
185  }
186 
187  const std::vector<LauAbscissas>& norm_abscissas = this->normAbscissas();
188  const std::vector<Double_t>& norm_weights = this->normWeights();
189 
190  // Now compute the integral
191  Double_t norm(0.0);
192  for (UInt_t i = 0; i < norm_weights.size(); i++) {
193  Double_t fun = this->currentPDFValue( norm_abscissas[i][0] );
194  Double_t intFactor = 0.5 * this->getRange();
195  norm += norm_weights[i]*intFactor*fun;
196  }
197 
198  return norm;
199 }
200 
202 {
203  static Double_t norm(0.0);
204  Double_t range = this->getRange();
205 
206  if (this->nNormPoints()==1){
207 
208  Double_t abscissa = this->getMinAbscissa();
209  Double_t funAbsMin = this->currentPDFValue(abscissa);
210 
211  abscissa = this->getMaxAbscissa();
212  Double_t funAbsMax = this->currentPDFValue(abscissa);
213 
214  norm = 0.5*range*(funAbsMin+funAbsMax);
215  return norm;
216 
217  } else {
218 
219  Double_t abscVal(0.0), tnm(0.0), sum(0.0), del(0.0);
220  Int_t it(0), j(0);
221 
222  for (it=1, j=1; j< this->nNormPoints()-1; j++) {it<<=1;}
223  tnm=it;
224  del=range/tnm;
225  abscVal= this->getMinAbscissa()+ 0.5*del;
226 
227  for (sum = 0.0, j=1; j<it; j++, abscVal+=del) {
228 
229  Double_t funVal = this->currentPDFValue(abscVal);
230  sum+=funVal;
231  }
232 
233  norm = 0.5*(norm + sum*range/tnm);
234  return norm;
235  }
236 }
237 
238 void LauDPDepCruijffPdf::scalePars(Double_t dpPos)
239 {
240  Int_t power = 1;
241  for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
242  Double_t coeff = (*iter);
243  meanVal_ += coeff * TMath::Power(dpPos,power);
244  ++power;
245  }
246 
247  power = 1;
248  for (std::vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
249  Double_t coeff = (*iter);
250  sigmaLVal_ += coeff * TMath::Power(dpPos,power);
251  ++power;
252  }
253 
254  power = 1;
255  for (std::vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
256  Double_t coeff = (*iter);
257  sigmaRVal_ += coeff * TMath::Power(dpPos,power);
258  ++power;
259  }
260 
261  power = 1;
262  for (std::vector<Double_t>::const_iterator iter = alphaLCoeffs_.begin(); iter != alphaLCoeffs_.end(); ++iter) {
263  Double_t coeff = (*iter);
264  alphaLVal_ += coeff * TMath::Power(dpPos,power);
265  ++power;
266  }
267 
268  power = 1;
269  for (std::vector<Double_t>::const_iterator iter = alphaRCoeffs_.begin(); iter != alphaRCoeffs_.end(); ++iter) {
270  Double_t coeff = (*iter);
271  alphaRVal_ += coeff * TMath::Power(dpPos,power);
272  ++power;
273  }
274 }
275 
277 {
278  // Get the up to date parameter values
284 
285  // Find out the DP position
286  Double_t dpPos(0.0);
287  if ( dpAxis_ == M12 ) {
288  dpPos = kinematics->getm12Sq();
289  } else if ( dpAxis_ == M13 ) {
290  dpPos = kinematics->getm13Sq();
291  } else if ( dpAxis_ == M23 ) {
292  dpPos = kinematics->getm23Sq();
293  } else {
294  dpPos = kinematics->distanceFromDPCentre();
295  }
296 
297  // Scale the parameters according to the DP position
298  this->scalePars( dpPos );
299 
300  // Calculate the PDF height
301 
302  LauAbscissas maxPoint(3);
303  maxPoint[0] = meanVal_;
304  maxPoint[1] = kinematics->getm13Sq();
305  maxPoint[2] = kinematics->getm23Sq();
306 
307  if (meanVal_ < this->getMinAbscissa()) {
308  maxPoint[0] = this->getMinAbscissa();
309  } else if (meanVal_ > this->getMaxAbscissa()) {
310  maxPoint[0] = this->getMaxAbscissa();
311  }
312 
313  this->calcLikelihoodInfo(maxPoint);
314  Double_t height = this->getUnNormLikelihood();
315 
316  // Multiply by a small factor to avoid problems from rounding errors
317  height *= 1.01;
318 
319  this->setMaxHeight(height);
320 }
321 
Double_t range() const
The range allowed for the parameter.
virtual Double_t integTrapezoid()
Integrate the PDF using the simple trapezoid method.
LauAbsRValue * sigmaL_
Sigma of left Gaussian.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:383
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:131
virtual Bool_t normWeightsDone() const
Check whether the normalisation weights have been calculated.
Definition: LauAbsPdf.hh:461
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:47
LauAbsRValue * alphaL_
Alpha of left Gaussian.
virtual const std::vector< LauAbscissas > & normAbscissas() const
Retrieve the abscissa points used for normalisation.
Definition: LauAbsPdf.hh:476
Double_t currentPDFValue(Double_t abscissa) const
Current PDF value.
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:210
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:339
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:227
virtual Int_t nNormPoints() const
Retrieve the number of points to integrate over when normalising.
Definition: LauAbsPdf.hh:296
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.
virtual Double_t unblindValue() const =0
The unblinded value of the parameter.
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.
Double_t calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const
Calculate the third invariant mass square from the two provided (e.g. mjkSq from mijSq and mikSq) ...
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:137
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:345
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:482
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:484
DPAxis
Define possibilties for the DP axes.
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:55
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:308
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:118
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:143
IntMethod
The possible numerical intergration methods.
Definition: LauAbsPdf.hh:62
Pure abstract base class for defining a parameter containing an R value.
Definition: LauAbsRValue.hh:43
LauAbsRValue * alphaR_
Alpha of right Gaussian.
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:59