laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauChebychevPdf.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 /*****************************************************************************
16  * Class based on RooFit/RooChebychev. *
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 
34 #include "TMath.h"
35 #include "TSystem.h"
36 
37 #include "LauChebychevPdf.hh"
38 
40 
41 LauChebychevPdf::LauChebychevPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
42  LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
43  coeffs_(params.size(),0)
44 {
45  // Constructor for the Chebychev Polynomial PDF.
46  //
47  // The parameters in params are the coefficients of the polynomial
48  // (polynomial can be anything from 0 to 7 orders).
49  // The last two arguments specify the range in which the PDF is defined, and the PDF
50  // will be normalised w.r.t. these limits.
51 
52  if ( (this->nParameters() > 7) ) {
53  std::cerr << "ERROR in LauChebychevPdf constructor: Too many coeffs - can only cope with order <=7 polynomial." << std::endl;
54  gSystem->Exit(EXIT_FAILURE);
55  }
56 
57  for ( UInt_t i(0); i < this->nParameters(); ++i ) {
58  TString name = "c";
59  name += i+1;
60  coeffs_[i] = this->findParameter( name );
61  }
62 
63  // Cache the normalisation factor
64  this->calcNorm();
65 }
66 
68 {
69  // Destructor
70 }
71 
72 LauChebychevPdf::LauChebychevPdf(const LauChebychevPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
73 {
74  // Copy constructor
75  this->setRandomFun(other.getRandomFun());
76  this->calcNorm();
77 }
78 
79 inline static double p0(double ,double a) { return a; }
80 inline static double p1(double t,double a,double b) { return a*t+b; }
81 inline static double p2(double t,double a,double b,double c) { return p1(t,p1(t,a,b),c); }
82 inline static double p3(double t,double a,double b,double c,double d) { return p2(t,p1(t,a,b),c,d); }
83 inline static double p4(double t,double a,double b,double c,double d,double e) { return p3(t,p1(t,a,b),c,d,e); }
84 
86 {
87  // Check that the given abscissa is within the allowed range
88  if (!this->checkRange(abscissas)) {
89  gSystem->Exit(EXIT_FAILURE);
90  }
91 
92  // Get our abscissa
93  Double_t abscissa = abscissas[0];
94 
95  Double_t xmin = this->getMinAbscissa();
96  Double_t xmax = this->getMaxAbscissa();
97  Double_t x(-1.0+2.0*(abscissa-xmin)/(xmax-xmin));
98  Double_t x2(x*x);
99  Double_t sum(0) ;
100  switch ( coeffs_.size() ) {
101  case 7: sum += coeffs_[6]->value()*x*p3(x2,64,-112,56,-7);
102  case 6: sum += coeffs_[5]->value()*p3(x2,32,-48,18,-1);
103  case 5: sum += coeffs_[4]->value()*x*p2(x2,16,-20,5);
104  case 4: sum += coeffs_[3]->value()*p2(x2,8,-8,1);
105  case 3: sum += coeffs_[2]->value()*x*p1(x2,4,-3);
106  case 2: sum += coeffs_[1]->value()*p1(x2,2,-1);
107  case 1: sum += coeffs_[0]->value()*x;
108  case 0: sum +=1;
109  }
110 
111  this->setUnNormPDFVal( sum );
112 
113  // if the parameters are floating then we
114  // need to recalculate the normalisation
115  if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
116  this->calcNorm();
117  }
118 }
119 
121 {
122  Double_t xmin = this->getMinAbscissa();
123  Double_t xmax = this->getMaxAbscissa();
124 
125  Double_t norm(0) ;
126  switch( coeffs_.size() ) {
127  case 7: case 6: norm += coeffs_[5]->value()*(-1 + 18./3. - 48./5. + 32./7.);
128  case 5: case 4: norm += coeffs_[3]->value()*( 1 - 8./3. + 8./5.);
129  case 3: case 2: norm += coeffs_[1]->value()*(-1 + 2./3.);
130  case 1: case 0: norm += 1;
131  }
132  norm *= xmax-xmin;
133 
134  this->setNorm(norm);
135 }
136 
138 {
139  // TODO - this method can hopefully be improved
140  // At present it scans through the range and then increases by a 20% safety factor
141  // Maybe there's a better way?
142 
143  if (this->heightUpToDate()) {
144  return;
145  }
146 
147  // Calculate the PDF height
148  LauAbscissas maxPoint(1);
149 
150  Double_t minAbs = this->getMinAbscissa();
151  Double_t maxAbs = this->getMaxAbscissa();
152  Double_t range = maxAbs - minAbs;
153  Double_t maxHeight(0.0);
154 
155  // Just scan through the range
156  for ( Double_t point = minAbs; point <= maxAbs; point += range/1000.0 ) {
157  maxPoint[0] = point;
158  this->calcLikelihoodInfo(maxPoint);
159  Double_t heightAtPoint = this->getUnNormLikelihood();
160  if ( heightAtPoint > maxHeight ) {
161  maxHeight = heightAtPoint;
162  }
163  }
164 
165 
166  // Mutliply by 120% to be on the safe side
167  maxHeight *= 1.2;
168 
169  this->setMaxHeight(maxHeight);
170 }
Double_t range() const
The range allowed for the parameter.
virtual void setUnNormPDFVal(Double_t unNormPDFVal)
Set the unnormalised likelihood.
Definition: LauAbsPdf.hh:369
static double p1(double t, double a, double b)
virtual Double_t getMinAbscissa() const
Retrieve the minimum value of the (primary) abscissa.
Definition: LauAbsPdf.hh:117
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
Class for defining a Chebychev Polynomial (1st kind) PDF.
ClassImp(LauAbsCoeffSet)
virtual void calcPDFHeight(const LauKinematics *kinematics)
Calculate the PDF height.
const TString & name() const
The parameter name.
std::vector< LauAbsRValue * > coeffs_
Coefficients of polynomial.
virtual Double_t getUnNormLikelihood() const
Retrieve the unnormalised likelihood value.
Definition: LauAbsPdf.hh:196
static double p0(double, double a)
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 Bool_t withinNormCalc() const
Check whether the calcNorm method is running.
Definition: LauAbsPdf.hh:423
static double p4(double t, double a, double b, double c, double d, double e)
virtual TRandom * getRandomFun() const
Retrieve the random function used for MC generation.
Definition: LauAbsPdf.hh:387
static double p3(double t, double a, double b, double c, double d)
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
static double p2(double t, double a, double b, double c)
virtual Bool_t withinGeneration() const
Check whether the generate method is running.
Definition: LauAbsPdf.hh:435
virtual void calcNorm()
Calculate the normalisation.
virtual Bool_t cachePDF() const
Check if the PDF is to be cached.
Definition: LauAbsPdf.hh:270
virtual ~LauChebychevPdf()
Destructor.
Class for defining the abstract interface for PDF classes.
Definition: LauAbsPdf.hh:41
Class for calculating 3-body kinematic quantities.
File containing declaration of LauChebychevPdf class.
LauChebychevPdf(const TString &theVarName, const std::vector< LauAbsRValue * > &params, Double_t minAbscissa, Double_t maxAbscissa)
Constructor.
virtual void setRandomFun(TRandom *randomFun)
Set the random function used for toy MC generation.
Definition: LauAbsPdf.hh:233
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