ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Gamma.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Gamma.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // ------------------------------------------------------------
29 // GEANT 4 class implementation
30 // ------------------------------------------------------------
31 
32 #include <cmath>
33 #include <string.h>
34 #include "Gamma.hh"
35 
37 
39 
40 //____________________________________________________________________________
41 double MyGamma::Gamma(double z)
42 {
43  // Computation of gamma(z) for all z>0.
44  //
45  // The algorithm is based on the article by C.Lanczos [1] as denoted in
46  // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
47  //
48  // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
49  //
50  //--- Nve 14-nov-1998 UU-SAP Utrecht
51 
52  if (z<=0) return 0;
53 
54  double v = LnGamma(z);
55  return std::exp(v);
56 }
57 
58 //____________________________________________________________________________
59 double MyGamma::Gamma(double a,double x)
60 {
61  // Computation of the incomplete gamma function P(a,x)
62  //
63  // The algorithm is based on the formulas and code as denoted in
64  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
65  //
66  //--- Nve 14-nov-1998 UU-SAP Utrecht
67 
68  if (a <= 0 || x <= 0) return 0;
69 
70  if (x < (a+1)) return GamSer(a,x);
71  else return GamCf(a,x);
72 }
73 
74 //____________________________________________________________________________
75 double MyGamma::GamCf(double a,double x)
76 {
77  // Computation of the incomplete gamma function P(a,x)
78  // via its continued fraction representation.
79  //
80  // The algorithm is based on the formulas and code as denoted in
81  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
82  //
83  //--- Nve 14-nov-1998 UU-SAP Utrecht
84 
85  int itmax = 100; // Maximum number of iterations
86  double eps = 3.e-7; // Relative accuracy
87  double fpmin = 1.e-30; // Smallest double value allowed here
88 
89  if (a <= 0 || x <= 0) return 0;
90 
91  double gln = LnGamma(a);
92  double b = x+1-a;
93  double c = 1/fpmin;
94  double d = 1/b;
95  double h = d;
96  double an,del;
97  for (int i=1; i<=itmax; i++) {
98  an = double(-i)*(double(i)-a);
99  b += 2;
100  d = an*d+b;
101  if (Abs(d) < fpmin) d = fpmin;
102  c = b+an/c;
103  if (Abs(c) < fpmin) c = fpmin;
104  d = 1/d;
105  del = d*c;
106  h = h*del;
107  if (Abs(del-1) < eps) break;
108  //if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
109  }
110  double v = Exp(-x+a*Log(x)-gln)*h;
111  return (1-v);
112 }
113 
114 //____________________________________________________________________________
115 double MyGamma::GamSer(double a,double x)
116 {
117  // Computation of the incomplete gamma function P(a,x)
118  // via its series representation.
119  //
120  // The algorithm is based on the formulas and code as denoted in
121  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
122  //
123  //--- Nve 14-nov-1998 UU-SAP Utrecht
124 
125  int itmax = 100; // Maximum number of iterations
126  double eps = 3.e-7; // Relative accuracy
127 
128  if (a <= 0 || x <= 0) return 0;
129 
130  double gln = LnGamma(a);
131  double ap = a;
132  double sum = 1/a;
133  double del = sum;
134  for (int n=1; n<=itmax; n++) {
135  ap += 1;
136  del = del*x/ap;
137  sum += del;
138  if (MyGamma::Abs(del) < Abs(sum*eps)) break;
139  //if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
140  }
141  double v = sum*Exp(-x+a*Log(x)-gln);
142  return v;
143 }
144 
145 
146 double MyGamma::LnGamma(double z)
147 {
148  // Computation of ln[gamma(z)] for all z>0.
149  //
150  // The algorithm is based on the article by C.Lanczos [1] as denoted in
151  // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
152  //
153  // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
154  //
155  // The accuracy of the result is better than 2e-10.
156  //
157  //--- Nve 14-nov-1998 UU-SAP Utrecht
158 
159  if (z<=0) return 0;
160 
161  // Coefficients for the series expansion
162  double c[7] = { 2.5066282746310005, 76.18009172947146, -86.50532032941677
163  ,24.01409824083091, -1.231739572450155, 0.1208650973866179e-2
164  ,-0.5395239384953e-5};
165 
166  double x = z;
167  double y = x;
168  double tmp = x+5.5;
169  tmp = (x+0.5)*Log(tmp)-tmp;
170  double ser = 1.000000000190015;
171  for (int i=1; i<7; i++) {
172  y += 1;
173  ser += c[i]/y;
174  }
175  double v = tmp+Log(c[0]*ser/x);
176  return v;
177 }