| Geant4 Cross Reference |
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 // G4Pow class implemenation
27 //
28 // Author: Vladimir Ivanchenko, 23.05.2009
29 // Revisions:
30 // - 23.03.2018, M.Novak: increased accuracy of A13 on the most critical
31 // [1/4,4] interval by introducing a denser(0.25) grid
32 // --------------------------------------------------------------------
33
34 #include "G4Pow.hh"
35 #include "G4Threading.hh"
36
37 G4Pow* G4Pow::fpInstance = nullptr;
38
39 // -------------------------------------------------------------------
40
41 G4Pow* G4Pow::GetInstance()
42 {
43 if(fpInstance == nullptr)
44 {
45 static G4Pow geant4pow;
46 fpInstance = &geant4pow;
47 }
48 return fpInstance;
49 }
50
51 // -------------------------------------------------------------------
52
53 G4Pow::G4Pow()
54 {
55 #ifdef G4MULTITHREADED
56 if(G4Threading::IsWorkerThread())
57 {
58 G4Exception("G4Pow::G4Pow()", "InvalidSetup", FatalException,
59 "Attempt to instantiate G4Pow in worker thread!");
60 }
61 #endif
62 const G4int maxZ = 512;
63 const G4int maxZfact = 170;
64 const G4int numLowA = 17;
65
66 maxA = -0.6 + maxZ;
67 maxLowA = 4.0;
68 maxA2 = 1.25 + max2 * 0.2;
69 maxAexp = -0.76 + maxZfact * 0.5;
70
71 ener.resize(max2 + 1, 1.0);
72 logen.resize(max2 + 1, 0.0);
73 lz2.resize(max2 + 1, 0.0);
74 pz13.resize(maxZ, 0.0);
75 lowa13.resize(numLowA, 0.0);
76 lz.resize(maxZ, 0.0);
77 fexp.resize(maxZfact, 0.0);
78 fact.resize(maxZfact, 0.0);
79 logfact.resize(maxZ, 0.0);
80
81 G4double f = 1.0;
82 G4double logf = 0.0;
83 fact[0] = 1.0;
84 fexp[0] = 1.0;
85
86 for(G4int i = 1; i <= max2; ++i)
87 {
88 ener[i] = powN(500., i);
89 logen[i] = G4Log(ener[i]);
90 lz2[i] = G4Log(1.0 + i * 0.2);
91 }
92
93 for(G4int i = 1; i < maxZ; ++i)
94 {
95 auto x = G4double(i);
96 pz13[i] = std::pow(x, onethird);
97 lz[i] = G4Log(x);
98 if(i < maxZfact)
99 {
100 f *= x;
101 fact[i] = f;
102 fexp[i] = G4Exp(0.5 * i);
103 }
104 logf += lz[i];
105 logfact[i] = logf;
106 }
107
108 for(G4int i = 4; i < numLowA; ++i)
109 {
110 lowa13[i] = std::pow(0.25 * i, onethird);
111 }
112 }
113
114 // -------------------------------------------------------------------
115
116 G4double G4Pow::A13(G4double A) const
117 {
118 G4double res = 0.;
119 if(A > 0.)
120 {
121 const bool invert = (A < 1.);
122 const G4double a = invert ? 1. / A : A;
123 res = (a < maxLowA) ? A13Low(a, invert) : A13High(a, invert);
124 }
125 return res;
126 }
127
128 // -------------------------------------------------------------------
129
130 G4double G4Pow::A13High(const G4double a, const bool invert) const
131 {
132 G4double res;
133 if(a < maxA)
134 {
135 const auto i = static_cast<G4int>(a + 0.5);
136 const G4double x = (a / i - 1.) * onethird;
137 res = pz13[i] * (1. + x - x * x * (1. - 1.666667 * x));
138 }
139 else
140 {
141 res = G4Exp(G4Log(a) * onethird);
142 }
143 res = invert ? 1. / res : res;
144 return res;
145 }
146
147 // -------------------------------------------------------------------
148
149 G4double G4Pow::A13Low(const G4double a, const G4bool invert) const
150 {
151 G4double res;
152 const auto i = static_cast<G4int>(4. * (a + 0.125));
153 const G4double y = 0.25 * i;
154 const G4double x = (a / y - 1.) * onethird;
155 res = lowa13[i] * (1. + x - x * x * (1. - 1.666667 * x));
156 res = invert ? 1. / res : res;
157 return res;
158 }
159
160 // -------------------------------------------------------------------
161
162 G4double G4Pow::powN(G4double x, G4int n) const
163 {
164 if(0.0 == x)
165 {
166 return 0.0;
167 }
168 if(std::abs(n) > 8)
169 {
170 return std::pow(x, G4double(n));
171 }
172 G4double res = 1.0;
173 if(n >= 0)
174 {
175 for(G4int i = 0; i < n; ++i)
176 {
177 res *= x;
178 }
179 }
180 else if(n < 0)
181 {
182 G4double y = 1.0 / x;
183 G4int nn = -n;
184 for(G4int i = 0; i < nn; ++i)
185 {
186 res *= y;
187 }
188 }
189 return res;
190 }
191