laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauIntegrals.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 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 #include <iostream>
16 using std::cout;
17 using std::endl;
18 
19 #include "TMath.h"
20 
21 #include "LauConstants.hh"
22 #include "LauIntegrals.hh"
23 
24 
25 ClassImp(LauIntegrals)
26 
27 
28 LauIntegrals::LauIntegrals(Double_t weightsPrecision) :
29  weightsPrecision_( weightsPrecision )
30 {
31 }
32 
34 {
35 }
36 
37 void LauIntegrals::calcGaussLegendreWeights(const Int_t numPoints, std::vector<Double_t>& abscissas, std::vector<Double_t>& weights)
38 {
39  // Calculate the Gauss-Legendre weights that will be used for the
40  // simple Gaussian integration method given the number of points
41 
42  abscissas.clear();
43  weights.clear();
44  abscissas.resize(numPoints);
45  weights.resize(numPoints);
46 
47  Int_t m = (numPoints+1)/2;
48  Double_t dnumPoints(numPoints);
49  Double_t dnumPointsPlusHalf = static_cast<Double_t>(numPoints + 0.5);
50  Int_t i(0), j(0);
51  Double_t p1(0.0), p2(0.0), p3(0.0), pp(0.0), z(0.0), zSq(0.0), z1(0.0), di(0.0);
52 
53  for (i = 1; i <= m; i++){
54 
55  di = static_cast<Double_t>(i);
56  z = TMath::Cos(LauConstants::pi*(di - 0.25)/dnumPointsPlusHalf);
57  zSq = z*z;
58 
59  // Starting with the above approximation for the ith root, we enter the
60  // main loop of refinement by Newton's method
61  do{
62  p1 = 1.0;
63  p2 = 0.0;
64 
65  // Calculate the Legendre polynomial at z - recurrence relation
66  for (j = 1; j <= numPoints; j++){
67  p3 = p2;
68  p2 = p1;
69  p1 = (( 2.0*j - 1.0) * z * p2 - (j - 1.0)*p3)/static_cast<Double_t>(j);
70  }
71  // p1 = Legendre polynomial. Compute its derivative, pp.
72 
73  pp = dnumPoints * (z*p1 - p2)/(zSq - 1.0);
74  z1 = z;
75  z = z1 - p1/pp; // Newton's method
76  } while (TMath::Abs(z-z1) > weightsPrecision_);
77 
78  // Scale the root to the desired interval
79  // Remember that the vector entries start with 0, hence i-1 (where first i value is 1)
80  abscissas[i-1] = z;
81  abscissas[numPoints-i] = abscissas[i-1]; // Symmetric abscissa
82  weights[i-1] = 2.0/((1.0 - zSq)*pp*pp);
83  weights[numPoints-i] = weights[i-1]; // Symmetric weight
84 
85  }
86 }
87 
88 void LauIntegrals::calcGaussHermiteWeights(const Int_t numPoints, std::vector<Double_t>& abscissas, std::vector<Double_t>& weights)
89 {
90  // Calculate the Gauss-Hermite weights that will be used for the
91  // simple Gaussian integration method for a function weighted by an exp(-x**2) term
92  // These weights are common to all integration intervals - therefore they
93  // should only be calculated once - when the constructor is called, for example.
94 
95  abscissas.clear();
96  weights.clear();
97  abscissas.resize(numPoints); weights.resize(numPoints);
98 
99  Int_t m = (numPoints+1)/2;
100  Double_t dnumPoints(numPoints);
101  Int_t i(0), j(0), its(0);
102  Double_t p1(0.0), p2(0.0), p3(0.0), pp(0.0), z(0.0), z1(0.0);
103  Double_t numPointsTerm(2*numPoints+1);
104 
105  // The roots are symmetric about the origin. Therefore, we only need to find half of them
106  // Loop over the desired roots.
107  for (i = 1; i <= m; i++) {
108 
109  if (i == 1) {
110  // Initial guess for largest root
111  z = TMath::Sqrt(numPointsTerm) - 1.85575*TMath::Power(numPointsTerm, -0.16667);
112  } else if (i == 2) {
113  // Initial guess for second largest root
114  z -= 1.14*TMath::Power(dnumPoints, 0.426)/z;
115  } else if (i == 3) {
116  // Initial guess for third largest root
117  z = 1.86*z - 0.86*abscissas[0];
118  } else if (i == 4) {
119  // Initial guess for fourth largest root
120  z = 1.91*z - 0.91*abscissas[1];
121  } else {
122  // Initial guess for other roots
123  z = 2.0*z - abscissas[i-3];
124  }
125 
126  for (its = 1; its <= 10; its++) {
128  p2 = 0.0;
129 
130  for (j = 1; j <= numPoints; j++) {
131  Double_t dj(j);
132  p3 = p2;
133  p2 = p1;
134  p1 = z*TMath::Sqrt(2.0/dj)*p2 - TMath::Sqrt((dj-1.0)/dj)*p3;
135  }
136  // p1 is now the desired Hermite polynomial. Compute its derivative, pp
137  pp = TMath::Sqrt(2.0*dnumPoints)*p2;
138  z1 = z;
139  if (pp > 1e-10) {z = z1 - p1/pp;}
140  if (TMath::Abs(z - z1) < weightsPrecision_) break;
141  }
142 
143  // Scale the root to the desired interval
144  // Remember that the vector entries start with 0, hence i-1 (where first i value is 1)
145  abscissas[i-1] = z;
146  abscissas[numPoints-i] = -z; // Symmetric abscissa
147  weights[i-1] = 2.0/(pp*pp);
148  weights[numPoints-i] = weights[i-1]; // Symmetric weight
149 
150  }
151 }
static double p1(double t, double a, double b)
Class for performing numerical integration routines.
Definition: LauIntegrals.hh:30
void calcGaussHermiteWeights(const Int_t numPoints, std::vector< Double_t > &abscissas, std::vector< Double_t > &weights)
Calculate the Gauss-Hermite weights.
Definition: LauIntegrals.cc:88
void calcGaussLegendreWeights(const Int_t numPoints, std::vector< Double_t > &abscissas, std::vector< Double_t > &weights)
Calculate the Gauss-Legendre weights.
Definition: LauIntegrals.cc:37
static double p3(double t, double a, double b, double c, double d)
const Double_t pi
Pi.
Definition: LauConstants.hh:89
static double p2(double t, double a, double b, double c)
Double_t weightsPrecision_
The precision to which the weights should be calculated.
Definition: LauIntegrals.hh:60
File containing LauConstants namespace.
File containing declaration of LauIntegrals class.
virtual ~LauIntegrals()
Destructor.
Definition: LauIntegrals.cc:33
const Double_t pim4
One over Pi to the one-fourth.
Definition: LauConstants.hh:99