laura is hosted by Hepforge, IPPP Durham
Laura++  v3r3
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 //inline static double p0(double ,double a) { return a; }
73 inline static double p1(double t,double a,double b) { return a*t+b; }
74 inline static double p2(double t,double a,double b,double c) { return p1(t,p1(t,a,b),c); }
75 inline static double p3(double t,double a,double b,double c,double d) { return p2(t,p1(t,a,b),c,d); }
76 //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); }
77 
79 {
80  // Check that the given abscissa is within the allowed range
81  if (!this->checkRange(abscissas)) {
82  gSystem->Exit(EXIT_FAILURE);
83  }
84 
85  // Get our abscissa
86  Double_t abscissa = abscissas[0];
87 
88  Double_t xmin = this->getMinAbscissa();
89  Double_t xmax = this->getMaxAbscissa();
90  Double_t x(-1.0+2.0*(abscissa-xmin)/(xmax-xmin));
91  Double_t x2(x*x);
92  Double_t sum(0) ;
93  switch ( coeffs_.size() ) {
94  case 7: sum += coeffs_[6]->unblindValue()*x*p3(x2,64,-112,56,-7);
95  case 6: sum += coeffs_[5]->unblindValue()*p3(x2,32,-48,18,-1);
96  case 5: sum += coeffs_[4]->unblindValue()*x*p2(x2,16,-20,5);
97  case 4: sum += coeffs_[3]->unblindValue()*p2(x2,8,-8,1);
98  case 3: sum += coeffs_[2]->unblindValue()*x*p1(x2,4,-3);
99  case 2: sum += coeffs_[1]->unblindValue()*p1(x2,2,-1);
100  case 1: sum += coeffs_[0]->unblindValue()*x;
101  case 0: sum +=1;
102  }
103 
104  this->setUnNormPDFVal( sum );
105 
106  // if the parameters are floating then we
107  // need to recalculate the normalisation
108  if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
109  this->calcNorm();
110  }
111 }
112 
114 {
115  Double_t xmin = this->getMinAbscissa();
116  Double_t xmax = this->getMaxAbscissa();
117 
118  Double_t norm(0) ;
119  switch( coeffs_.size() ) {
120  case 7: case 6: norm += coeffs_[5]->unblindValue()*(-1 + 18./3. - 48./5. + 32./7.);
121  case 5: case 4: norm += coeffs_[3]->unblindValue()*( 1 - 8./3. + 8./5.);
122  case 3: case 2: norm += coeffs_[1]->unblindValue()*(-1 + 2./3.);
123  case 1: case 0: norm += 1;
124  }
125  norm *= xmax-xmin;
126 
127  this->setNorm(norm);
128 }
129 
131 {
132  // TODO - this method can hopefully be improved
133  // At present it scans through the range and then increases by a 20% safety factor
134  // Maybe there's a better way?
135 
136  if (this->heightUpToDate()) {
137  return;
138  }
139 
140  // Calculate the PDF height
141  LauAbscissas maxPoint(1);
142 
143  Double_t minAbs = this->getMinAbscissa();
144  Double_t maxAbs = this->getMaxAbscissa();
145  Double_t range = maxAbs - minAbs;
146  Double_t maxHeight(0.0);
147 
148  // Just scan through the range
149  for ( Double_t point = minAbs; point <= maxAbs; point += range/1000.0 ) {
150  maxPoint[0] = point;
151  this->calcLikelihoodInfo(maxPoint);
152  Double_t heightAtPoint = this->getUnNormLikelihood();
153  if ( heightAtPoint > maxHeight ) {
154  maxHeight = heightAtPoint;
155  }
156  }
157 
158 
159  // Mutliply by 120% to be on the safe side
160  maxHeight *= 1.2;
161 
162  this->setMaxHeight(maxHeight);
163 }
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
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 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:276
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.
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