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 <iostream>
00038
00039 #include "G4LElastic.hh"
00040 #include "globals.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4SystemOfUnits.hh"
00043 #include "Randomize.hh"
00044 #include "G4ParticleTable.hh"
00045 #include "G4IonTable.hh"
00046
00047 G4LElastic::G4LElastic(const G4String& name)
00048 :G4HadronicInteraction(name)
00049 {
00050 SetMinEnergy(0.0);
00051 SetMaxEnergy(DBL_MAX);
00052 }
00053
00054
00055 void G4LElastic::ModelDescription(std::ostream& outFile) const
00056 {
00057 outFile << "G4LElastic is one of the Low Energy Parameterized (LEP)\n"
00058 << "models used to implement elastic hadron scattering from nuclei.\n"
00059 << "It is a re-engineered version of the GHEISHA code of\n"
00060 << "H. Fesefeldt. It performs simplified two-body elastic\n"
00061 << "scattering for all long-lived hadronic projectiles by using\n"
00062 << "a two-exponential parameterization in momentum transfer.\n"
00063 << "It is valid for incident hadrons of all energies.\n";
00064 }
00065
00066
00067 G4HadFinalState*
00068 G4LElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00069 {
00070 if(getenv("debug_LElastic")) verboseLevel = 5;
00071 theParticleChange.Clear();
00072 const G4HadProjectile* aParticle = &aTrack;
00073 G4double atno2 = targetNucleus.GetA_asInt();
00074 G4double zTarget = targetNucleus.GetZ_asInt();
00075 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00076 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00077
00078
00079
00080 G4DynamicParticle* aSecondary = 0;
00081 if (atno2 < 1.5) {
00082 const G4ParticleDefinition* aParticleType = aParticle->GetDefinition();
00083 if (aParticleType == G4PionPlus::PionPlus())
00084 aSecondary = LightMedia.PionPlusExchange(aParticle, targetNucleus);
00085 else if (aParticleType == G4PionMinus::PionMinus())
00086 aSecondary = LightMedia.PionMinusExchange(aParticle, targetNucleus);
00087 else if (aParticleType == G4KaonPlus::KaonPlus())
00088 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus);
00089 else if (aParticleType == G4KaonZeroShort::KaonZeroShort())
00090 aSecondary = LightMedia.KaonZeroShortExchange(aParticle,targetNucleus);
00091 else if (aParticleType == G4KaonZeroLong::KaonZeroLong())
00092 aSecondary = LightMedia.KaonZeroLongExchange(aParticle, targetNucleus);
00093 else if (aParticleType == G4KaonMinus::KaonMinus())
00094 aSecondary = LightMedia.KaonMinusExchange(aParticle, targetNucleus);
00095 else if (aParticleType == G4Proton::Proton())
00096 aSecondary = LightMedia.ProtonExchange(aParticle, targetNucleus);
00097 else if (aParticleType == G4AntiProton::AntiProton())
00098 aSecondary = LightMedia.AntiProtonExchange(aParticle, targetNucleus);
00099 else if (aParticleType == G4Neutron::Neutron())
00100 aSecondary = LightMedia.NeutronExchange(aParticle, targetNucleus);
00101 else if (aParticleType == G4AntiNeutron::AntiNeutron())
00102 aSecondary = LightMedia.AntiNeutronExchange(aParticle, targetNucleus);
00103 else if (aParticleType == G4Lambda::Lambda())
00104 aSecondary = LightMedia.LambdaExchange(aParticle, targetNucleus);
00105 else if (aParticleType == G4AntiLambda::AntiLambda())
00106 aSecondary = LightMedia.AntiLambdaExchange(aParticle, targetNucleus);
00107 else if (aParticleType == G4SigmaPlus::SigmaPlus())
00108 aSecondary = LightMedia.SigmaPlusExchange(aParticle, targetNucleus);
00109 else if (aParticleType == G4SigmaMinus::SigmaMinus())
00110 aSecondary = LightMedia.SigmaMinusExchange(aParticle, targetNucleus);
00111 else if (aParticleType == G4AntiSigmaPlus::AntiSigmaPlus())
00112 aSecondary = LightMedia.AntiSigmaPlusExchange(aParticle,targetNucleus);
00113 else if (aParticleType == G4AntiSigmaMinus::AntiSigmaMinus())
00114 aSecondary= LightMedia.AntiSigmaMinusExchange(aParticle,targetNucleus);
00115 else if (aParticleType == G4XiZero::XiZero())
00116 aSecondary = LightMedia.XiZeroExchange(aParticle, targetNucleus);
00117 else if (aParticleType == G4XiMinus::XiMinus())
00118 aSecondary = LightMedia.XiMinusExchange(aParticle, targetNucleus);
00119 else if (aParticleType == G4AntiXiZero::AntiXiZero())
00120 aSecondary = LightMedia.AntiXiZeroExchange(aParticle, targetNucleus);
00121 else if (aParticleType == G4AntiXiMinus::AntiXiMinus())
00122 aSecondary = LightMedia.AntiXiMinusExchange(aParticle, targetNucleus);
00123 else if (aParticleType == G4OmegaMinus::OmegaMinus())
00124 aSecondary = LightMedia.OmegaMinusExchange(aParticle, targetNucleus);
00125 else if (aParticleType == G4AntiOmegaMinus::AntiOmegaMinus())
00126 aSecondary= LightMedia.AntiOmegaMinusExchange(aParticle,targetNucleus);
00127 else if (aParticleType == G4KaonPlus::KaonPlus())
00128 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus);
00129 }
00130
00131
00132 if (aSecondary) {
00133 aSecondary->SetMomentum(aParticle->Get4Momentum().vect());
00134 theParticleChange.SetStatusChange(stopAndKill);
00135 theParticleChange.AddSecondary(aSecondary);
00136 }
00137
00138
00139 G4double p = aParticle->GetTotalMomentum()/GeV;
00140 if (verboseLevel > 1)
00141 G4cout << "G4LElastic::DoIt: Incident particle p=" << p << " GeV" << G4endl;
00142
00143 if (p < 0.01) return &theParticleChange;
00144
00145
00146
00147
00148
00149
00150 G4double ran = G4UniformRand();
00151 G4double aa, bb, cc, dd, rr;
00152 if (atno2 <= 62.) {
00153 aa = std::pow(atno2, 1.63);
00154 bb = 14.5*std::pow(atno2, 0.66);
00155 cc = 1.4*std::pow(atno2, 0.33);
00156 dd = 10.;
00157 }
00158 else {
00159 aa = std::pow(atno2, 1.33);
00160 bb = 60.*std::pow(atno2, 0.33);
00161 cc = 0.4*std::pow(atno2, 0.40);
00162 dd = 10.;
00163 }
00164 aa = aa/bb;
00165 cc = cc/dd;
00166 rr = (aa + cc)*ran;
00167 if (verboseLevel > 1) {
00168 G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
00169 G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
00170 }
00171 G4double t1 = -std::log(ran)/bb;
00172 G4double t2 = -std::log(ran)/dd;
00173 if (verboseLevel > 1) {
00174 G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) <<
00175 G4endl;
00176 G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) <<
00177 G4endl;
00178 }
00179 G4double eps = 0.001;
00180 G4int ind1 = 10;
00181 G4double t;
00182 G4int ier1;
00183 ier1 = Rtmi(&t, t1, t2, eps, ind1,
00184 aa, bb, cc, dd, rr);
00185 if (verboseLevel > 1) {
00186 G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
00187 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) <<
00188 G4endl;
00189 }
00190 if (ier1 != 0) t = 0.25*(3.*t1 + t2);
00191 if (verboseLevel > 1) {
00192 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) <<
00193 G4endl;
00194 }
00195 G4double phi = G4UniformRand()*twopi;
00196 rr = 0.5*t/(p*p);
00197 if (rr > 1.) rr = 0.;
00198 if (verboseLevel > 1)
00199 G4cout << "rr=" << rr << G4endl;
00200 G4double cost = 1. - rr;
00201 G4double sint = std::sqrt(std::max(rr*(2. - rr), 0.));
00202 if (sint == 0.) return &theParticleChange;
00203
00204 if (verboseLevel > 1)
00205 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
00206
00207
00208 G4double mass1 = aParticle->GetDefinition()->GetPDGMass();
00209 G4int Z=static_cast<G4int>(zTarget+.5);
00210 G4int A=static_cast<G4int>(atno2);
00211 if(G4UniformRand()<atno2-A) A++;
00212
00213 G4double mass2 =
00214 G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z)->GetPDGMass();
00215
00216
00217 G4double a = 1 + mass2/mass1;
00218 G4double b = -2.*p*cost;
00219 G4double c = p*p*(1 - mass2/mass1);
00220 G4double p1 = (-b+std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
00221 G4double px = p1*sint*std::sin(phi);
00222 G4double py = p1*sint*std::cos(phi);
00223 G4double pz = p1*cost;
00224
00225
00226 G4double etot = std::sqrt(mass1*mass1+p*p)+mass2;
00227 a = etot*etot-p*p*cost*cost;
00228 b = 2*p*p*(mass1*cost*cost-etot);
00229 c = p*p*p*p*sint*sint;
00230
00231 G4double de = (-b-std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
00232 G4double e1 = std::sqrt(p*p+mass1*mass1)-de;
00233 G4double p12=e1*e1-mass1*mass1;
00234 p1 = std::sqrt(std::max(1.*eV*eV,p12));
00235 px = p1*sint*std::sin(phi);
00236 py = p1*sint*std::cos(phi);
00237 pz = p1*cost;
00238
00239 if (verboseLevel > 1)
00240 {
00241 G4cout << "Relevant test "<<p<<" "<<p1<<" "<<cost<<" "<<de<<G4endl;
00242 G4cout << "p1/p = "<<p1/p<<" "<<mass1<<" "<<mass2<<" "<<a<<" "<<b<<" "<<c<<G4endl;
00243 G4cout << "rest = "<< b*b<<" "<<4.*a*c<<" "<<G4endl;
00244 G4cout << "make p1 = "<< p12<<" "<<e1*e1<<" "<<mass1*mass1<<" "<<G4endl;
00245 }
00246
00247 G4double pxinc = p*aParticle->Get4Momentum().vect().unit().x();
00248 G4double pyinc = p*aParticle->Get4Momentum().vect().unit().y();
00249 G4double pzinc = p*aParticle->Get4Momentum().vect().unit().z();
00250 if (verboseLevel > 1) {
00251 G4cout << "NOM SCAT " << px << " " << py << " " << pz << G4endl;
00252 G4cout << "INCIDENT " << pxinc << " " << pyinc << " " << pzinc << G4endl;
00253 }
00254
00255
00256 G4double pxnew, pynew, pznew;
00257 Defs1(p, px, py, pz, pxinc, pyinc, pzinc, &pxnew, &pynew, &pznew);
00258
00259 G4double pxre=pxinc-pxnew;
00260 G4double pyre=pyinc-pynew;
00261 G4double pzre=pzinc-pznew;
00262 G4ThreeVector it0(pxnew*GeV, pynew*GeV, pznew*GeV);
00263 if(p1>0)
00264 {
00265 pxnew = pxnew/p1;
00266 pynew = pynew/p1;
00267 pznew = pznew/p1;
00268 }
00269 else
00270 {
00271
00272
00273 pxnew = 0;
00274 pynew = 0;
00275 pznew = 1;
00276 }
00277 if (verboseLevel > 1) {
00278 G4cout << "DoIt: returning new momentum vector" << G4endl;
00279 G4cout << "DoIt: "<<pxinc << " " << pyinc << " " << pzinc <<" "<<p<< G4endl;
00280 G4cout << "DoIt: "<<pxnew << " " << pynew << " " << pznew <<" "<<p<< G4endl;
00281 }
00282
00283 if (aSecondary) {
00284 aSecondary->SetMomentumDirection(pxnew, pynew, pznew);
00285 } else {
00286 try
00287 {
00288 theParticleChange.SetMomentumChange(pxnew, pynew, pznew);
00289 theParticleChange.SetEnergyChange(std::sqrt(mass1*mass1+it0.mag2())-mass1);
00290 }
00291 catch(G4HadronicException)
00292 {
00293 std::cerr << "GHADException originating from components of G4LElastic"<<std::cout;
00294 throw;
00295 }
00296 G4ParticleDefinition * theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
00297 G4ThreeVector it(pxre*GeV, pyre*GeV, pzre*GeV);
00298 G4DynamicParticle * aSec =
00299 new G4DynamicParticle(theDef, it.unit(), it.mag2()/(2.*mass2));
00300 theParticleChange.AddSecondary(aSec);
00301
00302 }
00303 return &theParticleChange;
00304 }
00305
00306
00307
00308
00309
00310
00311
00312 G4int
00313 G4LElastic::Rtmi(G4double* x, G4double xli, G4double xri, G4double eps,
00314 G4int iend,
00315 G4double aa, G4double bb, G4double cc, G4double dd,
00316 G4double rr)
00317 {
00318 G4int ier = 0;
00319 G4double xl = xli;
00320 G4double xr = xri;
00321 *x = xl;
00322 G4double tol = *x;
00323 G4double f = Fctcos(tol, aa, bb, cc, dd, rr);
00324 if (f == 0.) return ier;
00325 G4double fl, fr;
00326 fl = f;
00327 *x = xr;
00328 tol = *x;
00329 f = Fctcos(tol, aa, bb, cc, dd, rr);
00330 if (f == 0.) return ier;
00331 fr = f;
00332
00333
00334 if (fl*fr >= 0.) {
00335 ier = 2;
00336 return ier;
00337 }
00338
00339
00340
00341 G4int i = 0;
00342 G4double tolf = 100.*eps;
00343
00344
00345 label4:
00346 i++;
00347
00348
00349 for (G4int k = 1; k <= iend; k++) {
00350 *x = 0.5*(xl + xr);
00351 tol = *x;
00352 f = Fctcos(tol, aa, bb, cc, dd, rr);
00353 if (f == 0.) return 0;
00354 if (f*fr < 0.) {
00355 tol = xl;
00356 xl = xr;
00357 xr = tol;
00358 tol = fl;
00359 fl = fr;
00360 fr = tol;
00361 }
00362 tol = f - fl;
00363 G4double a = f*tol;
00364 a = a + a;
00365 if (a < fr*(fr - fl) && i <= iend) goto label17;
00366 xr = *x;
00367 fr = f;
00368
00369
00370 tol = eps;
00371 a = std::abs(xr);
00372 if (a > 1.) tol = tol*a;
00373 if (std::abs(xr - xl) <= tol && std::abs(fr - fl) <= tolf) goto label14;
00374 }
00375
00376
00377
00378
00379
00380 ier = 1;
00381
00382 label14:
00383 if (std::abs(fr) > std::abs(fl)) {
00384 *x = xl;
00385 f = fl;
00386 }
00387 return ier;
00388
00389
00390 label17:
00391 G4double a = fr - f;
00392 G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
00393 G4double xm = *x;
00394 G4double fm = f;
00395 *x = xl - dx;
00396 tol = *x;
00397 f = Fctcos(tol, aa, bb, cc, dd, rr);
00398 if (f == 0.) return ier;
00399
00400
00401 tol = eps;
00402 a = std::abs(*x);
00403 if (a > 1) tol = tol*a;
00404 if (std::abs(dx) <= tol && std::abs(f) <= tolf) return ier;
00405
00406
00407 if (f*fl < 0.) {
00408 xr = *x;
00409 fr = f;
00410 }
00411 else {
00412 xl = *x;
00413 fl = f;
00414 xr = xm;
00415 fr = fm;
00416 }
00417 goto label4;
00418 }
00419
00420
00421
00422
00423 G4double
00424 G4LElastic::Fctcos(G4double t,
00425 G4double aa, G4double bb, G4double cc, G4double dd,
00426 G4double rr)
00427 {
00428 const G4double expxl = -82.;
00429 const G4double expxu = 82.;
00430
00431 G4double test1 = -bb*t;
00432 if (test1 > expxu) test1 = expxu;
00433 if (test1 < expxl) test1 = expxl;
00434
00435 G4double test2 = -dd*t;
00436 if (test2 > expxu) test2 = expxu;
00437 if (test2 < expxl) test2 = expxl;
00438
00439 return aa*std::exp(test1) + cc*std::exp(test2) - rr;
00440 }
00441
00442
00443 void
00444 G4LElastic::Defs1(G4double p, G4double px, G4double py, G4double pz,
00445 G4double pxinc, G4double pyinc, G4double pzinc,
00446 G4double* pxnew, G4double* pynew, G4double* pznew)
00447 {
00448
00449 G4double pt2 = pxinc*pxinc + pyinc*pyinc;
00450 if (pt2 > 0.) {
00451 G4double cost = pzinc/p;
00452 G4double sint1 = std::sqrt(std::abs((1. - cost )*(1.+cost)));
00453 G4double sint2 = std::sqrt(pt2)/p;
00454 G4double sint = 0.5*(sint1 + sint2);
00455 G4double ph = pi*0.5;
00456 if (pyinc < 0.) ph = pi*1.5;
00457 if (std::abs(pxinc) > 1.e-6) ph = std::atan2(pyinc, pxinc);
00458 G4double cosp = std::cos(ph);
00459 G4double sinp = std::sin(ph);
00460 if (verboseLevel > 1) {
00461 G4cout << "cost sint " << cost << " " << sint << G4endl;
00462 G4cout << "cosp sinp " << cosp << " " << sinp << G4endl;
00463 }
00464 *pxnew = cost*cosp*px - sinp*py + sint*cosp*pz;
00465 *pynew = cost*sinp*px + cosp*py + sint*sinp*pz;
00466 *pznew = -sint*px +cost*pz;
00467 }
00468 else {
00469 G4double cost=pzinc/p;
00470 *pxnew = cost*px;
00471 *pynew = py;
00472 *pznew = cost*pz;
00473 }
00474 }
00475
00476 const std::pair<G4double, G4double> G4LElastic::GetFatalEnergyCheckLevels() const
00477 {
00478 return std::pair<G4double, G4double>(5*perCent,250*GeV);
00479 }