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