34 #define INCLXX_IN_GEANT4_MODE 1
79 Nucleus const *
const nucleus)
const {
83 const G4int Z = (*particle)->getZ();
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) {
102 momentum = position * (p/
r);
103 (*particle)->setMomentum(momentum);
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;
115 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
116 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
118 (*particle)->setMomentum(newMomentum);
130 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
131 if(theMaxImpactParameterSquared<=0.)
133 const G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
134 return theMaxImpactParameter;
140 const G4double impactParameterSquared = positionTransverse.
mag2();
141 const G4double impactParameter = std::sqrt(impactParameterSquared);
146 G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
149 const G4double eccentricity = 1./std::cos(deltaTheta2);
154 const G4double impactParameterTangentSquared = radius*(radius-theMinimumDistance);
155 if(impactParameterSquared >= impactParameterTangentSquared) {
160 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity);
166 const G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
171 alpha = std::atan((1+std::cos(thetaIn))
172 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)))
175 newImpactParameter = radius * std::sin(thetaIn - alpha);
179 positionTransverse *= newImpactParameter/positionTransverse.
mag();
188 if(axisLength>1
E-20) {
189 rotationAxis /= axisLength;
206 }
else if(Zp==1 && Ap==3) {
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;
224 INCL_ERROR(
"Negative Coulomb radius! Using the sum of nuclear radii = " << radius <<
'\n');
228 ", Z=" << Zt <<
": " << radius <<
'\n');