#include <G4AntiNuclElastic.hh>
Inheritance diagram for G4AntiNuclElastic:
Public Member Functions | |
G4AntiNuclElastic () | |
virtual | ~G4AntiNuclElastic () |
virtual G4double | SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) |
G4double | SampleThetaCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) |
G4double | SampleThetaLab (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) |
G4double | CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum) |
G4double | CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2) |
G4double | CalculateAm (G4double momentum, G4double n, G4double Z) |
G4double | DampFactor (G4double z) |
G4double | BesselJzero (G4double z) |
G4double | BesselJone (G4double z) |
G4double | BesselOneByArg (G4double z) |
G4double | GetcosTeta1 (G4double plab, G4int A) |
G4ComponentAntiNuclNuclearXS * | GetComponentCrossSection () |
Definition at line 50 of file G4AntiNuclElastic.hh.
G4AntiNuclElastic::G4AntiNuclElastic | ( | ) |
Definition at line 55 of file G4AntiNuclElastic.cc.
References G4Alpha::Alpha(), G4AntiAlpha::AntiAlpha(), G4AntiDeuteron::AntiDeuteron(), G4AntiHe3::AntiHe3(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiTriton::AntiTriton(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), and G4Proton::Proton().
00056 : G4HadronElastic("AntiAElastic") 00057 { 00058 //V.Ivanchenko commented out 00059 //SetMinEnergy( 0.1*GeV ); 00060 //SetMaxEnergy( 10.*TeV ); 00061 00062 00063 theAProton = G4AntiProton::AntiProton(); 00064 theANeutron = G4AntiNeutron::AntiNeutron(); 00065 theADeuteron = G4AntiDeuteron::AntiDeuteron(); 00066 theATriton = G4AntiTriton::AntiTriton(); 00067 theAAlpha = G4AntiAlpha::AntiAlpha(); 00068 theAHe3 = G4AntiHe3::AntiHe3(); 00069 00070 theProton = G4Proton::Proton(); 00071 theNeutron = G4Neutron::Neutron(); 00072 theDeuteron = G4Deuteron::Deuteron(); 00073 theAlpha = G4Alpha::Alpha(); 00074 00075 00076 cs = new G4ComponentAntiNuclNuclearXS(); 00077 fParticle = 0; 00078 fWaveVector = 0.; 00079 fBeta = 0.; 00080 fZommerfeld = 0.; 00081 fAm = 0.; 00082 fTetaCMS = 0.; 00083 fRa = 0.; 00084 fRef = 0.; 00085 fceff = 0.; 00086 fptot = 0.; 00087 fTmax = 0.; 00088 fThetaLab = 0.; 00089 }
G4AntiNuclElastic::~G4AntiNuclElastic | ( | ) | [virtual] |
Definition at line 583 of file G4AntiNuclElastic.cc.
Referenced by BesselOneByArg().
00584 { 00585 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 00586 00587 modvalue = std::fabs(value); 00588 00589 if ( modvalue < 8.0 ) 00590 { 00591 value2 = value*value; 00592 fact1 = value*(72362614232.0 + value2*(-7895059235.0 00593 + value2*( 242396853.1 00594 + value2*(-2972611.439 00595 + value2*( 15704.48260 00596 + value2*(-30.16036606 ) ) ) ) ) ); 00597 00598 fact2 = 144725228442.0 + value2*(2300535178.0 00599 + value2*(18583304.74 00600 + value2*(99447.43394 00601 + value2*(376.9991397 00602 + value2*1.0 ) ) ) ); 00603 bessel = fact1/fact2; 00604 } 00605 else 00606 { 00607 arg = 8.0/modvalue; 00608 value2 = arg*arg; 00609 00610 shift = modvalue - 2.356194491; 00611 00612 fact1 = 1.0 + value2*( 0.183105e-2 00613 + value2*(-0.3516396496e-4 00614 + value2*(0.2457520174e-5 00615 + value2*(-0.240337019e-6 ) ) ) ); 00616 00617 fact2 = 0.04687499995 + value2*(-0.2002690873e-3 00618 + value2*( 0.8449199096e-5 00619 + value2*(-0.88228987e-6 00620 + value2*0.105787412e-6 ) ) ); 00621 00622 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 00623 if (value < 0.0) bessel = -bessel; 00624 } 00625 return bessel; 00626 }
Definition at line 532 of file G4AntiNuclElastic.cc.
Referenced by SampleInvariantT().
00533 { 00534 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 00535 00536 modvalue = std::fabs(value); 00537 00538 if ( value < 8.0 && value > -8.0 ) 00539 { 00540 value2 = value*value; 00541 00542 fact1 = 57568490574.0 + value2*(-13362590354.0 00543 + value2*( 651619640.7 00544 + value2*(-11214424.18 00545 + value2*( 77392.33017 00546 + value2*(-184.9052456 ) ) ) ) ); 00547 00548 fact2 = 57568490411.0 + value2*( 1029532985.0 00549 + value2*( 9494680.718 00550 + value2*(59272.64853 00551 + value2*(267.8532712 00552 + value2*1.0 ) ) ) ); 00553 00554 bessel = fact1/fact2; 00555 } 00556 else 00557 { 00558 arg = 8.0/modvalue; 00559 00560 value2 = arg*arg; 00561 00562 shift = modvalue-0.785398164; 00563 00564 fact1 = 1.0 + value2*(-0.1098628627e-2 00565 + value2*(0.2734510407e-4 00566 + value2*(-0.2073370639e-5 00567 + value2*0.2093887211e-6 ) ) ); 00568 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3 00569 + value2*(-0.6911147651e-5 00570 + value2*(0.7621095161e-6 00571 - value2*0.934945152e-7 ) ) ); 00572 00573 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 00574 } 00575 return bessel; 00576 }
Definition at line 630 of file G4AntiNuclElastic.cc.
References BesselJone().
Referenced by SampleInvariantT().
00631 { 00632 G4double x2, result; 00633 00634 if( std::fabs(x) < 0.01 ) 00635 { 00636 x *= 0.5; 00637 x2 = x*x; 00638 result = (2.- x2 + x2*x2/6.)/4.; 00639 } 00640 else 00641 { 00642 result = BesselJone(x)/x; 00643 } 00644 return result; 00645 }
Definition at line 516 of file G4AntiNuclElastic.cc.
References G4Pow::A13(), and G4Pow::GetInstance().
Referenced by SampleInvariantT().
00517 { 00518 G4double k = momentum/hbarc; 00519 G4double ch = 1.13 + 3.76*n*n; 00520 G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius; 00521 G4double zn2 = zn*zn; 00522 fAm = ch/zn2; 00523 00524 return fAm; 00525 }
G4double G4AntiNuclElastic::CalculateParticleBeta | ( | const G4ParticleDefinition * | particle, | |
G4double | momentum | |||
) |
Definition at line 493 of file G4AntiNuclElastic.cc.
References G4ParticleDefinition::GetPDGMass().
Referenced by SampleInvariantT().
00495 { 00496 G4double mass = particle->GetPDGMass(); 00497 G4double a = momentum/mass; 00498 fBeta = a/std::sqrt(1+a*a); 00499 00500 return fBeta; 00501 }
Definition at line 507 of file G4AntiNuclElastic.cc.
Referenced by SampleInvariantT().
00508 { 00509 fZommerfeld = fine_structure_const*Z1*Z2/beta; 00510 00511 return fZommerfeld; 00512 }
Definition at line 473 of file G4AntiNuclElastic.cc.
Referenced by SampleInvariantT().
00474 { 00475 G4double df; 00476 G4double f3 = 6.; // first factorials 00477 00478 if( std::fabs(x) < 0.01 ) 00479 { 00480 df=1./(1.+x*x/f3); 00481 } 00482 else 00483 { 00484 df = x/std::sinh(x); 00485 } 00486 return df; 00487 }
G4ComponentAntiNuclNuclearXS * G4AntiNuclElastic::GetComponentCrossSection | ( | ) | [inline] |
Definition at line 126 of file G4AntiNuclElastic.hh.
Referenced by G4HadronQElasticPhysics::ConstructProcess(), G4HadronHElasticPhysics::ConstructProcess(), G4HadronElasticPhysics::ConstructProcess(), and G4HadronDElasticPhysics::ConstructProcess().
Definition at line 649 of file G4AntiNuclElastic.cc.
References G4Pow::GetInstance(), and G4Pow::Z23().
Referenced by SampleInvariantT().
00650 { 00651 00652 // G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi; 00653 G4double p0 = 1.*hbarc/fermi; 00654 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab); 00655 G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab); 00657 if(cteta1 < -1.) cteta1 = -1.0; 00658 return cteta1; 00659 }
G4double G4AntiNuclElastic::SampleInvariantT | ( | const G4ParticleDefinition * | p, | |
G4double | plab, | |||
G4int | Z, | |||
G4int | A | |||
) | [virtual] |
Reimplemented from G4HadronElastic.
Definition at line 99 of file G4AntiNuclElastic.cc.
References BesselJzero(), BesselOneByArg(), CalculateAm(), CalculateParticleBeta(), CalculateZommerfeld(), DampFactor(), G4UniformRand, G4ComponentAntiNuclNuclearXS::GetAntiHadronNucleonTotCrSc(), G4ParticleDefinition::GetBaryonNumber(), GetcosTeta1(), G4ComponentAntiNuclNuclearXS::GetElasticElementCrossSection(), G4Pow::GetInstance(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ComponentAntiNuclNuclearXS::GetTotalElementCrossSection(), G4He3::He3(), CLHEP::detail::n, G4INCL::Math::pi, sqr(), G4Triton::Triton(), and G4Pow::Z13().
Referenced by SampleThetaCMS(), and SampleThetaLab().
00101 { 00102 G4double T; 00103 G4double Mproj = particle->GetPDGMass(); 00104 G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj)); 00105 G4double ctet1 = GetcosTeta1(Plab, A); 00106 00107 G4double energy=Pproj.e()-Mproj; 00108 00109 const G4ParticleDefinition* theParticle = particle; 00110 00111 G4ParticleDefinition * theDef = 0; 00112 00113 if(Z == 1 && A == 1) theDef = theProton; 00114 else if (Z == 1 && A == 2) theDef = theDeuteron; 00115 else if (Z == 1 && A == 3) theDef = G4Triton::Triton(); 00116 else if (Z == 2 && A == 3) theDef = G4He3::He3(); 00117 else if (Z == 2 && A == 4) theDef = theAlpha; 00118 00119 00120 G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z); 00121 00122 //transform to CMS 00123 00124 G4LorentzVector lv(0.0,0.0,0.0,TargMass); 00125 lv += Pproj; 00126 G4double S = lv.mag2()/GeV/GeV; 00127 00128 G4ThreeVector bst = lv.boostVector(); 00129 Pproj.boost(-bst); 00130 00131 G4ThreeVector p1 = Pproj.vect(); 00132 G4double ptot = p1.mag(); 00133 00134 fbst = bst; 00135 fptot= ptot; 00136 fTmax = 4.0*ptot*ptot; 00137 00138 if(Plab/std::abs(particle->GetBaryonNumber()) < 100.*MeV) // Uzhi 24 Nov. 2011 00139 {return fTmax*G4UniformRand();} // Uzhi 24 Nov. 2011 00140 00141 G4double Z1 = particle->GetPDGCharge(); 00142 G4double Z2 = Z; 00143 00144 G4double beta = CalculateParticleBeta(particle, ptot); 00145 G4double n = CalculateZommerfeld( beta, Z1, Z2 ); 00146 G4double Am = CalculateAm( ptot, n, Z2 ); 00147 fWaveVector = ptot; // /hbarc; 00148 00149 G4LorentzVector Fproj(0.,0.,0.,0.); 00150 G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1); 00151 XsCoulomb=XsCoulomb*0.38938e+6; 00152 00153 00154 G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A); 00155 G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A); 00156 00157 00158 XsElastHad/=millibarn; XstotalHad/=millibarn; 00159 00160 00161 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad); 00162 00163 // G4cout<<" XselastHadron " << XsElastHad << " XsCol "<< XsCoulomb <<G4endl; 00164 // G4cout <<" XsTotal" << XstotalHad <<G4endl; 00165 // G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl; 00166 00167 if(G4UniformRand() < CoulombProb) 00168 { // Simulation of Coulomb scattering 00169 00170 G4double phi = twopi * G4UniformRand(); 00171 G4double Ksi = G4UniformRand(); 00172 00173 G4double par1 = 2.*(1.+Am)/(1.+ctet1); 00174 00175 // ////sample ThetaCMS in Coulomb part 00176 00177 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi); 00178 00179 G4double PtZ=ptot*cosThetaCMS; 00180 Fproj.setPz(PtZ); 00181 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS); 00182 G4double PtX= PtProjCMS * std::cos(phi); 00183 G4double PtY= PtProjCMS * std::sin(phi); 00184 Fproj.setPx(PtX); 00185 Fproj.setPy(PtY); 00186 Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj)); 00187 T = -(Pproj-Fproj).mag2(); 00188 } else 00189 00190 { 00192 00193 // G4double Qmax = 2.*ptot*197.33; // in fm^-1 00194 G4double Qmax = 2.*3.0*197.33; // in fm^-1 00195 G4double Amag = 70*70; // A1 in Magora funct:A1*exp(-q*A2) 00196 G4double SlopeMag = 2.*3.0; // A2 in Magora funct:A1*exp(-q*A2) 00197 00198 G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy); 00199 00200 fRa = 1.113*G4Pow::GetInstance()->Z13(A) - 00201 0.227/G4Pow::GetInstance()->Z13(A); 00202 if(A == 3) fRa=1.81; 00203 if(A == 4) fRa=1.37; 00204 00205 if((A>=12.) && (A<27) ) fRa=fRa*0.85; 00206 if((A>=27.) && (A<48) ) fRa=fRa*0.90; 00207 if((A>=48.) && (A<65) ) fRa=fRa*0.95; 00208 00209 G4double Ref2 = 0; 00210 G4double ceff2 =0; 00211 G4double rho = 0; 00212 if ((theParticle == theAProton) || (theParticle == theANeutron)) 00213 { 00214 if(theDef == theProton) 00215 { 00216 // G4double Mp2=sqr(theDef->GetPDGMass()/GeV ); 00217 00218 // change 30 October 00219 00220 if(Plab < 610.) 00221 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.- 00222 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;} 00223 if((Plab < 5500.)&&(Plab >= 610.) ) 00224 { rho = 0.22; } 00225 if((Plab >= 5500.)&&(Plab < 12300.) ) 00226 { rho = -0.32; } 00227 if( Plab >= 12300.) 00228 { rho = 0.135-2.26/(std::sqrt(S)) ;} 00229 00230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*std::log(S) ; 00231 ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ; 00232 00233 /* 00234 Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55; 00235 if(S>1000.) Ref2=0.62+0.02*std::log(S) ; 00236 ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * std::log(S) ; 00237 if(S>1000.) ceff2 = 0.005 * std::log(S) + 0.29; 00238 */ 00239 00240 Ref2=Ref2*Ref2; 00241 ceff2 = ceff2*ceff2; 00242 00243 SlopeMag = 0.5; // Uzhi 00244 Amag= 1.; // Uzhi 00245 } 00246 00247 if(Z>2) 00248 { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa; 00249 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*std::exp(-0.03*sig_pbarp); 00250 } 00251 if( (Z==2)&&(A==4) ) 00252 { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp; 00253 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp); 00254 } 00255 if( (Z==1)&&(A==3) ) 00256 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp; 00257 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 00258 } 00259 if( (Z==2)&&(A==3) ) 00260 { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp; 00261 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 00262 } 00263 if( (Z==1)&&(A==2) ) 00264 { 00265 Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp; 00266 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp); 00267 } 00268 } 00269 00270 if (theParticle == theADeuteron) 00271 { 00272 sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.); 00273 Ref2 = XstotalHad/10./2./pi ; 00274 if(Z>2) 00275 { 00276 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * std::exp(-0.03*sig_pbarp); 00277 } 00278 if(theDef == theProton) 00279 { 00280 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*std::exp(-0.03*sig_pbarp); 00281 } 00282 if(theDef == theDeuteron) 00283 { 00284 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * std::exp(-0.03*sig_pbarp); 00285 } 00286 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) 00287 { 00288 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp); 00289 } 00290 if(theDef == theAlpha) 00291 { 00292 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp); 00293 } 00294 } 00295 00296 if( (theParticle ==theAHe3) || (theParticle ==theATriton) ) 00297 { 00298 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.); 00299 Ref2 = XstotalHad/10./2./pi ; 00300 if(Z>2) 00301 { 00302 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*std::exp(-0.03*sig_pbarp); 00303 } 00304 if(theDef == theProton) 00305 { 00306 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*std::exp(-0.03*sig_pbarp); 00307 } 00308 if(theDef == theDeuteron) 00309 { 00310 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * std::exp(-0.02*sig_pbarp); 00311 } 00312 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) 00313 { 00314 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * std::exp(-0.02*sig_pbarp); 00315 } 00316 if(theDef == theAlpha) 00317 { 00318 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp); 00319 } 00320 } 00321 00322 00323 if (theParticle == theAAlpha) 00324 { 00325 sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.); 00326 Ref2 = XstotalHad/10./2./pi ; 00327 if(Z>2) 00328 { 00329 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * std::exp(-0.03*sig_pbarp); 00330 } 00331 if(theDef == theProton) 00332 { 00333 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*std::exp(-0.03*sig_pbarp); 00334 } 00335 if(theDef == theDeuteron) 00336 { 00337 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * std::exp(-0.02*sig_pbarp); 00338 } 00339 if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) ) 00340 { 00341 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * std::exp(-0.03*sig_pbarp); 00342 } 00343 if(theDef == theAlpha) 00344 { 00345 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * std::exp(-0.03*sig_pbarp); 00346 } 00347 } 00348 00349 fRef=std::sqrt(Ref2); 00350 fceff = std::sqrt(ceff2); 00351 // G4cout<<" Ref "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl; 00352 00353 00354 G4double Q = 0.0 ; 00355 G4double BracFunct; 00356 do 00357 { 00358 Q = -std::log(1.-(1.- std::exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag; 00359 G4double x = fRef * Q; 00360 BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) ) 00361 * sqr(DampFactor(pi*fceff*Q))) /(Amag*std::exp(-SlopeMag*Q)); 00362 00363 BracFunct = BracFunct * Q * sqr(sqr(fRef)); 00364 } 00365 while (G4UniformRand()>BracFunct); 00366 00367 T= sqr(Q); 00368 T*=3.893913e+4; // fm -> MeV^2 00369 } 00370 00371 G4double cosTet=1.0-T/(2.*ptot*ptot); 00372 if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov. 00373 if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov. 00374 fTetaCMS=std::acos(cosTet); 00375 00376 return T; 00377 }
G4double G4AntiNuclElastic::SampleThetaCMS | ( | const G4ParticleDefinition * | p, | |
G4double | plab, | |||
G4int | Z, | |||
G4int | A | |||
) |
Definition at line 381 of file G4AntiNuclElastic.cc.
References G4cout, G4endl, G4UniformRand, SampleInvariantT(), and G4HadronicInteraction::verboseLevel.
00383 { 00384 G4double T; 00385 T = SampleInvariantT( p, plab, Z, A); 00386 00387 // NaN finder 00388 if(!(T < 0.0 || T >= 0.0)) 00389 { 00390 if (verboseLevel > 0) 00391 { 00392 G4cout << "G4DiffuseElastic:WARNING: A = " << A 00393 << " mom(GeV)= " << plab/GeV 00394 << " S-wave will be sampled" 00395 << G4endl; 00396 } 00397 T = G4UniformRand()*fTmax; 00398 00399 } 00400 00401 if(fptot > 0.) // Uzhi 24 Nov. 2011 00402 { 00403 G4double cosTet=1.0-T/(2.*fptot*fptot); 00404 if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov. 00405 if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov. 00406 fTetaCMS=std::acos(cosTet); 00407 return fTetaCMS; 00408 } else // Uzhi 24 Nov. 2011 00409 { // Uzhi 24 Nov. 2011 00410 return 2.*G4UniformRand()-1.; // Uzhi 24 Nov. 2011 00411 } // Uzhi 24 Nov. 2011 00412 }
G4double G4AntiNuclElastic::SampleThetaLab | ( | const G4ParticleDefinition * | p, | |
G4double | plab, | |||
G4int | Z, | |||
G4int | A | |||
) |
Definition at line 417 of file G4AntiNuclElastic.cc.
References G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), SampleInvariantT(), and G4HadronicInteraction::verboseLevel.
00419 { 00420 G4double T; 00421 T = SampleInvariantT( p, plab, Z, A); 00422 00423 // NaN finder 00424 if(!(T < 0.0 || T >= 0.0)) 00425 { 00426 if (verboseLevel > 0) 00427 { 00428 G4cout << "G4DiffuseElastic:WARNING: A = " << A 00429 << " mom(GeV)= " << plab/GeV 00430 << " S-wave will be sampled" 00431 << G4endl; 00432 } 00433 T = G4UniformRand()*fTmax; 00434 } 00435 00436 G4double phi = G4UniformRand()*twopi; 00437 00438 G4double cost(1.); 00439 if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} // Uzhi 24 Nov. 2011 00440 00441 G4double sint; 00442 if( cost >= 1.0 ) 00443 { 00444 cost = 1.0; 00445 sint = 0.0; 00446 } 00447 else if( cost <= -1.0) 00448 { 00449 cost = -1.0; 00450 sint = 0.0; 00451 } 00452 else 00453 { 00454 sint = std::sqrt((1.0-cost)*(1.0+cost)); 00455 } 00456 00457 G4double m1 = p->GetPDGMass(); 00458 G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost); 00459 v *= fptot; 00460 G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1)); 00461 00462 nlv.boost(fbst); 00463 00464 G4ThreeVector np = nlv.vect(); 00465 G4double theta = np.theta(); 00466 fThetaLab = theta; 00467 00468 return theta; 00469 }