00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef G4DiffuseElastic_h
00044 #define G4DiffuseElastic_h 1
00045
00046 #include <CLHEP/Units/PhysicalConstants.h>
00047 #include "globals.hh"
00048 #include "G4HadronElastic.hh"
00049 #include "G4HadProjectile.hh"
00050 #include "G4Nucleus.hh"
00051
00052 using namespace std;
00053
00054 class G4ParticleDefinition;
00055 class G4PhysicsTable;
00056 class G4PhysicsLogVector;
00057
00058 class G4DiffuseElastic : public G4HadronElastic
00059 {
00060 public:
00061
00062 G4DiffuseElastic();
00063
00064
00065
00066
00067
00068
00069
00070 virtual ~G4DiffuseElastic();
00071
00072 void Initialise();
00073
00074 void InitialiseOnFly(G4double Z, G4double A);
00075
00076 void BuildAngleTable();
00077
00078
00079
00080
00081 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
00082 G4double plab,
00083 G4int Z, G4int A);
00084
00085 void SetPlabLowLimit(G4double value);
00086
00087 void SetHEModelLowLimit(G4double value);
00088
00089 void SetQModelLowLimit(G4double value);
00090
00091 void SetLowestEnergyLimit(G4double value);
00092
00093 void SetRecoilKinEnergyLimit(G4double value);
00094
00095 G4double SampleT(const G4ParticleDefinition* aParticle,
00096 G4double p, G4double A);
00097
00098 G4double SampleTableT(const G4ParticleDefinition* aParticle,
00099 G4double p, G4double Z, G4double A);
00100
00101 G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
00102
00103 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
00104 G4double Z, G4double A);
00105
00106 G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
00107
00108 G4double SampleThetaLab(const G4HadProjectile* aParticle,
00109 G4double tmass, G4double A);
00110
00111 G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
00112 G4double theta,
00113 G4double momentum,
00114 G4double A );
00115
00116 G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
00117 G4double theta,
00118 G4double momentum,
00119 G4double A, G4double Z );
00120
00121 G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
00122 G4double theta,
00123 G4double momentum,
00124 G4double A, G4double Z );
00125
00126 G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
00127 G4double tMand,
00128 G4double momentum,
00129 G4double A, G4double Z );
00130
00131 G4double IntegralElasticProb( const G4ParticleDefinition* particle,
00132 G4double theta,
00133 G4double momentum,
00134 G4double A );
00135
00136
00137 G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
00138 G4double theta,
00139 G4double momentum,
00140 G4double Z );
00141
00142 G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
00143 G4double tMand,
00144 G4double momentum,
00145 G4double A, G4double Z );
00146
00147 G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
00148 G4double momentum, G4double Z );
00149
00150 G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
00151 G4double momentum, G4double Z,
00152 G4double theta1, G4double theta2 );
00153
00154
00155 G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
00156 G4double momentum );
00157
00158 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
00159
00160 G4double CalculateAm( G4double momentum, G4double n, G4double Z);
00161
00162 G4double CalculateNuclearRad( G4double A);
00163
00164 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
00165 G4double tmass, G4double thetaCMS);
00166
00167 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
00168 G4double tmass, G4double thetaLab);
00169
00170 void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
00171 G4double Z, G4double A);
00172
00173
00174
00175 G4double BesselJzero(G4double z);
00176 G4double BesselJone(G4double z);
00177 G4double DampFactor(G4double z);
00178 G4double BesselOneByArg(G4double z);
00179
00180 G4double GetDiffElasticProb(G4double theta);
00181 G4double GetDiffElasticSumProb(G4double theta);
00182 G4double GetDiffElasticSumProbA(G4double alpha);
00183 G4double GetIntegrandFunction(G4double theta);
00184
00185
00186 G4double GetNuclearRadius(){return fNuclearRadius;};
00187
00188 private:
00189
00190
00191 G4ParticleDefinition* theProton;
00192 G4ParticleDefinition* theNeutron;
00193 G4ParticleDefinition* theDeuteron;
00194 G4ParticleDefinition* theAlpha;
00195
00196 const G4ParticleDefinition* thePionPlus;
00197 const G4ParticleDefinition* thePionMinus;
00198
00199 G4double lowEnergyRecoilLimit;
00200 G4double lowEnergyLimitHE;
00201 G4double lowEnergyLimitQ;
00202 G4double lowestEnergyLimit;
00203 G4double plabLowLimit;
00204
00205 G4int fEnergyBin;
00206 G4int fAngleBin;
00207
00208 G4PhysicsLogVector* fEnergyVector;
00209 G4PhysicsTable* fAngleTable;
00210 std::vector<G4PhysicsTable*> fAngleBank;
00211
00212 std::vector<G4double> fElementNumberVector;
00213 std::vector<G4String> fElementNameVector;
00214
00215 const G4ParticleDefinition* fParticle;
00216 G4double fWaveVector;
00217 G4double fAtomicWeight;
00218 G4double fAtomicNumber;
00219 G4double fNuclearRadius;
00220 G4double fBeta;
00221 G4double fZommerfeld;
00222 G4double fAm;
00223 G4bool fAddCoulomb;
00224
00225 };
00226
00227
00228 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
00229 {
00230 lowEnergyRecoilLimit = value;
00231 }
00232
00233 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
00234 {
00235 plabLowLimit = value;
00236 }
00237
00238 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
00239 {
00240 lowEnergyLimitHE = value;
00241 }
00242
00243 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
00244 {
00245 lowEnergyLimitQ = value;
00246 }
00247
00248 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
00249 {
00250 lowestEnergyLimit = value;
00251 }
00252
00253
00255
00256
00257
00258
00259 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
00260 {
00261 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00262
00263 modvalue = fabs(value);
00264
00265 if ( value < 8.0 && value > -8.0 )
00266 {
00267 value2 = value*value;
00268
00269 fact1 = 57568490574.0 + value2*(-13362590354.0
00270 + value2*( 651619640.7
00271 + value2*(-11214424.18
00272 + value2*( 77392.33017
00273 + value2*(-184.9052456 ) ) ) ) );
00274
00275 fact2 = 57568490411.0 + value2*( 1029532985.0
00276 + value2*( 9494680.718
00277 + value2*(59272.64853
00278 + value2*(267.8532712
00279 + value2*1.0 ) ) ) );
00280
00281 bessel = fact1/fact2;
00282 }
00283 else
00284 {
00285 arg = 8.0/modvalue;
00286
00287 value2 = arg*arg;
00288
00289 shift = modvalue-0.785398164;
00290
00291 fact1 = 1.0 + value2*(-0.1098628627e-2
00292 + value2*(0.2734510407e-4
00293 + value2*(-0.2073370639e-5
00294 + value2*0.2093887211e-6 ) ) );
00295
00296 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
00297 + value2*(-0.6911147651e-5
00298 + value2*(0.7621095161e-6
00299 - value2*0.934945152e-7 ) ) );
00300
00301 bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
00302 }
00303 return bessel;
00304 }
00305
00307
00308
00309
00310
00311 inline G4double G4DiffuseElastic::BesselJone(G4double value)
00312 {
00313 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
00314
00315 modvalue = fabs(value);
00316
00317 if ( modvalue < 8.0 )
00318 {
00319 value2 = value*value;
00320
00321 fact1 = value*(72362614232.0 + value2*(-7895059235.0
00322 + value2*( 242396853.1
00323 + value2*(-2972611.439
00324 + value2*( 15704.48260
00325 + value2*(-30.16036606 ) ) ) ) ) );
00326
00327 fact2 = 144725228442.0 + value2*(2300535178.0
00328 + value2*(18583304.74
00329 + value2*(99447.43394
00330 + value2*(376.9991397
00331 + value2*1.0 ) ) ) );
00332 bessel = fact1/fact2;
00333 }
00334 else
00335 {
00336 arg = 8.0/modvalue;
00337
00338 value2 = arg*arg;
00339
00340 shift = modvalue - 2.356194491;
00341
00342 fact1 = 1.0 + value2*( 0.183105e-2
00343 + value2*(-0.3516396496e-4
00344 + value2*(0.2457520174e-5
00345 + value2*(-0.240337019e-6 ) ) ) );
00346
00347 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
00348 + value2*( 0.8449199096e-5
00349 + value2*(-0.88228987e-6
00350 + value2*0.105787412e-6 ) ) );
00351
00352 bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
00353
00354 if (value < 0.0) bessel = -bessel;
00355 }
00356 return bessel;
00357 }
00358
00360
00361
00362
00363 inline G4double G4DiffuseElastic::DampFactor(G4double x)
00364 {
00365 G4double df;
00366 G4double f2 = 2., f3 = 6., f4 = 24.;
00367
00368
00369
00370 if( std::fabs(x) < 0.01 )
00371 {
00372 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
00373 }
00374 else
00375 {
00376 df = x/std::sinh(x);
00377 }
00378 return df;
00379 }
00380
00381
00383
00384
00385
00386 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
00387 {
00388 G4double x2, result;
00389
00390 if( std::fabs(x) < 0.01 )
00391 {
00392 x *= 0.5;
00393 x2 = x*x;
00394 result = 2. - x2 + x2*x2/6.;
00395 }
00396 else
00397 {
00398 result = BesselJone(x)/x;
00399 }
00400 return result;
00401 }
00402
00404
00405
00406
00407 inline G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
00408 G4double momentum )
00409 {
00410 G4double mass = particle->GetPDGMass();
00411 G4double a = momentum/mass;
00412 fBeta = a/std::sqrt(1+a*a);
00413
00414 return fBeta;
00415 }
00416
00418
00419
00420
00421 inline G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
00422 {
00423 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
00424
00425 return fZommerfeld;
00426 }
00427
00429
00430
00431
00432 inline G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
00433 {
00434 G4double k = momentum/CLHEP::hbarc;
00435 G4double ch = 1.13 + 3.76*n*n;
00436 G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
00437 G4double zn2 = zn*zn;
00438 fAm = ch/zn2;
00439
00440 return fAm;
00441 }
00442
00444
00445
00446
00447 inline G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
00448 {
00449 G4double r0;
00450
00451 if( A < 50. )
00452 {
00453 if( A > 10. ) r0 = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi;
00454 else r0 = 1.1*CLHEP::fermi;
00455
00456 fNuclearRadius = r0*std::pow(A, 1./3.);
00457 }
00458 else
00459 {
00460 r0 = 1.7*CLHEP::fermi;
00461
00462 fNuclearRadius = r0*std::pow(A, 0.27);
00463 }
00464 return fNuclearRadius;
00465 }
00466
00468
00469
00470
00471 inline G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
00472 G4double theta,
00473 G4double momentum,
00474 G4double Z )
00475 {
00476 G4double sinHalfTheta = std::sin(0.5*theta);
00477 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00478 G4double beta = CalculateParticleBeta( particle, momentum);
00479 G4double z = particle->GetPDGCharge();
00480 G4double n = CalculateZommerfeld( beta, z, Z );
00481 G4double am = CalculateAm( momentum, n, Z);
00482 G4double k = momentum/CLHEP::hbarc;
00483 G4double ch = 0.5*n/k;
00484 G4double ch2 = ch*ch;
00485 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
00486
00487 return xsc;
00488 }
00489
00490
00492
00493
00494
00495 inline G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
00496 G4double momentum, G4double Z )
00497 {
00498 G4double beta = CalculateParticleBeta( particle, momentum);
00499 G4cout<<"beta = "<<beta<<G4endl;
00500 G4double z = particle->GetPDGCharge();
00501 G4double n = CalculateZommerfeld( beta, z, Z );
00502 G4cout<<"fZomerfeld = "<<n<<G4endl;
00503 G4double am = CalculateAm( momentum, n, Z);
00504 G4cout<<"cof Am = "<<am<<G4endl;
00505 G4double k = momentum/CLHEP::hbarc;
00506 G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
00507 G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
00508 G4double ch = n/k;
00509 G4double ch2 = ch*ch;
00510 G4double xsc = ch2*CLHEP::pi/(am +am*am);
00511
00512 return xsc;
00513 }
00514
00516
00517
00518
00519
00520 inline G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
00521 G4double momentum, G4double Z,
00522 G4double theta1, G4double theta2 )
00523 {
00524 G4double c1 = std::cos(theta1);
00525 G4cout<<"c1 = "<<c1<<G4endl;
00526 G4double c2 = std::cos(theta2);
00527 G4cout<<"c2 = "<<c2<<G4endl;
00528 G4double beta = CalculateParticleBeta( particle, momentum);
00529
00530 G4double z = particle->GetPDGCharge();
00531 G4double n = CalculateZommerfeld( beta, z, Z );
00532
00533 G4double am = CalculateAm( momentum, n, Z);
00534
00535 G4double k = momentum/CLHEP::hbarc;
00536
00537
00538 G4double ch = n/k;
00539 G4double ch2 = ch*ch;
00540 am *= 2.;
00541 G4double xsc = ch2*CLHEP::twopi*(c1-c2);
00542 xsc /= (1 - c1 + am)*(1 - c2 + am);
00543
00544 return xsc;
00545 }
00546
00547 #endif