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