33 #define INCLXX_IN_GEANT4_MODE 1
52 const G4bool success = coulombDeviation(p, n);
67 const G4bool success = coulombDeviation(c, n);
78 Nucleus const *
const nucleus)
const {
80 for(
ParticleIter particle=pL.begin(), e=pL.end(); particle!=e; ++particle) {
82 const G4int Z = (*particle)->getZ();
96 if(cosTheta < 0.999999) {
97 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
98 const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
99 if(eta > transmissionRadius-0.0001) {
101 momentum = position * (p/r);
102 (*particle)->setMomentum(momentum);
104 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
105 4. * std::pow(transmissionRadius*sinTheta,2)
106 * (1.-eta/transmissionRadius)));
107 const G4double bInf = std::sqrt(b0*(b0-eta));
108 const G4double thr = std::atan(eta/(2.*bInf));
109 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
110 b0/transmissionRadius;
111 if(uTemp>tcos) uTemp=tcos;
114 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
115 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
116 const ThreeVector newMomentum = momentum*c1 + position*c2;
117 (*particle)->setMomentum(newMomentum);
125 const G4double theMinimumDistance = minimumDistance(p, kinE, n);
129 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
130 if(theMaxImpactParameterSquared<=0.)
132 const G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
133 return theMaxImpactParameter;
139 const G4double impactParameterSquared = positionTransverse.
mag2();
140 const G4double impactParameter = std::sqrt(impactParameterSquared);
143 const G4double theMinimumDistance = minimumDistance(p, n);
145 const G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
146 const G4double eccentricity = 1./std::cos(deltaTheta2);
151 const G4double impactParameterTangentSquared = radius*(radius-theMinimumDistance);
152 if(impactParameterSquared >= impactParameterTangentSquared) {
157 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity);
163 const G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
169 alpha = std::atan((1+std::cos(thetaIn))
170 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)));
172 newImpactParameter = radius * std::sin(thetaIn - alpha);
176 positionTransverse *= newImpactParameter/positionTransverse.
mag();
182 ThreeVector rotationAxis = momentum.
vector(positionTransverse);
185 if(axisLength>1E-20) {
186 rotationAxis /= axisLength;
187 p->
rotate(alpha, rotationAxis);
193 G4double CoulombNonRelativistic::getCoulombRadius(ParticleSpecies
const &p, Nucleus
const *
const n)
const {
195 const G4int Zp = p.theZ;
196 const G4int Ap = p.theA;
197 const G4int Zt = n->getZ();
198 const G4int At = n->getA();
203 }
else if(Zp==1 && Ap==3) {
213 const G4double rp = 1.12*Ap13 - 0.94/Ap13;
214 const G4double rt = 1.12*At13 - 0.94/At13;
215 const G4double someRadius = rp+rt+3.2;
221 INCL_ERROR(
"Negative Coulomb radius! Using the sum of nuclear radii = " << radius << std::endl);
225 ", Z=" << Zt <<
": " << radius << std::endl);
228 return n->getUniverseRadius();
ParticleEntryAvatar * bringToSurface(Particle *const p, Nucleus *const n) const
Modify the momentum of the particle and position it on the surface of the nucleus.
Class for non-relativistic Coulomb distortion.
const G4double eSquared
Coulomb conversion factor [MeV*fm].
G4double dot(const ThreeVector &v) const
const G4INCL::ThreeVector & getMomentum() const
ThreeVector vector(const ThreeVector &v) const
G4double pow23(G4double x)
virtual void rotate(const G4double angle, const ThreeVector &axis)
Rotate the particle position and momentum.
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n) const
Return the maximum impact parameter for Coulomb-distorted trajectories.
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
virtual void setPosition(const G4INCL::ThreeVector &position)
G4int getZ() const
Returns the charge number.
UnorderedVector< IAvatar * > IAvatarList
ParticleEntryAvatar * bringToSurface(Particle *const p, Nucleus *const n) const
Position the particle on the surface of the nucleus.
ThreeVector getTransversePosition() const
Transverse component of the position w.r.t. the momentum.
NuclearDensity const * getDensity() const
Getter for theDensity.
virtual G4INCL::ParticleSpecies getSpecies() const
Get the particle species.
G4double getUniverseRadius() const
Getter for theUniverseRadius.
G4double getTransmissionRadius(Particle const *const p) const
The radius used for calculating the transmission coefficient.
G4double pow13(G4double x)
std::string getShortName(const ParticleType t)
Get the short INCL name of the particle.
ThreeVector getLongitudinalPosition() const
Longitudinal component of the position w.r.t. the momentum.
ParticleList::const_iterator ParticleIter
void distortOut(ParticleList const &pL, Nucleus const *const n) const
Modify the momenta of the outgoing particles.