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