laura is hosted by Hepforge, IPPP Durham
Laura++  v3r3
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCrystalBallPdf.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 /*****************************************************************************
16  * Class based on RooFit/RooCBShape. *
17  * Original copyright given below. *
18  *****************************************************************************
19  * Authors: *
20  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
21  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
22  * *
23  * Copyright (c) 2000-2005, Regents of the University of California *
24  * and Stanford University. All rights reserved. *
25  * *
26  * Redistribution and use in source and binary forms, *
27  * with or without modification, are permitted according to the terms *
28  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
29  *****************************************************************************/
30 
31 #include <iostream>
32 #include <vector>
33 using std::cout;
34 using std::cerr;
35 using std::endl;
36 using std::vector;
37 
38 #include "TMath.h"
39 #include "TSystem.h"
40 
41 #include "LauConstants.hh"
42 #include "LauCrystalBallPdf.hh"
43 
45 
46 
47 LauCrystalBallPdf::LauCrystalBallPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
48  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
49  mean_(0),
50  sigma_(0),
51  alpha_(0),
52  n_(0)
53 {
54  // Constructor for the Crystal Ball PDF, which is a gaussian and a decaying tail
55  // smoothly matched up. The tail goes as a 1/x^n
56  //
57  // The parameters in params are the mean and the sigma (half the width) of the gaussian,
58  // the distance from the mean in which the gaussian and the tail are matched up (which
59  // can be negative or positive), and the power "n" for the tail.
60  // The last two arguments specify the range in which the PDF is defined, and the PDF
61  // will be normalised w.r.t. these limits.
62 
63  mean_ = this->findParameter("mean");
64  sigma_ = this->findParameter("sigma");
65  alpha_ = this->findParameter("alpha");
66  n_ = this->findParameter("order");
67 
68  if ((this->nParameters() != 4) || (mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
69  cerr<<"ERROR in LauCrystalBallPdf constructor: LauCrystalBallPdf requires 4 parameters: \"mean\", \"sigma\", \"alpha\" and \"order\"."<<endl;
70  gSystem->Exit(EXIT_FAILURE);
71  }
72 
73  // Cache the normalisation factor
74  this->calcNorm();
75 }
76 
78 {
79  // Destructor
80 }
81 
83 {
84  // Check that the given abscissa is within the allowed range
85  if (!this->checkRange(abscissas)) {
86  gSystem->Exit(EXIT_FAILURE);
87  }
88 
89  // Get our abscissa
90  Double_t abscissa = abscissas[0];
91 
92  // Get the up to date parameter values
93  Double_t mean = mean_->unblindValue();
94  Double_t sigma = sigma_->unblindValue();
95  Double_t alpha = alpha_->unblindValue();
96  Double_t n = n_->unblindValue();
97 
98  Double_t result(0.0);
99  Double_t t = (abscissa - mean)/sigma;
100  if (alpha < 0.0) {
101  t = -t;
102  }
103 
104  Double_t absAlpha = TMath::Abs(alpha);
105 
106  if (t >= -absAlpha) {
107  result = TMath::Exp(-0.5*t*t);
108  } else {
109  Double_t a = TMath::Power(n/absAlpha,n)*TMath::Exp(-0.5*absAlpha*absAlpha);
110  Double_t b = n/absAlpha - absAlpha;
111 
112  result = a/TMath::Power(b - t, n);
113  }
114 
115  this->setUnNormPDFVal(result);
116 
117  // if the parameters are floating then we
118  // need to recalculate the normalisation
119  if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
120  this->calcNorm();
121  }
122 
123 }
124 
126 {
127  // Get the up to date parameter values
128  Double_t mean = mean_->unblindValue();
129  Double_t sigma = sigma_->unblindValue();
130  Double_t alpha = alpha_->unblindValue();
131  Double_t n = n_->unblindValue();
132 
133  Double_t result = 0.0;
134  Bool_t useLog = kFALSE;
135 
136  if ( TMath::Abs(n-1.0) < 1.0e-05 ) {
137  useLog = kTRUE;
138  }
139 
140  Double_t sig = TMath::Abs(sigma);
141 
142  Double_t tmin = (this->getMinAbscissa() - mean)/sig;
143  Double_t tmax = (this->getMaxAbscissa() - mean)/sig;
144 
145  if (alpha < 0) {
146  Double_t tmp = tmin;
147  tmin = -tmax;
148  tmax = -tmp;
149  }
150 
151  Double_t absAlpha = TMath::Abs(alpha);
152 
153  if ( tmin >= -absAlpha ) {
154  result += sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
155  - approxErf( tmin / LauConstants::root2 ) );
156  } else if ( tmax <= -absAlpha ) {
157  Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha);
158  Double_t b = n/absAlpha - absAlpha;
159 
160  if ( useLog == kTRUE ) {
161  result += a*sig*( TMath::Log(b - tmin) - TMath::Log(b - tmax) );
162  } else {
163  result += a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
164  - 1.0/(TMath::Power( b - tmax, n - 1.0)) );
165  }
166  } else {
167  Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha );
168  Double_t b = n/absAlpha - absAlpha;
169 
170  Double_t term1 = 0.0;
171  if ( useLog == kTRUE )
172  term1 = a*sig*( TMath::Log(b - tmin) - TMath::Log(n / absAlpha));
173  else
174  term1 = a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
175  - 1.0/(TMath::Power( n/absAlpha, n - 1.0)) );
176 
177  Double_t term2 = sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
178  - this->approxErf( -absAlpha / LauConstants::root2 ) );
179 
180  result += term1 + term2;
181  }
182  this->setNorm(result);
183 }
184 
185 void LauCrystalBallPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
186 {
187  if (this->heightUpToDate()) {
188  return;
189  }
190 
191  // Get the up to date parameter values
192  Double_t mean = mean_->unblindValue();
193 
194  // The Crystall Ball function is a Gaussian with an exponentially decaying tail
195  // Therefore, calculate the PDF height for the Gaussian function.
196  LauAbscissas abscissa(1);
197  abscissa[0] = mean;
198  this->calcLikelihoodInfo(abscissa);
199 
200  Double_t height = this->getUnNormLikelihood();
201  this->setMaxHeight(height);
202 }
203 
204 Double_t LauCrystalBallPdf::approxErf(Double_t arg) const
205 {
206  static const Double_t erflim = 5.0;
207  if ( arg > erflim ) {
208  return 1.0;
209  }
210  if ( arg < -erflim ) {
211  return -1.0;
212  }
213 
214  return TMath::Erf(arg);
215 }
216 
virtual void calcNorm()
Calculate the normalisation.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:369
virtual ~LauCrystalBallPdf()
Destructor.
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:117
Double_t approxErf(Double_t arg) const
Calculate the approximate error function of argument.
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
File containing declaration of LauCrystalBallPdf class.
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 Bool_t withinNormCalc() const
Check whether the calcNorm method is running.
Definition: LauAbsPdf.hh:423
const Double_t rootPiBy2
Square root of Pi divided by two.
const Double_t root2
Square root of two.
Class for defining a Crystal Ball PDF.
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
virtual Bool_t withinGeneration() const
Check whether the generate method is running.
Definition: LauAbsPdf.hh:435
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:276
File containing LauConstants namespace.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
Class for calculating 3-body kinematic quantities.
virtual void calcLikelihoodInfo(const LauAbscissas &abscissas)
Calculate the likelihood (and intermediate info) for a given abscissa.
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