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 << 26 // $Id: G4Pow.hh 93311 2015-10-16 10:16:37Z gcosmo $ >> 27 // >> 28 // >> 29 // ------------------------------------------------------------------- >> 30 // >> 31 // Class G4Pow 27 // 32 // 28 // Class description: 33 // Class description: 29 // 34 // 30 // Utility singleton class for the fast comput 35 // Utility singleton class for the fast computation of log and pow 31 // functions. Integer argument should be in th << 36 // functions. Integer argument should in the interval 0-512, no 32 // check is performed inside these methods for 37 // check is performed inside these methods for performance reasons. 33 // For factorial integer argument should be in 38 // For factorial integer argument should be in the interval 0-170 34 // Computations with double arguments are fast 39 // Computations with double arguments are fast for the interval 35 // 0.002-511.5 for all functions except expone << 40 // 0.002-511.5 for all functions except exponent, which is computed 36 // for the interval 0-84.4, standard library i 41 // for the interval 0-84.4, standard library is used in the opposite case 37 42 38 // Author: Vladimir Ivanchenko, 23.05.2009 << 43 // Author: Vladimir Ivanchenko 39 // ------------------------------------------- << 44 // 40 #ifndef G4Pow_hh << 45 // Creation date: 23.05.2009 41 #define G4Pow_hh 1 << 46 // ------------------------------------------------------------------- >> 47 >> 48 #ifndef G4Pow_h >> 49 #define G4Pow_h 1 42 50 43 #include "G4DataVector.hh" << 44 #include "G4Exp.hh" << 45 #include "G4Log.hh" << 46 #include "globals.hh" 51 #include "globals.hh" >> 52 #include "G4Log.hh" >> 53 #include "G4Exp.hh" >> 54 #include "G4DataVector.hh" 47 55 48 class G4Pow 56 class G4Pow 49 { 57 { 50 public: << 58 51 static G4Pow* GetInstance(); << 59 public: 52 ~G4Pow() = default; << 60 53 << 61 static G4Pow* GetInstance(); 54 // Fast computation of Z^1/3 << 62 ~G4Pow(); 55 // << 63 56 inline G4double Z13(G4int Z) const; << 64 // Fast computation of Z^1/3 57 G4double A13(G4double A) const; << 65 // 58 << 66 inline G4double Z13(G4int Z) const; 59 // Fast computation of Z^2/3 << 67 inline G4double A13(G4double A) const; 60 // << 68 61 inline G4double Z23(G4int Z) const; << 69 // Fast computation of Z^2/3 62 inline G4double A23(G4double A) const; << 70 // 63 << 71 inline G4double Z23(G4int Z) const; 64 // Fast computation of log(Z) << 72 inline G4double A23(G4double A) const; 65 // << 73 66 inline G4double logZ(G4int Z) const; << 74 // Fast computation of log(Z) 67 inline G4double logA(G4double A) const; << 75 // 68 inline G4double logX(G4double x) const; << 76 inline G4double logZ(G4int Z) const; 69 << 77 inline G4double logA(G4double A) const; 70 // Fast computation of log10(Z) << 78 inline G4double logX(G4double x) const; 71 // << 79 72 inline G4double log10Z(G4int Z) const; << 80 // Fast computation of log10(Z) 73 inline G4double log10A(G4double A) const; << 81 // 74 << 82 inline G4double log10Z(G4int Z) const; 75 // Fast computation of exp(X) << 83 inline G4double log10A(G4double A) const; 76 // << 84 77 inline G4double expA(G4double A) const; << 85 // Fast computation of exp(X) 78 << 86 // 79 // Fast computation of pow(Z,X) << 87 inline G4double expA(G4double A) const; 80 // << 88 81 inline G4double powZ(G4int Z, G4double y) co << 89 // Fast computation of pow(Z,X) 82 inline G4double powA(G4double A, G4double y) << 90 // 83 G4double powN(G4double x, G4int n) const; << 91 inline G4double powZ(G4int Z, G4double y) const; 84 << 92 inline G4double powA(G4double A, G4double y) const; 85 // Fast factorial << 93 G4double powN(G4double x, G4int n) const; 86 // << 94 87 inline G4double factorial(G4int Z) const; << 95 // Fast factorial 88 inline G4double logfactorial(G4int Z) const; << 96 // 89 << 97 inline G4double factorial(G4int Z) const; 90 private: << 98 inline G4double logfactorial(G4int Z) const; 91 G4Pow(); << 99 92 << 100 private: 93 G4double A13Low(const G4double, const G4bool << 101 94 G4double A13High(const G4double, const G4boo << 102 G4Pow(); 95 << 103 96 inline G4double logBase(G4double x) const; << 104 inline G4double logBase(G4double x) const; 97 << 105 98 static G4Pow* fpInstance; << 106 static G4Pow* fpInstance; 99 << 107 100 const G4double onethird = 1.0 / 3.0; << 108 const G4double onethird; 101 const G4int max2 = 5; << 109 const G4int max2; 102 << 110 103 G4double maxA; << 111 G4double maxA; 104 G4double maxLowA; << 112 G4double maxA2; 105 G4double maxA2; << 113 G4double maxAexp; 106 G4double maxAexp; << 114 107 << 115 G4DataVector ener; 108 G4DataVector ener; << 116 G4DataVector logen; 109 G4DataVector logen; << 117 G4DataVector pz13; 110 G4DataVector pz13; << 118 G4DataVector lz; 111 G4DataVector lowa13; << 119 G4DataVector lz2; 112 G4DataVector lz; << 120 G4DataVector fexp; 113 G4DataVector lz2; << 121 G4DataVector fact; 114 G4DataVector fexp; << 122 G4DataVector logfact; 115 G4DataVector fact; << 116 G4DataVector logfact; << 117 }; 123 }; 118 124 119 // ----------------------------- << 125 // ------------------------------------------------------------------- 120 // Inline methods implementation << 126 121 // ----------------------------- << 127 inline G4double G4Pow::Z13(G4int Z) const >> 128 { >> 129 return pz13[Z]; >> 130 } 122 131 123 inline G4double G4Pow::Z13(G4int Z) const { re << 132 inline G4double G4Pow::A13(G4double A) const >> 133 { >> 134 G4double res = 0.0; >> 135 if(A > 0.0) >> 136 { >> 137 G4double a = (1.0 <= A) ? A : 1.0/A; >> 138 if(1.0 > A) { a = 1.0/A; } >> 139 if(a <= maxA) >> 140 { >> 141 G4int i = G4int(a + 0.5); >> 142 G4double x = (a/G4double(i) - 1.0)*onethird; >> 143 res = pz13[i]*(1.0 + x - x*x*(1.0 - 1.66666666*x)); >> 144 if(1.0 > A) { res = 1.0/res; } >> 145 } >> 146 else >> 147 { >> 148 res = std::pow(A, onethird); >> 149 } >> 150 } >> 151 return res; >> 152 } 124 153 125 inline G4double G4Pow::Z23(G4int Z) const 154 inline G4double G4Pow::Z23(G4int Z) const 126 { 155 { 127 G4double x = Z13(Z); 156 G4double x = Z13(Z); 128 return x * x; << 157 return x*x; 129 } 158 } 130 159 131 inline G4double G4Pow::A23(G4double A) const 160 inline G4double G4Pow::A23(G4double A) const 132 { 161 { 133 G4double x = A13(A); 162 G4double x = A13(A); 134 return x * x; << 163 return x*x; 135 } 164 } 136 165 137 inline G4double G4Pow::logZ(G4int Z) const { r << 166 inline G4double G4Pow::logZ(G4int Z) const >> 167 { >> 168 return lz[Z]; >> 169 } 138 170 139 inline G4double G4Pow::logBase(G4double a) con 171 inline G4double G4Pow::logBase(G4double a) const 140 { 172 { 141 G4double res; 173 G4double res; 142 if(a <= maxA2) << 174 if(a <= maxA2) 143 { 175 { 144 G4int i = G4int(max2 * (a - 1) + 0.5); << 176 G4int i = G4int(max2*(a - 1) + 0.5); 145 if(i > max2) << 177 if(i > max2) { i = max2; } 146 { << 178 G4double x = a/(G4double(i)/max2 + 1) - 1; 147 i = max2; << 179 res = lz2[i] + x*(1.0 - (0.5 - onethird*x)*x); 148 } << 149 G4double x = a / (G4double(i) / max2 + 1) << 150 res = lz2[i] + x * (1.0 - (0.5 - on << 151 } 180 } 152 else if(a <= maxA) 181 else if(a <= maxA) 153 { 182 { 154 G4int i = G4int(a + 0.5); << 183 G4int i = G4int(a + 0.5); 155 G4double x = a / G4double(i) - 1; << 184 G4double x = a/G4double(i) - 1; 156 res = lz[i] + x * (1.0 - (0.5 - one << 185 res = lz[i] + x*(1.0 - (0.5 - onethird*x)*x); 157 } 186 } 158 else 187 else 159 { 188 { 160 res = G4Log(a); 189 res = G4Log(a); 161 } 190 } 162 return res; 191 return res; 163 } 192 } 164 193 165 inline G4double G4Pow::logA(G4double A) const 194 inline G4double G4Pow::logA(G4double A) const 166 { 195 { 167 return (1.0 <= A ? logBase(A) : -logBase(1. << 196 return (1.0 <= A ? logBase(A) : -logBase(1./A)); 168 } 197 } 169 198 170 inline G4double G4Pow::logX(G4double x) const 199 inline G4double G4Pow::logX(G4double x) const 171 { 200 { 172 G4double res = 0.0; 201 G4double res = 0.0; 173 G4double a = (1.0 <= x) ? x : 1.0 / x; << 202 G4double a = (1.0 <= x) ? x : 1.0/x; 174 203 175 if(a <= maxA) << 204 if(a <= maxA) 176 { 205 { 177 res = logBase(a); 206 res = logBase(a); 178 } 207 } 179 else if(a <= ener[2]) 208 else if(a <= ener[2]) 180 { 209 { 181 res = logen[1] + logBase(a / ener[1]); << 210 res = logen[1] + logBase(a/ener[1]); 182 } 211 } 183 else if(a <= ener[3]) 212 else if(a <= ener[3]) 184 { 213 { 185 res = logen[2] + logBase(a / ener[2]); << 214 res = logen[2] + logBase(a/ener[2]); 186 } 215 } 187 else 216 else 188 { 217 { 189 res = G4Log(a); 218 res = G4Log(a); 190 } 219 } 191 220 192 if(1.0 > x) << 221 if(1.0 > x) { res = -res; } 193 { << 194 res = -res; << 195 } << 196 return res; 222 return res; 197 } 223 } 198 224 199 inline G4double G4Pow::log10Z(G4int Z) const { << 225 inline G4double G4Pow::log10Z(G4int Z) const >> 226 { >> 227 return lz[Z]/lz[10]; >> 228 } 200 229 201 inline G4double G4Pow::log10A(G4double A) cons << 230 inline G4double G4Pow::log10A(G4double A) const >> 231 { >> 232 return logX(A)/lz[10]; >> 233 } 202 234 203 inline G4double G4Pow::expA(G4double A) const 235 inline G4double G4Pow::expA(G4double A) const 204 { 236 { 205 G4double res; 237 G4double res; 206 G4double a = (0.0 <= A) ? A : -A; 238 G4double a = (0.0 <= A) ? A : -A; 207 239 208 if(a <= maxAexp) 240 if(a <= maxAexp) 209 { 241 { 210 G4int i = G4int(2 * a + 0.5); << 242 G4int i = G4int(2*a + 0.5); 211 G4double x = a - i * 0.5; << 243 G4double x = a - i*0.5; 212 res = fexp[i] * (1.0 + x * (1.0 + 0 << 244 res = fexp[i]*(1.0 + x*(1.0 + 0.5*(1.0 + onethird*x)*x)); 213 } 245 } 214 else 246 else 215 { 247 { 216 res = G4Exp(a); 248 res = G4Exp(a); 217 } 249 } 218 if(0.0 > A) << 250 if(0.0 > A) { res = 1.0/res; } 219 { << 220 res = 1.0 / res; << 221 } << 222 return res; 251 return res; 223 } 252 } 224 253 225 inline G4double G4Pow::powZ(G4int Z, G4double 254 inline G4double G4Pow::powZ(G4int Z, G4double y) const 226 { 255 { 227 return expA(y * lz[Z]); << 256 return expA(y*lz[Z]); 228 } 257 } 229 258 230 inline G4double G4Pow::powA(G4double A, G4doub 259 inline G4double G4Pow::powA(G4double A, G4double y) const 231 { 260 { 232 return (0.0 == A ? 0.0 : expA(y * logX(A))); << 261 return (0.0 == A ? 0.0 : expA(y*logX(A))); 233 } 262 } 234 263 235 inline G4double G4Pow::factorial(G4int Z) cons << 264 inline G4double G4Pow::factorial(G4int Z) const >> 265 { >> 266 return fact[Z]; >> 267 } >> 268 >> 269 inline G4double G4Pow::logfactorial(G4int Z) const >> 270 { >> 271 return logfact[Z]; >> 272 } 236 273 237 inline G4double G4Pow::logfactorial(G4int Z) c << 274 // ------------------------------------------------------------------- 238 275 239 #endif 276 #endif 240 277