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(); particle!=pL.end(); ++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 G4double theMaxImpactParameter = maxImpactParameterParticle(p, kinE, n);
126 if(theMaxImpactParameter <= 0.)
130 return theMaxImpactParameter;
135 const G4double theMinimumDistance = minimumDistance(p, kinE, n);
137 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
138 if(theMaxImpactParameterSquared<=0.)
140 G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
141 return theMaxImpactParameter;
144 G4bool CoulombNonRelativistic::coulombDeviation(Particle *
const p, Nucleus
const *
const n)
const {
146 ThreeVector positionTransverse = p->getTransversePosition();
147 const G4double impactParameter = positionTransverse.mag();
150 const G4double theMinimumDistance = minimumDistance(p, n);
152 const G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
153 const G4double eccentricity = 1./std::cos(deltaTheta2);
157 ParticleSpecies aSpecies = p->getSpecies();
158 G4double kineticEnergy = p->getKineticEnergy();
162 if(impactParameter>maxImpactParameterParticle(aSpecies, kineticEnergy, n)) {
168 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity);
174 const G4double radius = n->getCoulombRadius(p->getSpecies());
175 G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
180 alpha = std::atan((1+std::cos(thetaIn))
181 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)));
183 newImpactParameter = radius * std::sin(thetaIn - alpha);
187 positionTransverse *= newImpactParameter/positionTransverse.mag();
188 const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse;
189 p->setPosition(theNewPosition);
192 const ThreeVector &momentum = p->getMomentum();
193 ThreeVector rotationAxis = momentum.vector(positionTransverse);
194 const G4double axisLength = rotationAxis.mag();
196 if(axisLength>1E-20) {
197 rotationAxis /= axisLength;
198 p->rotate(alpha, rotationAxis);