57 outFile <<
"G4LElastic is one of the Low Energy Parameterized (LEP)\n"
58 <<
"models used to implement elastic hadron scattering from nuclei.\n"
59 <<
"It is a re-engineered version of the GHEISHA code of\n"
60 <<
"H. Fesefeldt. It performs simplified two-body elastic\n"
61 <<
"scattering for all long-lived hadronic projectiles by using\n"
62 <<
"a two-exponential parameterization in momentum transfer.\n"
63 <<
"It is valid for incident hadrons of all energies.\n";
141 G4cout <<
"G4LElastic::DoIt: Incident particle p=" << p <<
" GeV" <<
G4endl;
153 aa = std::pow(atno2, 1.63);
154 bb = 14.5*std::pow(atno2, 0.66);
155 cc = 1.4*std::pow(atno2, 0.33);
159 aa = std::pow(atno2, 1.33);
160 bb = 60.*std::pow(atno2, 0.33);
161 cc = 0.4*std::pow(atno2, 0.40);
169 G4cout << aa <<
" " << bb <<
" " << cc <<
" " << dd <<
" " << rr <<
G4endl;
174 G4cout <<
"t1,Fctcos " << t1 <<
" " << Fctcos(t1, aa, bb, cc, dd, rr) <<
176 G4cout <<
"t2,Fctcos " << t2 <<
" " << Fctcos(t2, aa, bb, cc, dd, rr) <<
183 ier1 = Rtmi(&t, t1, t2, eps, ind1,
187 G4cout <<
"t, Fctcos " << t <<
" " << Fctcos(t, aa, bb, cc, dd, rr) <<
190 if (ier1 != 0) t = 0.25*(3.*t1 + t2);
192 G4cout <<
"t, Fctcos " << t <<
" " << Fctcos(t, aa, bb, cc, dd, rr) <<
197 if (rr > 1.) rr = 0.;
201 G4double sint = std::sqrt(std::max(rr*(2. - rr), 0.));
205 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
220 G4double p1 = (-b+std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
221 G4double px = p1*sint*std::sin(phi);
222 G4double py = p1*sint*std::cos(phi);
226 G4double etot = std::sqrt(mass1*mass1+p*p)+mass2;
227 a = etot*etot-p*p*cost*cost;
228 b = 2*p*p*(mass1*cost*cost-etot);
229 c = p*p*p*p*sint*sint;
231 G4double de = (-b-std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
232 G4double e1 = std::sqrt(p*p+mass1*mass1)-de;
234 p1 = std::sqrt(std::max(1.*eV*eV,p12));
235 px = p1*sint*std::sin(phi);
236 py = p1*sint*std::cos(phi);
241 G4cout <<
"Relevant test "<<p<<
" "<<p1<<
" "<<cost<<
" "<<de<<
G4endl;
242 G4cout <<
"p1/p = "<<p1/p<<
" "<<mass1<<
" "<<mass2<<
" "<<a<<
" "<<b<<
" "<<c<<
G4endl;
244 G4cout <<
"make p1 = "<< p12<<
" "<<e1*e1<<
" "<<mass1*mass1<<
" "<<
G4endl;
251 G4cout <<
"NOM SCAT " << px <<
" " << py <<
" " << pz <<
G4endl;
252 G4cout <<
"INCIDENT " << pxinc <<
" " << pyinc <<
" " << pzinc <<
G4endl;
257 Defs1(p, px, py, pz, pxinc, pyinc, pzinc, &pxnew, &pynew, &pznew);
278 G4cout <<
"DoIt: returning new momentum vector" <<
G4endl;
279 G4cout <<
"DoIt: "<<pxinc <<
" " << pyinc <<
" " << pzinc <<
" "<<p<<
G4endl;
280 G4cout <<
"DoIt: "<<pxnew <<
" " << pynew <<
" " << pznew <<
" "<<p<<
G4endl;
293 std::cerr <<
"GHADException originating from components of G4LElastic"<<std::cout;
323 G4double f = Fctcos(tol, aa, bb, cc, dd, rr);
324 if (f == 0.)
return ier;
329 f = Fctcos(tol, aa, bb, cc, dd, rr);
330 if (f == 0.)
return ier;
349 for (
G4int k = 1; k <= iend; k++) {
352 f = Fctcos(tol, aa, bb, cc, dd, rr);
353 if (f == 0.)
return 0;
365 if (a < fr*(fr - fl) && i <= iend)
goto label17;
372 if (a > 1.) tol = tol*a;
373 if (std::abs(xr - xl) <= tol && std::abs(fr - fl) <= tolf)
goto label14;
383 if (std::abs(fr) > std::abs(fl)) {
392 G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
397 f = Fctcos(tol, aa, bb, cc, dd, rr);
398 if (f == 0.)
return ier;
403 if (a > 1) tol = tol*a;
404 if (std::abs(dx) <= tol && std::abs(f) <= tolf)
return ier;
432 if (test1 > expxu) test1 = expxu;
433 if (test1 < expxl) test1 = expxl;
436 if (test2 > expxu) test2 = expxu;
437 if (test2 < expxl) test2 = expxl;
439 return aa*std::exp(test1) + cc*std::exp(test2) - rr;
449 G4double pt2 = pxinc*pxinc + pyinc*pyinc;
452 G4double sint1 = std::sqrt(std::abs((1. - cost )*(1.+cost)));
454 G4double sint = 0.5*(sint1 + sint2);
456 if (pyinc < 0.) ph =
pi*1.5;
457 if (std::abs(pxinc) > 1.e-6) ph = std::atan2(pyinc, pxinc);
461 G4cout <<
"cost sint " << cost <<
" " << sint <<
G4endl;
462 G4cout <<
"cosp sinp " << cosp <<
" " << sinp <<
G4endl;
464 *pxnew = cost*cosp*px - sinp*py + sint*cosp*pz;
465 *pynew = cost*sinp*px + cosp*py + sint*sinp*pz;
466 *pznew = -sint*px +cost*pz;
478 return std::pair<G4double, G4double>(5*perCent,250*GeV);
G4DLLIMPORT std::ostream G4cout
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiOmegaMinus * AntiOmegaMinus()
static G4AntiProton * AntiProton()
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiXiMinus * AntiXiMinus()
static G4AntiXiZero * AntiXiZero()
void SetMomentumDirection(const G4ThreeVector &aDirection)
void SetMomentum(const G4ThreeVector &momentum)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP)
void SetEnergyChange(G4double anEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
G4double GetTotalMomentum() const
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4HadFinalState theParticleChange
void SetMinEnergy(G4double anEnergy)
void SetMaxEnergy(const G4double anEnergy)
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
G4LElastic(const G4String &name="G4LElastic")
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
virtual void ModelDescription(std::ostream &outFile) const
static G4Lambda * Lambda()
static G4Neutron * Neutron()
static G4OmegaMinus * OmegaMinus()
G4double GetPDGMass() const
G4ParticleDefinition * FindIon(G4int atomicNumber, G4int atomicMass, G4double excitationEnergy)
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4Proton * Proton()
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4XiMinus * XiMinus()
static G4XiZero * XiZero()