ECCE @ EIC Software
Reference for
ECCE @ EIC
simulation and reconstruction software on GitHub
Home page
Related Pages
Modules
Namespaces
Classes
Files
External Links
File List
File Members
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
36
MyGamma::MyGamma
(){}
37
38
MyGamma::~MyGamma
(){}
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
}
geant4
tree
geant4-10.6-release
source
parameterisations
gflash
src
Gamma.cc
Built by
Jin Huang
. updated:
Wed Jun 29 2022 17:25:23
using
1.8.2 with
ECCE GitHub integration