| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4Pow class implemenation <<
27 // 26 //
28 // Author: Vladimir Ivanchenko, 23.05.2009 << 27 // -------------------------------------------------------------------
29 // Revisions: << 28 //
30 // - 23.03.2018, M.Novak: increased accuracy o << 29 // GEANT4 Class file
31 // [1/4,4] interval by << 30 //
32 // ------------------------------------------- << 31 //
>> 32 // File name: G4Pow
>> 33 //
>> 34 // Author: Vladimir Ivanchenko
>> 35 //
>> 36 // Creation date: 23.05.2009
>> 37 //
>> 38 // Modifications:
>> 39 // 08.01.2011 V.Ivanchenko extended maxZ from 256 to 512
>> 40 // 02.05.2013 V.Ivanchenko added expA and logX methods,
>> 41 // revised A13, logA, powZ, powA to improved accuracy
>> 42 // 23.03.2018 M.Novak increased accuracy of A13 on the most critical
>> 43 // [1/4,4] interval by introducing a denser(0.25) grid
>> 44 //
>> 45 // -------------------------------------------------------------------
33 46
34 #include "G4Pow.hh" 47 #include "G4Pow.hh"
35 #include "G4Threading.hh" 48 #include "G4Threading.hh"
36 49
37 G4Pow* G4Pow::fpInstance = nullptr; << 50 G4Pow* G4Pow::fpInstance = 0;
38 51
39 // ------------------------------------------- 52 // -------------------------------------------------------------------
40 53
41 G4Pow* G4Pow::GetInstance() 54 G4Pow* G4Pow::GetInstance()
42 { 55 {
43 if(fpInstance == nullptr) << 56 if (fpInstance == 0)
44 { 57 {
45 static G4Pow geant4pow; 58 static G4Pow geant4pow;
46 fpInstance = &geant4pow; 59 fpInstance = &geant4pow;
47 } 60 }
48 return fpInstance; 61 return fpInstance;
49 } 62 }
50 63
51 // ------------------------------------------- 64 // -------------------------------------------------------------------
52 65
53 G4Pow::G4Pow() 66 G4Pow::G4Pow()
>> 67 : onethird(1.0/3.0), max2(5)
54 { 68 {
55 #ifdef G4MULTITHREADED 69 #ifdef G4MULTITHREADED
56 if(G4Threading::IsWorkerThread()) 70 if(G4Threading::IsWorkerThread())
57 { 71 {
58 G4Exception("G4Pow::G4Pow()", "InvalidSetu << 72 G4Exception ("G4Pow::G4Pow()", "InvalidSetup", FatalException,
59 "Attempt to instantiate G4Pow << 73 "Attempt to instantiate G4Pow in worker thread!");
60 } 74 }
61 #endif 75 #endif
62 const G4int maxZ = 512; 76 const G4int maxZ = 512;
63 const G4int maxZfact = 170; 77 const G4int maxZfact = 170;
64 const G4int numLowA = 17; << 78 const G4int numLowA = 17;
65 79
66 maxA = -0.6 + maxZ; 80 maxA = -0.6 + maxZ;
67 maxLowA = 4.0; << 81 maxLowA = 4.0;
68 maxA2 = 1.25 + max2 * 0.2; << 82 maxA2 = 1.25 + max2*0.2;
69 maxAexp = -0.76 + maxZfact * 0.5; << 83 maxAexp = -0.76+ maxZfact*0.5;
70 << 84
71 ener.resize(max2 + 1, 1.0); << 85 ener.resize(max2+1,1.0);
72 logen.resize(max2 + 1, 0.0); << 86 logen.resize(max2+1,0.0);
73 lz2.resize(max2 + 1, 0.0); << 87 lz2.resize(max2+1,0.0);
74 pz13.resize(maxZ, 0.0); << 88 pz13.resize(maxZ,0.0);
75 lowa13.resize(numLowA, 0.0); << 89 lowa13.resize(numLowA,0.0);
76 lz.resize(maxZ, 0.0); << 90 lz.resize(maxZ,0.0);
77 fexp.resize(maxZfact, 0.0); << 91 fexp.resize(maxZfact,0.0);
78 fact.resize(maxZfact, 0.0); << 92 fact.resize(maxZfact,0.0);
79 logfact.resize(maxZ, 0.0); << 93 logfact.resize(maxZ,0.0);
80 94
81 G4double f = 1.0; << 95 G4double f = 1.0;
82 G4double logf = 0.0; 96 G4double logf = 0.0;
83 fact[0] = 1.0; << 97 fact[0] = 1.0;
84 fexp[0] = 1.0; << 98 fexp[0] = 1.0;
85 99
86 for(G4int i = 1; i <= max2; ++i) << 100 for(G4int i=1; i<=max2; ++i)
87 { 101 {
88 ener[i] = powN(500., i); << 102 ener[i] = powN(500.,i);
89 logen[i] = G4Log(ener[i]); << 103 logen[i]= G4Log(ener[i]);
90 lz2[i] = G4Log(1.0 + i * 0.2); << 104 lz2[i] = G4Log(1.0 + i*0.2);
91 } 105 }
92 106
93 for(G4int i = 1; i < maxZ; ++i) << 107 for(G4int i=1; i<maxZ; ++i)
94 { 108 {
95 auto x = G4double(i); << 109 G4double x = G4double(i);
96 pz13[i] = std::pow(x, onethird); << 110 pz13[i] = std::pow(x,onethird);
97 lz[i] = G4Log(x); << 111 lz[i] = G4Log(x);
98 if(i < maxZfact) 112 if(i < maxZfact)
99 { 113 {
100 f *= x; 114 f *= x;
101 fact[i] = f; 115 fact[i] = f;
102 fexp[i] = G4Exp(0.5 * i); << 116 fexp[i] = G4Exp(0.5*i);
103 } 117 }
104 logf += lz[i]; 118 logf += lz[i];
105 logfact[i] = logf; 119 logfact[i] = logf;
106 } 120 }
107 121
108 for(G4int i = 4; i < numLowA; ++i) << 122 for (G4int i=4; i<numLowA; ++i) {
109 { << 123 lowa13[i] = std::pow(0.25*i,onethird);
110 lowa13[i] = std::pow(0.25 * i, onethird); <<
111 } 124 }
>> 125
112 } 126 }
113 127
114 // ------------------------------------------- 128 // -------------------------------------------------------------------
115 129
116 G4double G4Pow::A13(G4double A) const << 130 G4Pow::~G4Pow()
117 { << 131 {}
>> 132
>> 133 // -------------------------------------------------------------------
>> 134
>> 135 G4double G4Pow::A13(G4double A) const {
118 G4double res = 0.; 136 G4double res = 0.;
119 if(A > 0.) << 137 if (A>0.) {
120 { << 138 const bool invert = (A<1.);
121 const bool invert = (A < 1.); << 139 const G4double a = invert ? 1./A : A;
122 const G4double a = invert ? 1. / A : A; << 140 res = (a<maxLowA) ? A13Low(a, invert) : A13High(a, invert);
123 res = (a < maxLowA) ? A13Low <<
124 } 141 }
125 return res; 142 return res;
126 } 143 }
127 144
128 // ------------------------------------------- 145 // -------------------------------------------------------------------
129 146
130 G4double G4Pow::A13High(const G4double a, cons << 147 G4double G4Pow::A13High(const G4double a, const bool invert) const {
131 { <<
132 G4double res; 148 G4double res;
133 if(a < maxA) << 149 if (a<maxA) {
134 { << 150 const G4int i = static_cast<G4int>(a+0.5);
135 const auto i = static_cast<G4int>(a + << 151 const G4double x = (a/i-1.)*onethird;
136 const G4double x = (a / i - 1.) * onethird << 152 res = pz13[i]*(1.+x-x*x*(1.-1.666667*x));
137 res = pz13[i] * (1. + x - x * << 153 } else {
138 } << 154 res = G4Exp(G4Log(a)*onethird);
139 else <<
140 { <<
141 res = G4Exp(G4Log(a) * onethird); <<
142 } 155 }
143 res = invert ? 1. / res : res; << 156 res = invert ? 1./res : res;
144 return res; 157 return res;
145 } 158 }
146 159
147 // ------------------------------------------- 160 // -------------------------------------------------------------------
148 161
149 G4double G4Pow::A13Low(const G4double a, const << 162 G4double G4Pow::A13Low(const G4double a, const bool invert) const {
150 { <<
151 G4double res; 163 G4double res;
152 const auto i = static_cast<G4int>(4. * ( << 164 const G4int i = static_cast<G4int>(4.*(a+0.125));
153 const G4double y = 0.25 * i; << 165 const G4double y = 0.25*i;
154 const G4double x = (a / y - 1.) * onethird; << 166 const G4double x = (a/y-1.)*onethird;
155 res = lowa13[i] * (1. + x - x * << 167 res = lowa13[i]*(1.+x-x*x*(1.-1.666667*x));
156 res = invert ? 1. / res : res; << 168 res = invert ? 1./res : res;
157 return res; 169 return res;
158 } 170 }
159 171
160 // ------------------------------------------- 172 // -------------------------------------------------------------------
161 173
162 G4double G4Pow::powN(G4double x, G4int n) cons 174 G4double G4Pow::powN(G4double x, G4int n) const
163 { 175 {
164 if(0.0 == x) << 176 if(0.0 == x) { return 0.0; }
165 { << 177 if(std::abs(n) > 8) { return std::pow(x, G4double(n)); }
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; 178 G4double res = 1.0;
173 if(n >= 0) << 179 if(n >= 0) { for(G4int i=0; i<n; ++i) { res *= x; } }
174 { <<
175 for(G4int i = 0; i < n; ++i) <<
176 { <<
177 res *= x; <<
178 } <<
179 } <<
180 else if(n < 0) 180 else if(n < 0)
181 { 181 {
182 G4double y = 1.0 / x; << 182 G4double y = 1.0/x;
183 G4int nn = -n; << 183 G4int nn = -n;
184 for(G4int i = 0; i < nn; ++i) << 184 for(G4int i=0; i<nn; ++i) { res *= y; }
185 { <<
186 res *= y; <<
187 } <<
188 } 185 }
189 return res; 186 return res;
190 } 187 }
191 188