laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1p1
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 
82 LauCrystalBallPdf::LauCrystalBallPdf(const LauCrystalBallPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
83 {
84  // Copy constructor
85  this->setRandomFun(other.getRandomFun());
86  this->calcNorm();
87 }
88 
90 {
91  // Check that the given abscissa is within the allowed range
92  if (!this->checkRange(abscissas)) {
93  gSystem->Exit(EXIT_FAILURE);
94  }
95 
96  // Get our abscissa
97  Double_t abscissa = abscissas[0];
98 
99  // Get the up to date parameter values
100  Double_t mean = mean_->value();
101  Double_t sigma = sigma_->value();
102  Double_t alpha = alpha_->value();
103  Double_t n = n_->value();
104 
105  Double_t result(0.0);
106  Double_t t = (abscissa - mean)/sigma;
107  if (alpha < 0.0) {
108  t = -t;
109  }
110 
111  Double_t absAlpha = TMath::Abs(alpha);
112 
113  if (t >= -absAlpha) {
114  result = TMath::Exp(-0.5*t*t);
115  } else {
116  Double_t a = TMath::Power(n/absAlpha,n)*TMath::Exp(-0.5*absAlpha*absAlpha);
117  Double_t b = n/absAlpha - absAlpha;
118 
119  result = a/TMath::Power(b - t, n);
120  }
121 
122  this->setUnNormPDFVal(result);
123 
124  // if the parameters are floating then we
125  // need to recalculate the normalisation
126  if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
127  this->calcNorm();
128  }
129 
130 }
131 
133 {
134  // Get the up to date parameter values
135  Double_t mean = mean_->value();
136  Double_t sigma = sigma_->value();
137  Double_t alpha = alpha_->value();
138  Double_t n = n_->value();
139 
140  Double_t result = 0.0;
141  Bool_t useLog = kFALSE;
142 
143  if ( TMath::Abs(n-1.0) < 1.0e-05 ) {
144  useLog = kTRUE;
145  }
146 
147  Double_t sig = TMath::Abs(sigma);
148 
149  Double_t tmin = (this->getMinAbscissa() - mean)/sig;
150  Double_t tmax = (this->getMaxAbscissa() - mean)/sig;
151 
152  if (alpha < 0) {
153  Double_t tmp = tmin;
154  tmin = -tmax;
155  tmax = -tmp;
156  }
157 
158  Double_t absAlpha = TMath::Abs(alpha);
159 
160  if ( tmin >= -absAlpha ) {
161  result += sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
162  - approxErf( tmin / LauConstants::root2 ) );
163  } else if ( tmax <= -absAlpha ) {
164  Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha);
165  Double_t b = n/absAlpha - absAlpha;
166 
167  if ( useLog == kTRUE ) {
168  result += a*sig*( TMath::Log(b - tmin) - TMath::Log(b - tmax) );
169  } else {
170  result += a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
171  - 1.0/(TMath::Power( b - tmax, n - 1.0)) );
172  }
173  } else {
174  Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha );
175  Double_t b = n/absAlpha - absAlpha;
176 
177  Double_t term1 = 0.0;
178  if ( useLog == kTRUE )
179  term1 = a*sig*( TMath::Log(b - tmin) - TMath::Log(n / absAlpha));
180  else
181  term1 = a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
182  - 1.0/(TMath::Power( n/absAlpha, n - 1.0)) );
183 
184  Double_t term2 = sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
185  - this->approxErf( -absAlpha / LauConstants::root2 ) );
186 
187  result += term1 + term2;
188  }
189  this->setNorm(result);
190 }
191 
192 void LauCrystalBallPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
193 {
194  if (this->heightUpToDate()) {
195  return;
196  }
197 
198  // Get the up to date parameter values
199  Double_t mean = mean_->value();
200 
201  // The Crystall Ball function is a Gaussian with an exponentially decaying tail
202  // Therefore, calculate the PDF height for the Gaussian function.
203  LauAbscissas abscissa(1);
204  abscissa[0] = mean;
205  this->calcLikelihoodInfo(abscissa);
206 
207  Double_t height = this->getUnNormLikelihood();
208  this->setMaxHeight(height);
209 }
210 
211 Double_t LauCrystalBallPdf::approxErf(Double_t arg) const
212 {
213  static const Double_t erflim = 5.0;
214  if ( arg > erflim ) {
215  return 1.0;
216  }
217  if ( arg < -erflim ) {
218  return -1.0;
219  }
220 
221  return TMath::Erf(arg);
222 }
223 
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
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:387
const Double_t rootPiBy2
Square root of Pi divided by two.
LauCrystalBallPdf(const TString &theVarName, const std::vector< LauAbsRValue * > &params, Double_t minAbscissa, Double_t maxAbscissa)
Constructor.
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:270
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.
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
std::vector< Double_t > LauAbscissas
The type used for containing multiple abscissa values.
Definition: LauAbsPdf.hh:45