Geant4-11
Public Member Functions | Private Member Functions | Private Attributes
G4INCL::CoulombNonRelativistic Class Reference

#include <G4INCLCoulombNonRelativistic.hh>

Inheritance diagram for G4INCL::CoulombNonRelativistic:
G4INCL::ICoulomb

Public Member Functions

IAvatarList bringToSurface (Cluster *const c, Nucleus *const n) const
 Modify the momentum of the incoming cluster and position it on the surface of the nucleus. More...
 
ParticleEntryAvatar * bringToSurface (Particle *const p, Nucleus *const n) const
 Modify the momentum of the particle and position it on the surface of the nucleus. More...
 
 CoulombNonRelativistic ()
 
void distortOut (ParticleList const &pL, Nucleus const *const n) const
 Modify the momenta of the outgoing particles. More...
 
G4double maxImpactParameter (ParticleSpecies const &p, const G4double kinE, Nucleus const *const n) const
 Return the maximum impact parameter for Coulomb-distorted trajectories. More...
 
virtual ~CoulombNonRelativistic ()
 

Private Member Functions

G4bool coulombDeviation (Particle *const p, Nucleus const *const n) const
 Perform Coulomb deviation. More...
 
G4double getCoulombRadius (ParticleSpecies const &p, Nucleus const *const n) const
 Get the Coulomb radius for a given particle. More...
 
G4double minimumDistance (Particle const *const p, Nucleus const *const n) const
 Return the minimum distance of approach in a head-on collision (b=0). More...
 
G4double minimumDistance (ParticleSpecies const &p, const G4double kineticEnergy, Nucleus const *const n) const
 Return the minimum distance of approach in a head-on collision (b=0). More...
 

Private Attributes

CoulombNone theCoulombNoneSlave
 Internal CoulombNone slave to generate the avatars. More...
 

Detailed Description

Definition at line 56 of file G4INCLCoulombNonRelativistic.hh.

Constructor & Destructor Documentation

◆ CoulombNonRelativistic()

G4INCL::CoulombNonRelativistic::CoulombNonRelativistic ( )
inline

Definition at line 58 of file G4INCLCoulombNonRelativistic.hh.

58{}

◆ ~CoulombNonRelativistic()

virtual G4INCL::CoulombNonRelativistic::~CoulombNonRelativistic ( )
inlinevirtual

Definition at line 59 of file G4INCLCoulombNonRelativistic.hh.

59{}

Member Function Documentation

◆ bringToSurface() [1/2]

IAvatarList G4INCL::CoulombNonRelativistic::bringToSurface ( Cluster *const  c,
Nucleus *const  n 
) const
virtual

Modify the momentum of the incoming cluster and position it on the surface of the nucleus.

This method performs non-relativistic distortion. The momenta of the particles that compose the cluster are also distorted.

Parameters
cincoming cluster
ndistorting nucleus

Implements G4INCL::ICoulomb.

Definition at line 63 of file G4INCLCoulombNonRelativistic.cc.

63 {
64 // Neutral clusters?!
65// assert(c->getZ()>0);
66
67 // Perform the actual Coulomb deviation
68 const G4bool success = coulombDeviation(c, n);
69 if(!success) {
70 return IAvatarList();
71 }
72
73 // Rely on the CoulombNone slave to compute the straight-line intersection
74 // and actually bring the particle to the surface of the nucleus
75 return theCoulombNoneSlave.bringToSurface(c,n);
76 }
G4bool
bool G4bool
Definition: G4Types.hh:86
G4INCL::CoulombNonRelativistic::theCoulombNoneSlave
CoulombNone theCoulombNoneSlave
Internal CoulombNone slave to generate the avatars.
Definition: G4INCLCoulombNonRelativistic.hh:160
G4INCL::CoulombNonRelativistic::coulombDeviation
G4bool coulombDeviation(Particle *const p, Nucleus const *const n) const
Perform Coulomb deviation.
Definition: G4INCLCoulombNonRelativistic.cc:137
G4INCL::CoulombNone::bringToSurface
ParticleEntryAvatar * bringToSurface(Particle *const p, Nucleus *const n) const
Position the particle on the surface of the nucleus.
Definition: G4INCLCoulombNone.cc:50
G4INCL::IAvatarList
UnorderedVector< IAvatar * > IAvatarList
Definition: G4INCLIAvatar.hh:107

References G4INCL::CoulombNone::bringToSurface(), coulombDeviation(), CLHEP::detail::n, and theCoulombNoneSlave.

◆ bringToSurface() [2/2]

ParticleEntryAvatar * G4INCL::CoulombNonRelativistic::bringToSurface ( Particle *const  p,
Nucleus *const  n 
) const
virtual

Modify the momentum of the particle and position it on the surface of the nucleus.

This method performs non-relativistic distortion.

Parameters
pincoming particle
ndistorting nucleus

Implements G4INCL::ICoulomb.

Definition at line 50 of file G4INCLCoulombNonRelativistic.cc.

50 {
51 // No distortion for neutral particles
52 if(p->getZ()!=0) {
53 const G4bool success = coulombDeviation(p, n);
54 if(!success) // transparent
55 return NULL;
56 }
57
58 // Rely on the CoulombNone slave to compute the straight-line intersection
59 // and actually bring the particle to the surface of the nucleus
60 return theCoulombNoneSlave.bringToSurface(p,n);
61 }

References G4INCL::CoulombNone::bringToSurface(), coulombDeviation(), G4INCL::Particle::getZ(), CLHEP::detail::n, and theCoulombNoneSlave.

◆ coulombDeviation()

G4bool G4INCL::CoulombNonRelativistic::coulombDeviation ( Particle *const  p,
Nucleus const *const  n 
) const
private

Perform Coulomb deviation.

Modifies the entrance angle of the particle and its impact parameter. Can be applied to Particles and Clusters.

The trajectory for an asymptotic impact parameter $b$ is parametrised as follows:

\[ r(\theta) = \frac{(1-e^2)r_0/2}{1-e \sin(\theta-\theta_R/2)}, \]

here $e$ is the hyperbola eccentricity:

\[ e = \sqrt{1+4b^2/r_0^2}; \]

$\theta_R$ is the Rutherford scattering angle:

\[ \theta_R = \pi - 2\arctan\left(\frac{2b}{r_0}\right) \]

$\theta$ ranges from $\pi$ (initial state) to $\theta_R$ (scattered particle) and $r_0$ is the minimum distance of approach in a head-on collision (see the minimumDistance() method).

Parameters
ppointer to the Particle
npointer to the Nucleus
Returns
false if below the barrier

Definition at line 137 of file G4INCLCoulombNonRelativistic.cc.

137 {
138 // Determine the rotation angle and the new impact parameter
139 ThreeVector positionTransverse = p->getTransversePosition();
140 const G4double impactParameterSquared = positionTransverse.mag2();
141 const G4double impactParameter = std::sqrt(impactParameterSquared);
142
143 // Some useful variables
144 const G4double theMinimumDistance = minimumDistance(p, n);
145 // deltaTheta2 = (pi - Rutherford scattering angle)/2
146 G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
147 if(deltaTheta2<0.)
148 deltaTheta2 += Math::pi;
149 const G4double eccentricity = 1./std::cos(deltaTheta2);
150
151 G4double newImpactParameter, alpha; // Parameters that must be determined by the deviation
152
153 const G4double radius = getCoulombRadius(p->getSpecies(), n);
154 const G4double impactParameterTangentSquared = radius*(radius-theMinimumDistance);
155 if(impactParameterSquared >= impactParameterTangentSquared) {
156 // The particle trajectory misses the Coulomb sphere
157 // In this case the new impact parameter is the minimum distance of
158 // approach of the hyperbola
159// assert(std::abs(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))>=eccentricity);
160 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity); // the minimum distance of approach
161 alpha = Math::piOverTwo - deltaTheta2; // half the Rutherford scattering angle
162 } else {
163 // The particle trajectory intersects the Coulomb sphere
164
165 // Compute the entrance angle
166 const G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
167 / eccentricity;
168 const G4double thetaIn = Math::twoPi - Math::arcCos(argument) - deltaTheta2;
169
170 // Velocity angle at the entrance point
171 alpha = std::atan((1+std::cos(thetaIn))
172 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)))
173 * Math::sign(theMinimumDistance);
174 // New impact parameter
175 newImpactParameter = radius * std::sin(thetaIn - alpha);
176 }
177
178 // Modify the impact parameter of the particle
179 positionTransverse *= newImpactParameter/positionTransverse.mag();
180 const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse;
181 p->setPosition(theNewPosition);
182
183 // Determine the rotation axis for the incoming particle
184 const ThreeVector &momentum = p->getMomentum();
185 ThreeVector rotationAxis = momentum.vector(positionTransverse);
186 const G4double axisLength = rotationAxis.mag();
187 // Apply the rotation
188 if(axisLength>1E-20) {
189 rotationAxis /= axisLength;
190 p->rotatePositionAndMomentum(alpha, rotationAxis);
191 }
192
193 return true;
194 }
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:80
G4double
double G4double
Definition: G4Types.hh:83
G4INCL::CoulombNonRelativistic::getCoulombRadius
G4double getCoulombRadius(ParticleSpecies const &p, Nucleus const *const n) const
Get the Coulomb radius for a given particle.
Definition: G4INCLCoulombNonRelativistic.cc:196
G4INCL::CoulombNonRelativistic::minimumDistance
G4double minimumDistance(ParticleSpecies const &p, const G4double kineticEnergy, Nucleus const *const n) const
Return the minimum distance of approach in a head-on collision (b=0).
Definition: G4INCLCoulombNonRelativistic.hh:98
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:68
G4INCL::Math::twoPi
const G4double twoPi
Definition: G4INCLGlobals.hh:69
G4INCL::Math::arcCos
G4double arcCos(const G4double x)
Calculates arccos with some tolerance on illegal arguments.
Definition: G4INCLGlobals.cc:103
G4INCL::Math::piOverTwo
const G4double piOverTwo
Definition: G4INCLGlobals.hh:71
G4INCL::Math::sign
G4int sign(const T t)
Definition: G4INCLGlobals.hh:107

References alpha, G4INCL::Math::arcCos(), getCoulombRadius(), G4INCL::Particle::getLongitudinalPosition(), G4INCL::Particle::getMomentum(), G4INCL::Particle::getSpecies(), G4INCL::Particle::getTransversePosition(), G4INCL::ThreeVector::mag(), G4INCL::ThreeVector::mag2(), minimumDistance(), CLHEP::detail::n, G4INCL::Math::pi, G4INCL::Math::piOverTwo, G4INCL::Particle::rotatePositionAndMomentum(), G4INCL::Particle::setPosition(), G4INCL::Math::sign(), G4INCL::Math::twoPi, and G4INCL::ThreeVector::vector().

Referenced by bringToSurface().

◆ distortOut()

void G4INCL::CoulombNonRelativistic::distortOut ( ParticleList const &  pL,
Nucleus const *const  n 
) const
virtual

Modify the momenta of the outgoing particles.

This method performs non-relativistic distortion.

Parameters
pLlist of outgoing particles
ndistorting nucleus

Implements G4INCL::ICoulomb.

Definition at line 78 of file G4INCLCoulombNonRelativistic.cc.

79 {
80
81 for(ParticleIter particle=pL.begin(), e=pL.end(); particle!=e; ++particle) {
82
83 const G4int Z = (*particle)->getZ();
84 if(Z == 0) continue;
85
86 const G4double tcos=1.-0.000001;
87
88 const G4double et1 = PhysicalConstants::eSquared * nucleus->getZ();
89 const G4double transmissionRadius =
90 nucleus->getDensity()->getTransmissionRadius(*particle);
91
92 const ThreeVector position = (*particle)->getPosition();
93 ThreeVector momentum = (*particle)->getMomentum();
94 const G4double r = position.mag();
95 const G4double p = momentum.mag();
96 const G4double cosTheta = position.dot(momentum)/(r*p);
97 if(cosTheta < 0.999999) {
98 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
99 const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
100 if(eta > transmissionRadius-0.0001) {
101 // If below the Coulomb barrier, radial emission:
102 momentum = position * (p/r);
103 (*particle)->setMomentum(momentum);
104 } else {
105 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
106 4. * std::pow(transmissionRadius*sinTheta,2)
107 * (1.-eta/transmissionRadius)));
108 const G4double bInf = std::sqrt(b0*(b0-eta));
109 const G4double thr = std::atan(eta/(2.*bInf));
110 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
111 b0/transmissionRadius;
112 if(uTemp>tcos) uTemp=tcos;
113 const G4double thd = Math::arcCos(cosTheta)-Math::piOverTwo + thr +
114 Math::arcCos(uTemp);
115 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
116 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
117 const ThreeVector newMomentum = momentum*c1 + position*c2;
118 (*particle)->setMomentum(newMomentum);
119 }
120 }
121 }
122 }
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
G4INCL::PhysicalConstants::eSquared
const G4double eSquared
Coulomb conversion factor [MeV*fm].
Definition: G4INCLGlobals.hh:64
G4INCL::ParticleIter
ParticleList::const_iterator ParticleIter
Definition: G4INCLParticle.hh:72

References G4INCL::Math::arcCos(), G4INCL::PhysicalConstants::eSquared, G4INCL::Nucleus::getDensity(), G4INCL::NuclearDensity::getTransmissionRadius(), G4INCL::Particle::getZ(), G4INCL::ThreeVector::mag(), G4INCL::Math::piOverTwo, and Z.

◆ getCoulombRadius()

G4double G4INCL::CoulombNonRelativistic::getCoulombRadius ( ParticleSpecies const &  p,
Nucleus const *const  n 
) const
private

Get the Coulomb radius for a given particle.

That's the radius of the sphere that the Coulomb trajectory of the incoming particle should intersect. The intersection point is used to determine the effective impact parameter of the trajectory and the new entrance angle.

If the particle is not a Cluster, the Coulomb radius reduces to the surface radius. We use a parametrisation for d, t, He3 and alphas. For heavier clusters we fall back to the surface radius.

Parameters
pthe particle species
nthe deflecting nucleus
Returns
Coulomb radius

Definition at line 196 of file G4INCLCoulombNonRelativistic.cc.

196 {
197 if(p.theType == Composite) {
198 const G4int Zp = p.theZ;
199 const G4int Ap = p.theA;
200 const G4int Zt = n->getZ();
201 const G4int At = n->getA();
202 G4double barr, radius = 0.;
203 if(Zp==1 && Ap==2) { // d
204 barr = 0.2565*Math::pow23((G4double)At)-0.78;
205 radius = PhysicalConstants::eSquared*Zp*Zt/barr - 2.5;
206 } else if(Zp==1 && Ap==3) { // t
207 barr = 0.5*(0.5009*Math::pow23((G4double)At)-1.16);
208 radius = PhysicalConstants::eSquared*Zt/barr - 0.5;
209 } else if(Zp==2) { // alpha, He3
210 barr = 0.5939*Math::pow23((G4double)At)-1.64;
211 radius = PhysicalConstants::eSquared*Zp*Zt/barr - 0.5;
212 } else if(Zp>2) {
213 // Coulomb radius from the Shen model
214 const G4double Ap13 = Math::pow13((G4double)Ap);
215 const G4double At13 = Math::pow13((G4double)At);
216 const G4double rp = 1.12*Ap13 - 0.94/Ap13;
217 const G4double rt = 1.12*At13 - 0.94/At13;
218 const G4double someRadius = rp+rt+3.2;
219 const G4double theShenBarrier = PhysicalConstants::eSquared*Zp*Zt/someRadius - rt*rp/(rt+rp);
220 radius = PhysicalConstants::eSquared*Zp*Zt/theShenBarrier;
221 }
222 if(radius<=0.) {
223 radius = ParticleTable::getLargestNuclearRadius(Ap,Zp) + ParticleTable::getLargestNuclearRadius(At, Zt);
224 INCL_ERROR("Negative Coulomb radius! Using the sum of nuclear radii = " << radius << '\n');
225 }
226 INCL_DEBUG("Coulomb radius for particle "
227 << ParticleTable::getShortName(p) << " in nucleus A=" << At <<
228 ", Z=" << Zt << ": " << radius << '\n');
229 return radius;
230 } else
231 return n->getUniverseRadius();
232 }
INCL_ERROR
#define INCL_ERROR(x)
Definition: G4INCLLogger.hh:225
INCL_DEBUG
#define INCL_DEBUG(x)
Definition: G4INCLLogger.hh:240
G4INCL::Math::pow13
G4double pow13(G4double x)
Definition: G4INCLGlobals.hh:88
G4INCL::Math::pow23
G4double pow23(G4double x)
Definition: G4INCLGlobals.hh:96
G4INCL::ParticleTable::getLargestNuclearRadius
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
Definition: G4INCLParticleTable.cc:993
G4INCL::ParticleTable::getShortName
std::string getShortName(const ParticleType t)
Get the short INCL name of the particle.
Definition: G4INCLParticleTable.cc:627

References G4INCL::Composite, G4INCL::PhysicalConstants::eSquared, G4INCL::ParticleTable::getLargestNuclearRadius(), G4INCL::ParticleTable::getShortName(), INCL_DEBUG, INCL_ERROR, CLHEP::detail::n, G4INCL::Math::pow13(), G4INCL::Math::pow23(), G4INCL::ParticleSpecies::theA, G4INCL::ParticleSpecies::theType, and G4INCL::ParticleSpecies::theZ.

Referenced by coulombDeviation().

◆ maxImpactParameter()

G4double G4INCL::CoulombNonRelativistic::maxImpactParameter ( ParticleSpecies const &  p,
const G4double  kinE,
Nucleus const *const  n 
) const
virtual

Return the maximum impact parameter for Coulomb-distorted trajectories.

Implements G4INCL::ICoulomb.

Definition at line 124 of file G4INCLCoulombNonRelativistic.cc.

125 {
126 const G4double theMinimumDistance = minimumDistance(p, kinE, n);
127 G4double rMax = n->getUniverseRadius();
128 if(p.theType == Composite)
129 rMax += 2.*ParticleTable::getLargestNuclearRadius(p.theA, p.theZ);
130 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
131 if(theMaxImpactParameterSquared<=0.)
132 return 0.;
133 const G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
134 return theMaxImpactParameter;
135 }

References G4INCL::Composite, G4INCL::ParticleTable::getLargestNuclearRadius(), minimumDistance(), CLHEP::detail::n, G4INCL::ParticleSpecies::theA, G4INCL::ParticleSpecies::theType, and G4INCL::ParticleSpecies::theZ.

◆ minimumDistance() [1/2]

G4double G4INCL::CoulombNonRelativistic::minimumDistance ( Particle const *const  p,
Nucleus const *const  n 
) const
inlineprivate

Return the minimum distance of approach in a head-on collision (b=0).

Definition at line 110 of file G4INCLCoulombNonRelativistic.hh.

110 {
111 return minimumDistance(p->getSpecies(), p->getKineticEnergy(), n);
112 }

References G4INCL::Particle::getKineticEnergy(), G4INCL::Particle::getSpecies(), minimumDistance(), and CLHEP::detail::n.

◆ minimumDistance() [2/2]

G4double G4INCL::CoulombNonRelativistic::minimumDistance ( ParticleSpecies const &  p,
const G4double  kineticEnergy,
Nucleus const *const  n 
) const
inlineprivate

Return the minimum distance of approach in a head-on collision (b=0).

Definition at line 98 of file G4INCLCoulombNonRelativistic.hh.

98 {
99 const G4double particleMass = ParticleTable::getTableSpeciesMass(p);
100 const G4double nucleusMass = n->getTableMass();
101 const G4double reducedMass = particleMass*nucleusMass/(particleMass+nucleusMass);
102 const G4double kineticEnergyInCM = kineticEnergy * reducedMass / particleMass;
103 const G4double theMinimumDistance = PhysicalConstants::eSquared * p.theZ * n->getZ() * particleMass
104 / (kineticEnergyInCM * reducedMass);
105 INCL_DEBUG("Minimum distance of approach due to Coulomb = " << theMinimumDistance << '\n');
106 return theMinimumDistance;
107 }
G4INCL::ParticleTable::getTableSpeciesMass
G4double getTableSpeciesMass(const ParticleSpecies &p)
Definition: G4INCLParticleTable.cc:843

References G4INCL::PhysicalConstants::eSquared, G4INCL::ParticleTable::getTableSpeciesMass(), INCL_DEBUG, CLHEP::detail::n, and G4INCL::ParticleSpecies::theZ.

Referenced by coulombDeviation(), maxImpactParameter(), and minimumDistance().

Field Documentation

◆ theCoulombNoneSlave

CoulombNone G4INCL::CoulombNonRelativistic::theCoulombNoneSlave
private

Internal CoulombNone slave to generate the avatars.

Definition at line 160 of file G4INCLCoulombNonRelativistic.hh.

Referenced by bringToSurface().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3