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 // 27 // 28 // Author: V. Grichine (Vladimir,Grichine@cern.ch) 29 // 30 // 31 // G4 Model: diffuse optical elastic scattering with 4-momentum balance 32 // 33 // Class Description 34 // Final state production model for hadron nuclear elastic scattering; 35 // Class Description - End 36 // 37 // 38 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering 39 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering 40 // 12.06.11 V. Grichine, new interface to G4hadronElastic 41 // 24.11.17 W. Pokorski, code cleanup and performance improvements 42 43 44 #ifndef G4DiffuseElasticV2_h 45 #define G4DiffuseElasticV2_h 1 46 47 #include <CLHEP/Units/PhysicalConstants.h> 48 #include "globals.hh" 49 #include "G4HadronElastic.hh" 50 #include "G4HadProjectile.hh" 51 #include "G4Nucleus.hh" 52 53 #include "G4Pow.hh" 54 55 #include <vector> 56 57 class G4ParticleDefinition; 58 class G4PhysicsTable; 59 class G4PhysicsLogVector; 60 61 class G4DiffuseElasticV2 : public G4HadronElastic // G4HadronicInteraction 62 { 63 public: 64 65 G4DiffuseElasticV2(); 66 67 virtual ~G4DiffuseElasticV2(); 68 69 virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/, 70 G4Nucleus & /*targetNucleus*/); 71 72 void Initialise(); 73 74 void InitialiseOnFly(G4double Z, G4double A); 75 76 void BuildAngleTable(); 77 78 virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 79 G4double plab, 80 G4int Z, G4int A); 81 82 G4double NeutronTuniform(G4int Z); 83 84 void SetPlabLowLimit(G4double value); 85 86 void SetHEModelLowLimit(G4double value); 87 88 void SetQModelLowLimit(G4double value); 89 90 void SetLowestEnergyLimit(G4double value); 91 92 void SetRecoilKinEnergyLimit(G4double value); 93 94 G4double SampleTableT(const G4ParticleDefinition* aParticle, 95 G4double p, G4double Z, G4double A); 96 97 G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A); 98 99 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 100 G4double Z, G4double A); 101 102 G4double GetScatteringAngle(G4int iMomentum, unsigned long iAngle, G4double position); 103 104 G4double SampleThetaLab(const G4HadProjectile* aParticle, 105 G4double tmass, G4double A); 106 107 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 ); 108 109 G4double CalculateAm( G4double momentum, G4double n, G4double Z); 110 111 G4double CalculateNuclearRad( G4double A); 112 113 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 114 G4double tmass, G4double thetaCMS); 115 116 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 117 G4double tmass, G4double thetaLab); 118 119 120 G4double BesselJzero(G4double z); 121 G4double BesselJone(G4double z); 122 G4double DampFactor(G4double z); 123 G4double BesselOneByArg(G4double z); 124 125 G4double GetDiffElasticSumProbA(G4double alpha); 126 G4double GetIntegrandFunction(G4double theta); 127 128 129 G4double GetNuclearRadius(){return fNuclearRadius;}; 130 131 private: 132 133 134 G4ParticleDefinition* theProton; 135 G4ParticleDefinition* theNeutron; 136 137 G4double lowEnergyRecoilLimit; 138 G4double lowEnergyLimitHE; 139 G4double lowEnergyLimitQ; 140 G4double lowestEnergyLimit; 141 G4double plabLowLimit; 142 143 G4int fEnergyBin; 144 unsigned long fAngleBin; 145 146 G4PhysicsLogVector* fEnergyVector; 147 148 std::vector<std::vector<std::vector<double>*>*> fEnergyAngleVectorBank; 149 std::vector<std::vector<std::vector<double>*>*> fEnergySumVectorBank; 150 151 std::vector<std::vector<double>*>* fEnergyAngleVector; 152 std::vector<std::vector<double>*>* fEnergySumVector; 153 154 155 std::vector<G4double> fElementNumberVector; 156 std::vector<G4String> fElementNameVector; 157 158 const G4ParticleDefinition* fParticle; 159 G4double fWaveVector; 160 G4double fAtomicWeight; 161 G4double fAtomicNumber; 162 G4double fNuclearRadius; 163 G4double fBeta; 164 G4double fZommerfeld; 165 G4double fAm; 166 G4bool fAddCoulomb; 167 168 }; 169 170 inline G4bool G4DiffuseElasticV2::IsApplicable(const G4HadProjectile & projectile, 171 G4Nucleus & nucleus) 172 { 173 if( ( projectile.GetDefinition() == G4Proton::Proton() || 174 projectile.GetDefinition() == G4Neutron::Neutron() || 175 projectile.GetDefinition() == G4PionPlus::PionPlus() || 176 projectile.GetDefinition() == G4PionMinus::PionMinus() || 177 projectile.GetDefinition() == G4KaonPlus::KaonPlus() || 178 projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) && 179 180 nucleus.GetZ_asInt() >= 2 ) return true; 181 else return false; 182 } 183 184 inline void G4DiffuseElasticV2::SetRecoilKinEnergyLimit(G4double value) 185 { 186 lowEnergyRecoilLimit = value; 187 } 188 189 inline void G4DiffuseElasticV2::SetPlabLowLimit(G4double value) 190 { 191 plabLowLimit = value; 192 } 193 194 inline void G4DiffuseElasticV2::SetHEModelLowLimit(G4double value) 195 { 196 lowEnergyLimitHE = value; 197 } 198 199 inline void G4DiffuseElasticV2::SetQModelLowLimit(G4double value) 200 { 201 lowEnergyLimitQ = value; 202 } 203 204 inline void G4DiffuseElasticV2::SetLowestEnergyLimit(G4double value) 205 { 206 lowestEnergyLimit = value; 207 } 208 209 210 ///////////////////////////////////////////////////////////// 211 // 212 // Bessel J0 function based on rational approximation from 213 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 214 215 inline G4double G4DiffuseElasticV2::BesselJzero(G4double value) 216 { 217 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 218 219 modvalue = std::fabs(value); 220 221 if ( value < 8.0 && value > -8.0 ) 222 { 223 value2 = value*value; 224 225 fact1 = 57568490574.0 + value2*(-13362590354.0 226 + value2*( 651619640.7 227 + value2*(-11214424.18 228 + value2*( 77392.33017 229 + value2*(-184.9052456 ) ) ) ) ); 230 231 fact2 = 57568490411.0 + value2*( 1029532985.0 232 + value2*( 9494680.718 233 + value2*(59272.64853 234 + value2*(267.8532712 235 + value2*1.0 ) ) ) ); 236 237 bessel = fact1/fact2; 238 } 239 else 240 { 241 arg = 8.0/modvalue; 242 243 value2 = arg*arg; 244 245 shift = modvalue-0.785398164; 246 247 fact1 = 1.0 + value2*(-0.1098628627e-2 248 + value2*(0.2734510407e-4 249 + value2*(-0.2073370639e-5 250 + value2*0.2093887211e-6 ) ) ); 251 252 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3 253 + value2*(-0.6911147651e-5 254 + value2*(0.7621095161e-6 255 - value2*0.934945152e-7 ) ) ); 256 257 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 ); 258 } 259 return bessel; 260 } 261 262 ///////////////////////////////////////////////////////////// 263 // 264 // Bessel J1 function based on rational approximation from 265 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 266 267 inline G4double G4DiffuseElasticV2::BesselJone(G4double value) 268 { 269 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 270 271 modvalue = std::fabs(value); 272 273 if ( modvalue < 8.0 ) 274 { 275 value2 = value*value; 276 277 fact1 = value*(72362614232.0 + value2*(-7895059235.0 278 + value2*( 242396853.1 279 + value2*(-2972611.439 280 + value2*( 15704.48260 281 + value2*(-30.16036606 ) ) ) ) ) ); 282 283 fact2 = 144725228442.0 + value2*(2300535178.0 284 + value2*(18583304.74 285 + value2*(99447.43394 286 + value2*(376.9991397 287 + value2*1.0 ) ) ) ); 288 bessel = fact1/fact2; 289 } 290 else 291 { 292 arg = 8.0/modvalue; 293 294 value2 = arg*arg; 295 296 shift = modvalue - 2.356194491; 297 298 fact1 = 1.0 + value2*( 0.183105e-2 299 + value2*(-0.3516396496e-4 300 + value2*(0.2457520174e-5 301 + value2*(-0.240337019e-6 ) ) ) ); 302 303 fact2 = 0.04687499995 + value2*(-0.2002690873e-3 304 + value2*( 0.8449199096e-5 305 + value2*(-0.88228987e-6 306 + value2*0.105787412e-6 ) ) ); 307 308 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 309 310 if (value < 0.0) bessel = -bessel; 311 } 312 return bessel; 313 } 314 315 //////////////////////////////////////////////////////////////////// 316 // 317 // damp factor in diffraction x/sh(x), x was already *pi 318 319 inline G4double G4DiffuseElasticV2::DampFactor(G4double x) 320 { 321 G4double df; 322 G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials 323 324 // x *= pi; 325 326 if( std::fabs(x) < 0.01 ) 327 { 328 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4); 329 } 330 else 331 { 332 df = x/std::sinh(x); 333 } 334 return df; 335 } 336 337 338 //////////////////////////////////////////////////////////////////// 339 // 340 // return J1(x)/x with special case for small x 341 342 inline G4double G4DiffuseElasticV2::BesselOneByArg(G4double x) 343 { 344 G4double x2, result; 345 346 if( std::fabs(x) < 0.01 ) 347 { 348 x *= 0.5; 349 x2 = x*x; 350 result = 2. - x2 + x2*x2/6.; 351 } 352 else 353 { 354 result = BesselJone(x)/x; 355 } 356 return result; 357 } 358 359 360 //////////////////////////////////////////////////////////////////// 361 // 362 // return Zommerfeld parameter for Coulomb scattering 363 364 inline G4double G4DiffuseElasticV2::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 ) 365 { 366 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta; 367 368 return fZommerfeld; 369 } 370 371 //////////////////////////////////////////////////////////////////// 372 // 373 // return Wentzel correction for Coulomb scattering 374 375 inline G4double G4DiffuseElasticV2::CalculateAm( G4double momentum, G4double n, G4double Z) 376 { 377 G4double k = momentum/CLHEP::hbarc; 378 G4double ch = 1.13 + 3.76*n*n; 379 G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius; 380 G4double zn2 = zn*zn; 381 fAm = ch/zn2; 382 383 return fAm; 384 } 385 386 //////////////////////////////////////////////////////////////////// 387 // 388 // calculate nuclear radius for different atomic weights using different approximations 389 390 inline G4double G4DiffuseElasticV2::CalculateNuclearRad( G4double A) 391 { 392 G4double R, r0, a11, a12, a13, a2, a3; 393 394 a11 = 1.26; // 1.08, 1.16 395 a12 = 1.; // 1.08, 1.16 396 a13 = 1.12; // 1.08, 1.16 397 a2 = 1.1; 398 a3 = 1.; 399 400 // Special rms radii for light nucleii 401 402 if (A < 50.) 403 { 404 if (std::abs(A-1.) < 0.5) return 0.89*CLHEP::fermi; // p 405 else if(std::abs(A-2.) < 0.5) return 2.13*CLHEP::fermi; // d 406 else if( // std::abs(Z-1.) < 0.5 && 407 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t 408 409 // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3 410 else if( // std::abs(Z-2.) < 0.5 && 411 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4 412 413 else if( // std::abs(Z-3.) < 0.5 414 std::abs(A-7.) < 0.5 ) return 2.40*CLHEP::fermi; // Li7 415 else if( // std::abs(Z-4.) < 0.5 416 std::abs(A-9.) < 0.5) return 2.51*CLHEP::fermi; // Be9 417 418 else if( 10. < A && A <= 16. ) r0 = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi; 419 else if( 15. < A && A <= 20. ) r0 = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; 420 else if( 20. < A && A <= 30. ) r0 = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; 421 else r0 = a2*CLHEP::fermi; 422 423 R = r0*G4Pow::GetInstance()->A13(A); 424 } 425 else 426 { 427 r0 = a3*CLHEP::fermi; 428 429 R = r0*G4Pow::GetInstance()->powA(A, 0.27); 430 } 431 fNuclearRadius = R; 432 433 return R; 434 } 435 436 437 #endif 438