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 #include "G4DiffuseElastic.hh"
00038 #include "G4ParticleTable.hh"
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4IonTable.hh"
00041 #include "G4NucleiProperties.hh"
00042
00043 #include "Randomize.hh"
00044 #include "G4Integrator.hh"
00045 #include "globals.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048
00049 #include "G4Proton.hh"
00050 #include "G4Neutron.hh"
00051 #include "G4Deuteron.hh"
00052 #include "G4Alpha.hh"
00053 #include "G4PionPlus.hh"
00054 #include "G4PionMinus.hh"
00055
00056 #include "G4Element.hh"
00057 #include "G4ElementTable.hh"
00058 #include "G4PhysicsTable.hh"
00059 #include "G4PhysicsLogVector.hh"
00060 #include "G4PhysicsFreeVector.hh"
00061
00063
00064
00065
00066
00067 G4DiffuseElastic::G4DiffuseElastic()
00068 : G4HadronElastic("DiffuseElastic"), fParticle(0)
00069 {
00070 SetMinEnergy( 0.01*GeV );
00071 SetMaxEnergy( 1.*TeV );
00072 verboseLevel = 0;
00073 lowEnergyRecoilLimit = 100.*keV;
00074 lowEnergyLimitQ = 0.0*GeV;
00075 lowEnergyLimitHE = 0.0*GeV;
00076 lowestEnergyLimit= 0.0*keV;
00077 plabLowLimit = 20.0*MeV;
00078
00079 theProton = G4Proton::Proton();
00080 theNeutron = G4Neutron::Neutron();
00081 theDeuteron = G4Deuteron::Deuteron();
00082 theAlpha = G4Alpha::Alpha();
00083 thePionPlus = G4PionPlus::PionPlus();
00084 thePionMinus= G4PionMinus::PionMinus();
00085
00086 fEnergyBin = 200;
00087 fAngleBin = 200;
00088
00089 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
00090 fAngleTable = 0;
00091
00092 fParticle = 0;
00093 fWaveVector = 0.;
00094 fAtomicWeight = 0.;
00095 fAtomicNumber = 0.;
00096 fNuclearRadius = 0.;
00097 fBeta = 0.;
00098 fZommerfeld = 0.;
00099 fAm = 0.;
00100 fAddCoulomb = false;
00101 }
00102
00104
00105
00106
00107 G4DiffuseElastic::~G4DiffuseElastic()
00108 {
00109 if(fEnergyVector) delete fEnergyVector;
00110
00111 if( fAngleTable )
00112 {
00113 fAngleTable->clearAndDestroy();
00114 delete fAngleTable ;
00115 }
00116 }
00117
00119
00120
00121
00122 void G4DiffuseElastic::Initialise()
00123 {
00124
00125
00126
00127 const G4ElementTable* theElementTable = G4Element::GetElementTable();
00128
00129 size_t jEl, numOfEl = G4Element::GetNumberOfElements();
00130
00131 for(jEl = 0 ; jEl < numOfEl; ++jEl)
00132 {
00133 fAtomicNumber = (*theElementTable)[jEl]->GetZ();
00134 fAtomicWeight = (*theElementTable)[jEl]->GetN();
00135 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
00136
00137 if(verboseLevel > 0)
00138 {
00139 G4cout<<"G4DiffuseElastic::Initialise() the element: "
00140 <<(*theElementTable)[jEl]->GetName()<<G4endl;
00141 }
00142 fElementNumberVector.push_back(fAtomicNumber);
00143 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
00144
00145 BuildAngleTable();
00146 fAngleBank.push_back(fAngleTable);
00147 }
00148 return;
00149 }
00150
00152
00153
00154
00155 G4double
00156 G4DiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
00157 G4double theta,
00158 G4double momentum,
00159 G4double A )
00160 {
00161 fParticle = particle;
00162 fWaveVector = momentum/hbarc;
00163 fAtomicWeight = A;
00164 fAddCoulomb = false;
00165 fNuclearRadius = CalculateNuclearRad(A);
00166
00167 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
00168
00169 return sigma;
00170 }
00171
00173
00174
00175
00176 G4double
00177 G4DiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle,
00178 G4double tMand,
00179 G4double plab,
00180 G4double A, G4double Z )
00181 {
00182 G4double m1 = particle->GetPDGMass();
00183 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00184
00185 G4int iZ = static_cast<G4int>(Z+0.5);
00186 G4int iA = static_cast<G4int>(A+0.5);
00187 G4ParticleDefinition * theDef = 0;
00188
00189 if (iZ == 1 && iA == 1) theDef = theProton;
00190 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00191 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00192 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00193 else if (iZ == 2 && iA == 4) theDef = theAlpha;
00194 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00195
00196 G4double tmass = theDef->GetPDGMass();
00197
00198 G4LorentzVector lv(0.0,0.0,0.0,tmass);
00199 lv += lv1;
00200
00201 G4ThreeVector bst = lv.boostVector();
00202 lv1.boost(-bst);
00203
00204 G4ThreeVector p1 = lv1.vect();
00205 G4double ptot = p1.mag();
00206 G4double ptot2 = ptot*ptot;
00207 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00208
00209 if( cost >= 1.0 ) cost = 1.0;
00210 else if( cost <= -1.0) cost = -1.0;
00211
00212 G4double thetaCMS = std::acos(cost);
00213
00214 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
00215
00216 sigma *= pi/ptot2;
00217
00218 return sigma;
00219 }
00220
00222
00223
00224
00225
00226 G4double
00227 G4DiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
00228 G4double theta,
00229 G4double momentum,
00230 G4double A, G4double Z )
00231 {
00232 fParticle = particle;
00233 fWaveVector = momentum/hbarc;
00234 fAtomicWeight = A;
00235 fAtomicNumber = Z;
00236 fNuclearRadius = CalculateNuclearRad(A);
00237 fAddCoulomb = false;
00238
00239 G4double z = particle->GetPDGCharge();
00240
00241 G4double kRt = fWaveVector*fNuclearRadius*theta;
00242 G4double kRtC = 1.9;
00243
00244 if( z && (kRt > kRtC) )
00245 {
00246 fAddCoulomb = true;
00247 fBeta = CalculateParticleBeta( particle, momentum);
00248 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
00249 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
00250 }
00251 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
00252
00253 return sigma;
00254 }
00255
00257
00258
00259
00260
00261 G4double
00262 G4DiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle,
00263 G4double tMand,
00264 G4double plab,
00265 G4double A, G4double Z )
00266 {
00267 G4double m1 = particle->GetPDGMass();
00268
00269 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00270
00271 G4int iZ = static_cast<G4int>(Z+0.5);
00272 G4int iA = static_cast<G4int>(A+0.5);
00273
00274 G4ParticleDefinition* theDef = 0;
00275
00276 if (iZ == 1 && iA == 1) theDef = theProton;
00277 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00278 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00279 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00280 else if (iZ == 2 && iA == 4) theDef = theAlpha;
00281 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00282
00283 G4double tmass = theDef->GetPDGMass();
00284
00285 G4LorentzVector lv(0.0,0.0,0.0,tmass);
00286 lv += lv1;
00287
00288 G4ThreeVector bst = lv.boostVector();
00289 lv1.boost(-bst);
00290
00291 G4ThreeVector p1 = lv1.vect();
00292 G4double ptot = p1.mag();
00293 G4double ptot2 = ptot*ptot;
00294 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00295
00296 if( cost >= 1.0 ) cost = 1.0;
00297 else if( cost <= -1.0) cost = -1.0;
00298
00299 G4double thetaCMS = std::acos(cost);
00300
00301 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
00302
00303 sigma *= pi/ptot2;
00304
00305 return sigma;
00306 }
00307
00309
00310
00311
00312
00313 G4double
00314 G4DiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
00315 G4double tMand,
00316 G4double plab,
00317 G4double A, G4double Z )
00318 {
00319 G4double m1 = particle->GetPDGMass();
00320 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
00321
00322 G4int iZ = static_cast<G4int>(Z+0.5);
00323 G4int iA = static_cast<G4int>(A+0.5);
00324 G4ParticleDefinition * theDef = 0;
00325
00326 if (iZ == 1 && iA == 1) theDef = theProton;
00327 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
00328 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
00329 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
00330 else if (iZ == 2 && iA == 4) theDef = theAlpha;
00331 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
00332
00333 G4double tmass = theDef->GetPDGMass();
00334
00335 G4LorentzVector lv(0.0,0.0,0.0,tmass);
00336 lv += lv1;
00337
00338 G4ThreeVector bst = lv.boostVector();
00339 lv1.boost(-bst);
00340
00341 G4ThreeVector p1 = lv1.vect();
00342 G4double ptot = p1.mag();
00343 G4double ptot2 = ptot*ptot;
00344 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
00345
00346 if( cost >= 1.0 ) cost = 1.0;
00347 else if( cost <= -1.0) cost = -1.0;
00348
00349 G4double thetaCMS = std::acos(cost);
00350
00351 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
00352
00353 sigma *= pi/ptot2;
00354
00355 return sigma;
00356 }
00357
00359
00360
00361
00362 G4double
00363 G4DiffuseElastic::GetDiffElasticProb(
00364 G4double theta
00365
00366
00367 )
00368 {
00369 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00370 G4double delta, diffuse, gamma;
00371 G4double e1, e2, bone, bone2;
00372
00373
00374
00375
00376
00377 G4double kr = fWaveVector*fNuclearRadius;
00378 G4double kr2 = kr*kr;
00379 G4double krt = kr*theta;
00380
00381 bzero = BesselJzero(krt);
00382 bzero2 = bzero*bzero;
00383 bone = BesselJone(krt);
00384 bone2 = bone*bone;
00385 bonebyarg = BesselOneByArg(krt);
00386 bonebyarg2 = bonebyarg*bonebyarg;
00387
00388 if (fParticle == theProton)
00389 {
00390 diffuse = 0.63*fermi;
00391 gamma = 0.3*fermi;
00392 delta = 0.1*fermi*fermi;
00393 e1 = 0.3*fermi;
00394 e2 = 0.35*fermi;
00395 }
00396 else
00397 {
00398 diffuse = 0.63*fermi;
00399 gamma = 0.3*fermi;
00400 delta = 0.1*fermi*fermi;
00401 e1 = 0.3*fermi;
00402 e2 = 0.35*fermi;
00403 }
00404 G4double lambda = 15.;
00405
00406
00407
00408 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
00409 G4double kgamma2 = kgamma*kgamma;
00410
00411
00412
00413
00414
00415 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));
00416
00417 damp = DampFactor(pikdt);
00418 damp2 = damp*damp;
00419
00420 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
00421 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00422
00423
00424 sigma = kgamma2;
00425
00426 sigma *= bzero2;
00427 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
00428 sigma += kr2*bonebyarg2;
00429 sigma *= damp2;
00430
00431 return sigma;
00432 }
00433
00435
00436
00437
00438
00439 G4double
00440 G4DiffuseElastic::GetDiffElasticSumProb(
00441 G4double theta
00442
00443
00444 )
00445 {
00446 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00447 G4double delta, diffuse, gamma;
00448 G4double e1, e2, bone, bone2;
00449
00450
00451
00452
00453
00454 G4double kr = fWaveVector*fNuclearRadius;
00455 G4double kr2 = kr*kr;
00456 G4double krt = kr*theta;
00457
00458 bzero = BesselJzero(krt);
00459 bzero2 = bzero*bzero;
00460 bone = BesselJone(krt);
00461 bone2 = bone*bone;
00462 bonebyarg = BesselOneByArg(krt);
00463 bonebyarg2 = bonebyarg*bonebyarg;
00464
00465 if (fParticle == theProton)
00466 {
00467 diffuse = 0.63*fermi;
00468
00469 gamma = 0.3*fermi;
00470 delta = 0.1*fermi*fermi;
00471 e1 = 0.3*fermi;
00472 e2 = 0.35*fermi;
00473 }
00474 else
00475 {
00476 diffuse = 0.63*fermi;
00477 gamma = 0.3*fermi;
00478 delta = 0.1*fermi*fermi;
00479 e1 = 0.3*fermi;
00480 e2 = 0.35*fermi;
00481 }
00482 G4double lambda = 15.;
00483
00484 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
00485
00486
00487
00488 if(fAddCoulomb)
00489 {
00490 G4double sinHalfTheta = std::sin(0.5*theta);
00491 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00492
00493 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm);
00494
00495 }
00496
00497 G4double kgamma2 = kgamma*kgamma;
00498
00499
00500
00501
00502
00503
00504 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));
00505
00506
00507
00508 damp = DampFactor(pikdt);
00509 damp2 = damp*damp;
00510
00511 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
00512 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00513
00514 sigma = kgamma2;
00515
00516 sigma *= bzero2;
00517 sigma += mode2k2*bone2;
00518 sigma += e2dk3t*bzero*bone;
00519
00520
00521 sigma += kr2*bonebyarg2;
00522
00523 sigma *= damp2;
00524
00525 return sigma;
00526 }
00527
00528
00530
00531
00532
00533
00534 G4double
00535 G4DiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
00536 {
00537 G4double theta;
00538
00539 theta = std::sqrt(alpha);
00540
00541
00542
00543 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
00544 G4double delta, diffuse, gamma;
00545 G4double e1, e2, bone, bone2;
00546
00547
00548
00549
00550
00551 G4double kr = fWaveVector*fNuclearRadius;
00552 G4double kr2 = kr*kr;
00553 G4double krt = kr*theta;
00554
00555 bzero = BesselJzero(krt);
00556 bzero2 = bzero*bzero;
00557 bone = BesselJone(krt);
00558 bone2 = bone*bone;
00559 bonebyarg = BesselOneByArg(krt);
00560 bonebyarg2 = bonebyarg*bonebyarg;
00561
00562 if (fParticle == theProton)
00563 {
00564 diffuse = 0.63*fermi;
00565
00566 gamma = 0.3*fermi;
00567 delta = 0.1*fermi*fermi;
00568 e1 = 0.3*fermi;
00569 e2 = 0.35*fermi;
00570 }
00571 else
00572 {
00573 diffuse = 0.63*fermi;
00574 gamma = 0.3*fermi;
00575 delta = 0.1*fermi*fermi;
00576 e1 = 0.3*fermi;
00577 e2 = 0.35*fermi;
00578 }
00579 G4double lambda = 15.;
00580
00581 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
00582
00583
00584
00585 if(fAddCoulomb)
00586 {
00587 G4double sinHalfTheta = theta*0.5;
00588 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
00589
00590 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm);
00591
00592 }
00593
00594 G4double kgamma2 = kgamma*kgamma;
00595
00596
00597
00598
00599
00600
00601 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));
00602
00603
00604
00605 damp = DampFactor(pikdt);
00606 damp2 = damp*damp;
00607
00608 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
00609 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
00610
00611 sigma = kgamma2;
00612
00613 sigma *= bzero2;
00614 sigma += mode2k2*bone2;
00615 sigma += e2dk3t*bzero*bone;
00616
00617
00618 sigma += kr2*bonebyarg2;
00619
00620 sigma *= damp2;
00621
00622 return sigma;
00623 }
00624
00625
00627
00628
00629
00630 G4double
00631 G4DiffuseElastic::GetIntegrandFunction( G4double alpha )
00632 {
00633 G4double result;
00634
00635 result = GetDiffElasticSumProbA(alpha);
00636
00637
00638
00639 return result;
00640 }
00641
00643
00644
00645
00646 G4double
00647 G4DiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle,
00648 G4double theta,
00649 G4double momentum,
00650 G4double A )
00651 {
00652 G4double result;
00653 fParticle = particle;
00654 fWaveVector = momentum/hbarc;
00655 fAtomicWeight = A;
00656
00657 fNuclearRadius = CalculateNuclearRad(A);
00658
00659
00660 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00661
00662
00663 result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
00664
00665 return result;
00666 }
00667
00669
00670
00671
00672 G4double G4DiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
00673 {
00674 G4double theta = SampleThetaCMS( aParticle, p, A);
00675 G4double t = 2*p*p*( 1 - std::cos(theta) );
00676 return t;
00677 }
00678
00680
00681
00682
00683
00684 G4double
00685 G4DiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle,
00686 G4double momentum, G4double A)
00687 {
00688 G4int i, iMax = 100;
00689 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
00690
00691 fParticle = particle;
00692 fWaveVector = momentum/hbarc;
00693 fAtomicWeight = A;
00694
00695 fNuclearRadius = CalculateNuclearRad(A);
00696
00697 thetaMax = 10.174/fWaveVector/fNuclearRadius;
00698
00699 if (thetaMax > pi) thetaMax = pi;
00700
00701 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00702
00703
00704 norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
00705
00706 norm *= G4UniformRand();
00707
00708 for(i = 1; i <= iMax; i++)
00709 {
00710 theta1 = (i-1)*thetaMax/iMax;
00711 theta2 = i*thetaMax/iMax;
00712 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
00713
00714 if ( sum >= norm )
00715 {
00716 result = 0.5*(theta1 + theta2);
00717 break;
00718 }
00719 }
00720 if (i > iMax ) result = 0.5*(theta1 + theta2);
00721
00722 G4double sigma = pi*thetaMax/iMax;
00723
00724 result += G4RandGauss::shoot(0.,sigma);
00725
00726 if(result < 0.) result = 0.;
00727 if(result > thetaMax) result = thetaMax;
00728
00729 return result;
00730 }
00731
00735
00736
00737
00738 G4double G4DiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
00739 G4int Z, G4int A)
00740 {
00741 fParticle = aParticle;
00742 G4double m1 = fParticle->GetPDGMass();
00743 G4double totElab = std::sqrt(m1*m1+p*p);
00744 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00745 G4LorentzVector lv1(p,0.0,0.0,totElab);
00746 G4LorentzVector lv(0.0,0.0,0.0,mass2);
00747 lv += lv1;
00748
00749 G4ThreeVector bst = lv.boostVector();
00750 lv1.boost(-bst);
00751
00752 G4ThreeVector p1 = lv1.vect();
00753 G4double momentumCMS = p1.mag();
00754
00755 G4double t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) );
00756 return t;
00757 }
00758
00760
00761
00762
00763 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
00764 G4double Z, G4double A)
00765 {
00766 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A);
00767
00768 G4double t = p*p*alpha;
00769 return t;
00770 }
00771
00773
00774
00775
00776
00777 G4double
00778 G4DiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
00779 G4double momentum, G4double Z, G4double A)
00780 {
00781 size_t iElement;
00782 G4int iMomentum, iAngle;
00783 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
00784 G4double m1 = particle->GetPDGMass();
00785
00786 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
00787 {
00788 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
00789 }
00790 if ( iElement == fElementNumberVector.size() )
00791 {
00792 InitialiseOnFly(Z,A);
00793
00794
00795
00796
00797
00798
00799 }
00800
00801
00802 fAngleTable = fAngleBank[iElement];
00803
00804 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
00805
00806 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
00807 {
00808 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
00809 }
00810 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
00811 if ( iMomentum < 0 ) iMomentum = 0;
00812
00813
00814
00815 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
00816 {
00817 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
00818
00819
00820
00821 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00822 {
00823 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00824 }
00825 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00826
00827
00828
00829 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
00830
00831
00832 }
00833 else
00834 {
00835
00836 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
00837
00838
00839
00840 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00841 {
00842
00843 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00844 }
00845 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00846
00847
00848
00849 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
00850
00851
00852 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00853
00854
00855
00856 iMomentum--;
00857
00858
00859
00860
00861
00862 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
00863 {
00864
00865 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
00866 }
00867 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
00868
00869 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
00870
00871
00872
00873 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
00874
00875
00876
00877 W = 1.0/(E2 - E1);
00878 W1 = (E2 - kinE)*W;
00879 W2 = (kinE - E1)*W;
00880
00881 randAngle = W1*theta1 + W2*theta2;
00882
00883
00884
00885 }
00886
00887
00888
00889 return randAngle;
00890 }
00891
00893
00894
00895
00896 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
00897 {
00898 fAtomicNumber = Z;
00899 fAtomicWeight = A;
00900
00901 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
00902
00903 if( verboseLevel > 0 )
00904 {
00905 G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "
00906 <<Z<<"; and A = "<<A<<G4endl;
00907 }
00908 fElementNumberVector.push_back(fAtomicNumber);
00909
00910 BuildAngleTable();
00911
00912 fAngleBank.push_back(fAngleTable);
00913
00914 return;
00915 }
00916
00918
00919
00920
00921
00922 void G4DiffuseElastic::BuildAngleTable()
00923 {
00924 G4int i, j;
00925 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
00926 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
00927
00928 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
00929
00930 fAngleTable = new G4PhysicsTable(fEnergyBin);
00931
00932 for( i = 0; i < fEnergyBin; i++)
00933 {
00934 kinE = fEnergyVector->GetLowEdgeEnergy(i);
00935 partMom = std::sqrt( kinE*(kinE + 2*m1) );
00936
00937 fWaveVector = partMom/hbarc;
00938
00939 G4double kR = fWaveVector*fNuclearRadius;
00940 G4double kR2 = kR*kR;
00941 G4double kRmax = 18.6;
00942 G4double kRcoul = 1.9;
00943
00944
00945
00946 alphaMax = kRmax*kRmax/kR2;
00947
00948 if (alphaMax > 4.) alphaMax = 4.;
00949
00950 alphaCoulomb = kRcoul*kRcoul/kR2;
00951
00952 if( z )
00953 {
00954 a = partMom/m1;
00955 fBeta = a/std::sqrt(1+a*a);
00956 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
00957 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
00958 }
00959 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
00960
00961
00962
00963 G4double delth = alphaMax/fAngleBin;
00964
00965 sum = 0.;
00966
00967
00968 fAddCoulomb = true;
00969
00970
00971 for(j = fAngleBin-1; j >= 1; j--)
00972 {
00973
00974
00975
00976
00977
00978
00979 alpha1 = delth*(j-1);
00980
00981 alpha2 = alpha1 + delth;
00982
00983
00984 if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
00985
00986 delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
00987
00988
00989 sum += delta;
00990
00991 angleVector->PutValue( j-1 , alpha1, sum );
00992
00993 }
00994 fAngleTable->insertAt(i,angleVector);
00995
00996
00997
00998 }
00999 return;
01000 }
01001
01003
01004
01005
01006 G4double
01007 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
01008 {
01009 G4double x1, x2, y1, y2, randAngle;
01010
01011 if( iAngle == 0 )
01012 {
01013 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01014
01015 }
01016 else
01017 {
01018 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
01019 {
01020 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
01021 }
01022 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
01023 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
01024
01025 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
01026 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
01027
01028 if ( x1 == x2 ) randAngle = x2;
01029 else
01030 {
01031 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
01032 else
01033 {
01034 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
01035 }
01036 }
01037 }
01038 return randAngle;
01039 }
01040
01041
01042
01044
01045
01046
01047
01048
01049 G4double
01050 G4DiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
01051 G4double tmass, G4double A)
01052 {
01053 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01054 G4double m1 = theParticle->GetPDGMass();
01055 G4double plab = aParticle->GetTotalMomentum();
01056 G4LorentzVector lv1 = aParticle->Get4Momentum();
01057 G4LorentzVector lv(0.0,0.0,0.0,tmass);
01058 lv += lv1;
01059
01060 G4ThreeVector bst = lv.boostVector();
01061 lv1.boost(-bst);
01062
01063 G4ThreeVector p1 = lv1.vect();
01064 G4double ptot = p1.mag();
01065 G4double tmax = 4.0*ptot*ptot;
01066 G4double t = 0.0;
01067
01068
01069
01070
01071
01072
01073 t = SampleT( theParticle, ptot, A);
01074
01075
01076 if(!(t < 0.0 || t >= 0.0))
01077 {
01078 if (verboseLevel > 0)
01079 {
01080 G4cout << "G4DiffuseElastic:WARNING: A = " << A
01081 << " mom(GeV)= " << plab/GeV
01082 << " S-wave will be sampled"
01083 << G4endl;
01084 }
01085 t = G4UniformRand()*tmax;
01086 }
01087 if(verboseLevel>1)
01088 {
01089 G4cout <<" t= " << t << " tmax= " << tmax
01090 << " ptot= " << ptot << G4endl;
01091 }
01092
01093
01094 G4double phi = G4UniformRand()*twopi;
01095 G4double cost = 1. - 2.0*t/tmax;
01096 G4double sint;
01097
01098 if( cost >= 1.0 )
01099 {
01100 cost = 1.0;
01101 sint = 0.0;
01102 }
01103 else if( cost <= -1.0)
01104 {
01105 cost = -1.0;
01106 sint = 0.0;
01107 }
01108 else
01109 {
01110 sint = std::sqrt((1.0-cost)*(1.0+cost));
01111 }
01112 if (verboseLevel>1)
01113 {
01114 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
01115 }
01116 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01117 v1 *= ptot;
01118 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01119
01120 nlv1.boost(bst);
01121
01122 G4ThreeVector np1 = nlv1.vect();
01123
01124
01125
01126 G4double theta = np1.theta();
01127
01128 return theta;
01129 }
01130
01132
01133
01134
01135
01136
01137 G4double
01138 G4DiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
01139 G4double tmass, G4double thetaCMS)
01140 {
01141 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01142 G4double m1 = theParticle->GetPDGMass();
01143
01144 G4LorentzVector lv1 = aParticle->Get4Momentum();
01145 G4LorentzVector lv(0.0,0.0,0.0,tmass);
01146
01147 lv += lv1;
01148
01149 G4ThreeVector bst = lv.boostVector();
01150
01151 lv1.boost(-bst);
01152
01153 G4ThreeVector p1 = lv1.vect();
01154 G4double ptot = p1.mag();
01155
01156 G4double phi = G4UniformRand()*twopi;
01157 G4double cost = std::cos(thetaCMS);
01158 G4double sint;
01159
01160 if( cost >= 1.0 )
01161 {
01162 cost = 1.0;
01163 sint = 0.0;
01164 }
01165 else if( cost <= -1.0)
01166 {
01167 cost = -1.0;
01168 sint = 0.0;
01169 }
01170 else
01171 {
01172 sint = std::sqrt((1.0-cost)*(1.0+cost));
01173 }
01174 if (verboseLevel>1)
01175 {
01176 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
01177 }
01178 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01179 v1 *= ptot;
01180 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
01181
01182 nlv1.boost(bst);
01183
01184 G4ThreeVector np1 = nlv1.vect();
01185
01186
01187 G4double thetaLab = np1.theta();
01188
01189 return thetaLab;
01190 }
01192
01193
01194
01195
01196
01197 G4double
01198 G4DiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
01199 G4double tmass, G4double thetaLab)
01200 {
01201 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
01202 G4double m1 = theParticle->GetPDGMass();
01203 G4double plab = aParticle->GetTotalMomentum();
01204 G4LorentzVector lv1 = aParticle->Get4Momentum();
01205 G4LorentzVector lv(0.0,0.0,0.0,tmass);
01206
01207 lv += lv1;
01208
01209 G4ThreeVector bst = lv.boostVector();
01210
01211
01212
01213
01214
01215
01216 G4double phi = G4UniformRand()*twopi;
01217 G4double cost = std::cos(thetaLab);
01218 G4double sint;
01219
01220 if( cost >= 1.0 )
01221 {
01222 cost = 1.0;
01223 sint = 0.0;
01224 }
01225 else if( cost <= -1.0)
01226 {
01227 cost = -1.0;
01228 sint = 0.0;
01229 }
01230 else
01231 {
01232 sint = std::sqrt((1.0-cost)*(1.0+cost));
01233 }
01234 if (verboseLevel>1)
01235 {
01236 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
01237 }
01238 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
01239 v1 *= plab;
01240 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
01241
01242 nlv1.boost(-bst);
01243
01244 G4ThreeVector np1 = nlv1.vect();
01245
01246
01247 G4double thetaCMS = np1.theta();
01248
01249 return thetaCMS;
01250 }
01251
01253
01254
01255
01256
01257 void G4DiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
01258 G4double Z, G4double A)
01259 {
01260 fAtomicNumber = Z;
01261 fAtomicWeight = A;
01262 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
01263
01264
01265
01266 G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
01267 <<Z<<"; and A = "<<A<<G4endl;
01268
01269 fElementNumberVector.push_back(fAtomicNumber);
01270
01271
01272
01273
01274 G4int i=0, j;
01275 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
01276 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
01277 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
01278 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
01279 G4double epsilon = 0.001;
01280
01281 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
01282
01283 fAngleTable = new G4PhysicsTable(fEnergyBin);
01284
01285 fWaveVector = partMom/hbarc;
01286
01287 G4double kR = fWaveVector*fNuclearRadius;
01288 G4double kR2 = kR*kR;
01289 G4double kRmax = 10.6;
01290 G4double kRcoul = 1.2;
01291
01292 alphaMax = kRmax*kRmax/kR2;
01293
01294 if (alphaMax > 4.) alphaMax = 4.;
01295
01296 alphaCoulomb = kRcoul*kRcoul/kR2;
01297
01298 if( z )
01299 {
01300 a = partMom/m1;
01301 fBeta = a/std::sqrt(1+a*a);
01302 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
01303 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
01304 }
01305 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
01306
01307
01308
01309
01310 fAddCoulomb = false;
01311
01312 for(j = 1; j < fAngleBin; j++)
01313 {
01314
01315
01316
01317 alpha1 = alphaMax*(j-1)/fAngleBin;
01318 alpha2 = alphaMax*( j )/fAngleBin;
01319
01320 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
01321
01322 deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01323 deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
01324 deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
01325 alpha1, alpha2,epsilon);
01326
01327
01328
01329
01330 sumL10 += deltaL10;
01331 sumL96 += deltaL96;
01332 sumAG += deltaAG;
01333
01334 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
01335 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
01336
01337 angleVector->PutValue( j-1 , alpha1, sumL10 );
01338 }
01339 fAngleTable->insertAt(i,angleVector);
01340 fAngleBank.push_back(fAngleTable);
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 return;
01354 }
01355
01356
01357