File indexing completed on 2025-01-18 09:58:06
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042 #ifndef G4DiffuseElastic_h
0043 #define G4DiffuseElastic_h 1
0044
0045 #include <CLHEP/Units/PhysicalConstants.h>
0046 #include "globals.hh"
0047 #include "G4HadronElastic.hh"
0048 #include "G4HadProjectile.hh"
0049 #include "G4Nucleus.hh"
0050
0051 #include "G4Pow.hh"
0052
0053
0054 class G4ParticleDefinition;
0055 class G4PhysicsTable;
0056 class G4PhysicsLogVector;
0057
0058 class G4DiffuseElastic : public G4HadronElastic
0059 {
0060 public:
0061
0062 G4DiffuseElastic();
0063
0064
0065
0066
0067
0068
0069
0070 virtual ~G4DiffuseElastic();
0071
0072 virtual G4bool IsApplicable(const G4HadProjectile &,
0073 G4Nucleus & );
0074
0075 void Initialise();
0076
0077 void InitialiseOnFly(G4double Z, G4double A);
0078
0079 void BuildAngleTable();
0080
0081
0082
0083
0084 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
0085 G4double plab,
0086 G4int Z, G4int A);
0087
0088 G4double NeutronTuniform(G4int Z);
0089
0090 void SetPlabLowLimit(G4double value);
0091
0092 void SetHEModelLowLimit(G4double value);
0093
0094 void SetQModelLowLimit(G4double value);
0095
0096 void SetLowestEnergyLimit(G4double value);
0097
0098 void SetRecoilKinEnergyLimit(G4double value);
0099
0100 G4double SampleT(const G4ParticleDefinition* aParticle,
0101 G4double p, G4double A);
0102
0103 G4double SampleTableT(const G4ParticleDefinition* aParticle,
0104 G4double p, G4double Z, G4double A);
0105
0106 G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
0107
0108 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
0109 G4double Z, G4double A);
0110
0111 G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
0112
0113 G4double SampleThetaLab(const G4HadProjectile* aParticle,
0114 G4double tmass, G4double A);
0115
0116 G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
0117 G4double theta,
0118 G4double momentum,
0119 G4double A );
0120
0121 G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
0122 G4double theta,
0123 G4double momentum,
0124 G4double A, G4double Z );
0125
0126 G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
0127 G4double theta,
0128 G4double momentum,
0129 G4double A, G4double Z );
0130
0131 G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
0132 G4double tMand,
0133 G4double momentum,
0134 G4double A, G4double Z );
0135
0136 G4double IntegralElasticProb( const G4ParticleDefinition* particle,
0137 G4double theta,
0138 G4double momentum,
0139 G4double A );
0140
0141
0142 G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
0143 G4double theta,
0144 G4double momentum,
0145 G4double Z );
0146
0147 G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
0148 G4double tMand,
0149 G4double momentum,
0150 G4double A, G4double Z );
0151
0152 G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
0153 G4double momentum, G4double Z );
0154
0155 G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
0156 G4double momentum, G4double Z,
0157 G4double theta1, G4double theta2 );
0158
0159
0160 G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
0161 G4double momentum );
0162
0163 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
0164
0165 G4double CalculateAm( G4double momentum, G4double n, G4double Z);
0166
0167 G4double CalculateNuclearRad( G4double A);
0168
0169 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
0170 G4double tmass, G4double thetaCMS);
0171
0172 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
0173 G4double tmass, G4double thetaLab);
0174
0175 void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
0176 G4double Z, G4double A);
0177
0178
0179
0180 G4double BesselJzero(G4double z);
0181 G4double BesselJone(G4double z);
0182 G4double DampFactor(G4double z);
0183 G4double BesselOneByArg(G4double z);
0184
0185 G4double GetDiffElasticProb(G4double theta);
0186 G4double GetDiffElasticSumProb(G4double theta);
0187 G4double GetDiffElasticSumProbA(G4double alpha);
0188 G4double GetIntegrandFunction(G4double theta);
0189
0190
0191 G4double GetNuclearRadius(){return fNuclearRadius;};
0192
0193 private:
0194
0195
0196 G4ParticleDefinition* theProton;
0197 G4ParticleDefinition* theNeutron;
0198 G4ParticleDefinition* theDeuteron;
0199 G4ParticleDefinition* theAlpha;
0200
0201 const G4ParticleDefinition* thePionPlus;
0202 const G4ParticleDefinition* thePionMinus;
0203
0204 G4double lowEnergyRecoilLimit;
0205 G4double lowEnergyLimitHE;
0206 G4double lowEnergyLimitQ;
0207 G4double lowestEnergyLimit;
0208 G4double plabLowLimit;
0209
0210 G4int fEnergyBin;
0211 G4int fAngleBin;
0212
0213 G4PhysicsLogVector* fEnergyVector;
0214 G4PhysicsTable* fAngleTable;
0215 std::vector<G4PhysicsTable*> fAngleBank;
0216
0217 std::vector<G4double> fElementNumberVector;
0218 std::vector<G4String> fElementNameVector;
0219
0220 const G4ParticleDefinition* fParticle;
0221 G4double fWaveVector;
0222 G4double fAtomicWeight;
0223 G4double fAtomicNumber;
0224 G4double fNuclearRadius;
0225 G4double fBeta;
0226 G4double fZommerfeld;
0227 G4double fAm;
0228 G4bool fAddCoulomb;
0229
0230 };
0231
0232 inline G4bool G4DiffuseElastic::IsApplicable(const G4HadProjectile & projectile,
0233 G4Nucleus & nucleus)
0234 {
0235 if( ( projectile.GetDefinition() == G4Proton::Proton() ||
0236 projectile.GetDefinition() == G4Neutron::Neutron() ||
0237 projectile.GetDefinition() == G4PionPlus::PionPlus() ||
0238 projectile.GetDefinition() == G4PionMinus::PionMinus() ||
0239 projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
0240 projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
0241
0242 nucleus.GetZ_asInt() >= 2 ) return true;
0243 else return false;
0244 }
0245
0246 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
0247 {
0248 lowEnergyRecoilLimit = value;
0249 }
0250
0251 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
0252 {
0253 plabLowLimit = value;
0254 }
0255
0256 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
0257 {
0258 lowEnergyLimitHE = value;
0259 }
0260
0261 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
0262 {
0263 lowEnergyLimitQ = value;
0264 }
0265
0266 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
0267 {
0268 lowestEnergyLimit = value;
0269 }
0270
0271
0272
0273
0274
0275
0276
0277 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
0278 {
0279 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0280
0281 modvalue = std::fabs(value);
0282
0283 if ( value < 8.0 && value > -8.0 )
0284 {
0285 value2 = value*value;
0286
0287 fact1 = 57568490574.0 + value2*(-13362590354.0
0288 + value2*( 651619640.7
0289 + value2*(-11214424.18
0290 + value2*( 77392.33017
0291 + value2*(-184.9052456 ) ) ) ) );
0292
0293 fact2 = 57568490411.0 + value2*( 1029532985.0
0294 + value2*( 9494680.718
0295 + value2*(59272.64853
0296 + value2*(267.8532712
0297 + value2*1.0 ) ) ) );
0298
0299 bessel = fact1/fact2;
0300 }
0301 else
0302 {
0303 arg = 8.0/modvalue;
0304
0305 value2 = arg*arg;
0306
0307 shift = modvalue-0.785398164;
0308
0309 fact1 = 1.0 + value2*(-0.1098628627e-2
0310 + value2*(0.2734510407e-4
0311 + value2*(-0.2073370639e-5
0312 + value2*0.2093887211e-6 ) ) );
0313
0314 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
0315 + value2*(-0.6911147651e-5
0316 + value2*(0.7621095161e-6
0317 - value2*0.934945152e-7 ) ) );
0318
0319 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
0320 }
0321 return bessel;
0322 }
0323
0324
0325
0326
0327
0328
0329 inline G4double G4DiffuseElastic::BesselJone(G4double value)
0330 {
0331 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0332
0333 modvalue = std::fabs(value);
0334
0335 if ( modvalue < 8.0 )
0336 {
0337 value2 = value*value;
0338
0339 fact1 = value*(72362614232.0 + value2*(-7895059235.0
0340 + value2*( 242396853.1
0341 + value2*(-2972611.439
0342 + value2*( 15704.48260
0343 + value2*(-30.16036606 ) ) ) ) ) );
0344
0345 fact2 = 144725228442.0 + value2*(2300535178.0
0346 + value2*(18583304.74
0347 + value2*(99447.43394
0348 + value2*(376.9991397
0349 + value2*1.0 ) ) ) );
0350 bessel = fact1/fact2;
0351 }
0352 else
0353 {
0354 arg = 8.0/modvalue;
0355
0356 value2 = arg*arg;
0357
0358 shift = modvalue - 2.356194491;
0359
0360 fact1 = 1.0 + value2*( 0.183105e-2
0361 + value2*(-0.3516396496e-4
0362 + value2*(0.2457520174e-5
0363 + value2*(-0.240337019e-6 ) ) ) );
0364
0365 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
0366 + value2*( 0.8449199096e-5
0367 + value2*(-0.88228987e-6
0368 + value2*0.105787412e-6 ) ) );
0369
0370 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
0371
0372 if (value < 0.0) bessel = -bessel;
0373 }
0374 return bessel;
0375 }
0376
0377
0378
0379
0380
0381 inline G4double G4DiffuseElastic::DampFactor(G4double x)
0382 {
0383 G4double df;
0384 G4double f2 = 2., f3 = 6., f4 = 24.;
0385
0386
0387
0388 if( std::fabs(x) < 0.01 )
0389 {
0390 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
0391 }
0392 else
0393 {
0394 df = x/std::sinh(x);
0395 }
0396 return df;
0397 }
0398
0399
0400
0401
0402
0403
0404 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
0405 {
0406 G4double x2, result;
0407
0408 if( std::fabs(x) < 0.01 )
0409 {
0410 x *= 0.5;
0411 x2 = x*x;
0412 result = 2. - x2 + x2*x2/6.;
0413 }
0414 else
0415 {
0416 result = BesselJone(x)/x;
0417 }
0418 return result;
0419 }
0420
0421
0422
0423
0424
0425 inline G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
0426 G4double momentum )
0427 {
0428 G4double mass = particle->GetPDGMass();
0429 G4double a = momentum/mass;
0430 fBeta = a/std::sqrt(1+a*a);
0431
0432 return fBeta;
0433 }
0434
0435
0436
0437
0438
0439 inline G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
0440 {
0441 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
0442
0443 return fZommerfeld;
0444 }
0445
0446
0447
0448
0449
0450 inline G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
0451 {
0452 G4double k = momentum/CLHEP::hbarc;
0453 G4double ch = 1.13 + 3.76*n*n;
0454 G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
0455 G4double zn2 = zn*zn;
0456 fAm = ch/zn2;
0457
0458 return fAm;
0459 }
0460
0461
0462
0463
0464
0465 inline G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
0466 {
0467 G4double R, r0, a11, a12, a13, a2, a3;
0468
0469 a11 = 1.26;
0470 a12 = 1.;
0471 a13 = 1.12;
0472 a2 = 1.1;
0473 a3 = 1.;
0474
0475
0476
0477 if (A < 50.)
0478 {
0479 if (std::abs(A-1.) < 0.5) return 0.89*CLHEP::fermi;
0480 else if(std::abs(A-2.) < 0.5) return 2.13*CLHEP::fermi;
0481 else if(
0482 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi;
0483
0484
0485 else if(
0486 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi;
0487
0488 else if(
0489 std::abs(A-7.) < 0.5 ) return 2.40*CLHEP::fermi;
0490 else if(
0491 std::abs(A-9.) < 0.5) return 2.51*CLHEP::fermi;
0492
0493 else if( 10. < A && A <= 16. ) r0 = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0494 else if( 15. < A && A <= 20. ) r0 = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0495 else if( 20. < A && A <= 30. ) r0 = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0496 else r0 = a2*CLHEP::fermi;
0497
0498 R = r0*G4Pow::GetInstance()->A13(A);
0499 }
0500 else
0501 {
0502 r0 = a3*CLHEP::fermi;
0503
0504 R = r0*G4Pow::GetInstance()->powA(A, 0.27);
0505 }
0506 fNuclearRadius = R;
0507 return R;
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522
0523 }
0524
0525
0526
0527
0528
0529 inline G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
0530 G4double theta,
0531 G4double momentum,
0532 G4double Z )
0533 {
0534 G4double sinHalfTheta = std::sin(0.5*theta);
0535 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
0536 G4double beta = CalculateParticleBeta( particle, momentum);
0537 G4double z = particle->GetPDGCharge();
0538 G4double n = CalculateZommerfeld( beta, z, Z );
0539 G4double am = CalculateAm( momentum, n, Z);
0540 G4double k = momentum/CLHEP::hbarc;
0541 G4double ch = 0.5*n/k;
0542 G4double ch2 = ch*ch;
0543 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
0544
0545 return xsc;
0546 }
0547
0548
0549
0550
0551
0552
0553 inline G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
0554 G4double momentum, G4double Z )
0555 {
0556 G4double beta = CalculateParticleBeta( particle, momentum);
0557 G4cout<<"beta = "<<beta<<G4endl;
0558 G4double z = particle->GetPDGCharge();
0559 G4double n = CalculateZommerfeld( beta, z, Z );
0560 G4cout<<"fZomerfeld = "<<n<<G4endl;
0561 G4double am = CalculateAm( momentum, n, Z);
0562 G4cout<<"cof Am = "<<am<<G4endl;
0563 G4double k = momentum/CLHEP::hbarc;
0564 G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
0565 G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
0566 G4double ch = n/k;
0567 G4double ch2 = ch*ch;
0568 G4double xsc = ch2*CLHEP::pi/(am +am*am);
0569
0570 return xsc;
0571 }
0572
0573
0574
0575
0576
0577
0578 inline G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
0579 G4double momentum, G4double Z,
0580 G4double theta1, G4double theta2 )
0581 {
0582 G4double c1 = std::cos(theta1);
0583 G4cout<<"c1 = "<<c1<<G4endl;
0584 G4double c2 = std::cos(theta2);
0585 G4cout<<"c2 = "<<c2<<G4endl;
0586 G4double beta = CalculateParticleBeta( particle, momentum);
0587
0588 G4double z = particle->GetPDGCharge();
0589 G4double n = CalculateZommerfeld( beta, z, Z );
0590
0591 G4double am = CalculateAm( momentum, n, Z);
0592
0593 G4double k = momentum/CLHEP::hbarc;
0594
0595
0596 G4double ch = n/k;
0597 G4double ch2 = ch*ch;
0598 am *= 2.;
0599 G4double xsc = ch2*CLHEP::twopi*(c1-c2);
0600 xsc /= (1 - c1 + am)*(1 - c2 + am);
0601
0602 return xsc;
0603 }
0604
0605 #endif