00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00044 #include "G4INCLCoulombNonRelativistic.hh"
00045 #include "G4INCLGlobals.hh"
00046
00047 namespace G4INCL {
00048
00049 ParticleEntryAvatar *CoulombNonRelativistic::bringToSurface(Particle * const p, Nucleus * const n) const {
00050
00051 if(p->getZ()!=0) {
00052 const G4bool success = coulombDeviation(p, n);
00053 if(!success)
00054 return NULL;
00055 }
00056
00057
00058
00059 return theCoulombNoneSlave.bringToSurface(p,n);
00060 }
00061
00062 IAvatarList CoulombNonRelativistic::bringToSurface(Cluster * const c, Nucleus * const n) const {
00063
00064
00065
00066
00067 const G4bool success = coulombDeviation(c, n);
00068 if(!success) {
00069 return IAvatarList();
00070 }
00071
00072
00073
00074 return theCoulombNoneSlave.bringToSurface(c,n);
00075 }
00076
00077 void CoulombNonRelativistic::distortOut(ParticleList const &pL,
00078 Nucleus const * const nucleus) const {
00079
00080 for(ParticleIter particle=pL.begin(); particle!=pL.end(); ++particle) {
00081
00082 const G4int Z = (*particle)->getZ();
00083 if(Z == 0) continue;
00084
00085 const G4double tcos=1.-0.000001;
00086
00087 const G4double et1 = PhysicalConstants::eSquared * nucleus->getZ();
00088 const G4double transmissionRadius =
00089 nucleus->getDensity()->getTransmissionRadius(*particle);
00090
00091 const ThreeVector position = (*particle)->getPosition();
00092 ThreeVector momentum = (*particle)->getMomentum();
00093 const G4double r = position.mag();
00094 const G4double p = momentum.mag();
00095 const G4double cosTheta = position.dot(momentum)/(r*p);
00096 if(cosTheta < 0.999999) {
00097 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
00098 const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
00099 if(eta > transmissionRadius-0.0001) {
00100
00101 momentum = position * (p/r);
00102 (*particle)->setMomentum(momentum);
00103 } else {
00104 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
00105 4. * std::pow(transmissionRadius*sinTheta,2)
00106 * (1.-eta/transmissionRadius)));
00107 const G4double bInf = std::sqrt(b0*(b0-eta));
00108 const G4double thr = std::atan(eta/(2.*bInf));
00109 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
00110 b0/transmissionRadius;
00111 if(uTemp>tcos) uTemp=tcos;
00112 const G4double thd = std::acos(cosTheta)-Math::piOverTwo + thr +
00113 std::acos(uTemp);
00114 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
00115 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
00116 const ThreeVector newMomentum = momentum*c1 + position*c2;
00117 (*particle)->setMomentum(newMomentum);
00118 }
00119 }
00120 }
00121 }
00122
00123 G4double CoulombNonRelativistic::maxImpactParameter(ParticleSpecies const &p, const G4double kinE,
00124 Nucleus const * const n) const {
00125 G4double theMaxImpactParameter = maxImpactParameterParticle(p, kinE, n);
00126 if(theMaxImpactParameter <= 0.)
00127 return 0.;
00128 if(p.theType == Composite)
00129 theMaxImpactParameter += 2.*ParticleTable::getNuclearRadius(p.theA, p.theZ);
00130 return theMaxImpactParameter;
00131 }
00132
00133 G4double CoulombNonRelativistic::maxImpactParameterParticle(ParticleSpecies const &p, const G4double kinE,
00134 Nucleus const * const n) const {
00135 const G4double theMinimumDistance = minimumDistance(p, kinE, n);
00136 const G4double rMax = n->getCoulombRadius(p);
00137 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
00138 if(theMaxImpactParameterSquared<=0.)
00139 return 0.;
00140 G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
00141 return theMaxImpactParameter;
00142 }
00143
00144 G4bool CoulombNonRelativistic::coulombDeviation(Particle * const p, Nucleus const * const n) const {
00145
00146 ThreeVector positionTransverse = p->getTransversePosition();
00147 const G4double impactParameter = positionTransverse.mag();
00148
00149
00150 const G4double theMinimumDistance = minimumDistance(p, n);
00151
00152 const G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
00153 const G4double eccentricity = 1./std::cos(deltaTheta2);
00154
00155 G4double newImpactParameter, alpha;
00156
00157 ParticleSpecies aSpecies = p->getSpecies();
00158 G4double kineticEnergy = p->getKineticEnergy();
00159
00160
00161
00162 if(impactParameter>maxImpactParameterParticle(aSpecies, kineticEnergy, n)) {
00163
00164
00165
00166
00167
00168 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity);
00169 alpha = Math::piOverTwo - deltaTheta2;
00170 } else {
00171
00172
00173
00174 const G4double radius = n->getCoulombRadius(p->getSpecies());
00175 G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
00176 / eccentricity;
00177 const G4double thetaIn = Math::twoPi - std::acos(argument) - deltaTheta2;
00178
00179
00180 alpha = std::atan((1+std::cos(thetaIn))
00181 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)));
00182
00183 newImpactParameter = radius * std::sin(thetaIn - alpha);
00184 }
00185
00186
00187 positionTransverse *= newImpactParameter/positionTransverse.mag();
00188 const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse;
00189 p->setPosition(theNewPosition);
00190
00191
00192 const ThreeVector &momentum = p->getMomentum();
00193 ThreeVector rotationAxis = momentum.vector(positionTransverse);
00194 const G4double axisLength = rotationAxis.mag();
00195
00196 if(axisLength>1E-20) {
00197 rotationAxis /= axisLength;
00198 p->rotate(alpha, rotationAxis);
00199 }
00200
00201 return true;
00202 }
00203
00204 }