Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics 28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA 29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA 31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege 32 // Davide Mancusi, CEA-Saclay, France << 33 // 32 // 34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1 35 34 36 #include "globals.hh" 35 #include "globals.hh" 37 36 38 /** \file G4INCLCoulombNonRelativistic.cc 37 /** \file G4INCLCoulombNonRelativistic.cc 39 * \brief Class for non-relativistic Coulomb d 38 * \brief Class for non-relativistic Coulomb distortion. 40 * 39 * 41 * \date 14 February 2011 40 * \date 14 February 2011 42 * \author Davide Mancusi 41 * \author Davide Mancusi 43 */ 42 */ 44 43 45 #include "G4INCLCoulombNonRelativistic.hh" 44 #include "G4INCLCoulombNonRelativistic.hh" 46 #include "G4INCLGlobals.hh" 45 #include "G4INCLGlobals.hh" 47 46 48 namespace G4INCL { 47 namespace G4INCL { 49 48 50 ParticleEntryAvatar *CoulombNonRelativistic: 49 ParticleEntryAvatar *CoulombNonRelativistic::bringToSurface(Particle * const p, Nucleus * const n) const { 51 // No distortion for neutral particles 50 // No distortion for neutral particles 52 if(p->getZ()!=0) { 51 if(p->getZ()!=0) { 53 const G4bool success = coulombDeviation( 52 const G4bool success = coulombDeviation(p, n); 54 if(!success) // transparent 53 if(!success) // transparent 55 return NULL; 54 return NULL; 56 } 55 } 57 56 58 // Rely on the CoulombNone slave to comput 57 // Rely on the CoulombNone slave to compute the straight-line intersection 59 // and actually bring the particle to the 58 // and actually bring the particle to the surface of the nucleus 60 return theCoulombNoneSlave.bringToSurface( 59 return theCoulombNoneSlave.bringToSurface(p,n); 61 } 60 } 62 61 63 IAvatarList CoulombNonRelativistic::bringToS 62 IAvatarList CoulombNonRelativistic::bringToSurface(Cluster * const c, Nucleus * const n) const { 64 // Neutral clusters?! 63 // Neutral clusters?! 65 // assert(c->getZ()>0); 64 // assert(c->getZ()>0); 66 65 67 // Perform the actual Coulomb deviation 66 // Perform the actual Coulomb deviation 68 const G4bool success = coulombDeviation(c, 67 const G4bool success = coulombDeviation(c, n); 69 if(!success) { 68 if(!success) { 70 return IAvatarList(); 69 return IAvatarList(); 71 } 70 } 72 71 73 // Rely on the CoulombNone slave to comput 72 // Rely on the CoulombNone slave to compute the straight-line intersection 74 // and actually bring the particle to the 73 // and actually bring the particle to the surface of the nucleus 75 return theCoulombNoneSlave.bringToSurface( 74 return theCoulombNoneSlave.bringToSurface(c,n); 76 } 75 } 77 76 78 void CoulombNonRelativistic::distortOut(Part 77 void CoulombNonRelativistic::distortOut(ParticleList const &pL, 79 Nucleus const * const nucleus) const { 78 Nucleus const * const nucleus) const { 80 79 81 for(ParticleIter particle=pL.begin(), e=pL << 80 for(ParticleIter particle=pL.begin(); particle!=pL.end(); ++particle) { 82 81 83 const G4int Z = (*particle)->getZ(); 82 const G4int Z = (*particle)->getZ(); 84 if(Z == 0) continue; 83 if(Z == 0) continue; 85 84 86 const G4double tcos=1.-0.000001; 85 const G4double tcos=1.-0.000001; 87 86 88 const G4double et1 = PhysicalConstants:: 87 const G4double et1 = PhysicalConstants::eSquared * nucleus->getZ(); 89 const G4double transmissionRadius = 88 const G4double transmissionRadius = 90 nucleus->getDensity()->getTransmission 89 nucleus->getDensity()->getTransmissionRadius(*particle); 91 90 92 const ThreeVector position = (*particle) 91 const ThreeVector position = (*particle)->getPosition(); 93 ThreeVector momentum = (*particle)->getM 92 ThreeVector momentum = (*particle)->getMomentum(); 94 const G4double r = position.mag(); 93 const G4double r = position.mag(); 95 const G4double p = momentum.mag(); 94 const G4double p = momentum.mag(); 96 const G4double cosTheta = position.dot(m 95 const G4double cosTheta = position.dot(momentum)/(r*p); 97 if(cosTheta < 0.999999) { 96 if(cosTheta < 0.999999) { 98 const G4double sinTheta = std::sqrt(1. 97 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta); 99 const G4double eta = et1 * Z / (*parti 98 const G4double eta = et1 * Z / (*particle)->getKineticEnergy(); 100 if(eta > transmissionRadius-0.0001) { 99 if(eta > transmissionRadius-0.0001) { 101 // If below the Coulomb barrier, rad 100 // If below the Coulomb barrier, radial emission: 102 momentum = position * (p/r); 101 momentum = position * (p/r); 103 (*particle)->setMomentum(momentum); 102 (*particle)->setMomentum(momentum); 104 } else { 103 } else { 105 const G4double b0 = 0.5 * (eta + std 104 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta + 106 4. * std::pow(transmissionRadi 105 4. * std::pow(transmissionRadius*sinTheta,2) 107 * (1.-eta/transmissionRadius)) 106 * (1.-eta/transmissionRadius))); 108 const G4double bInf = std::sqrt(b0*( 107 const G4double bInf = std::sqrt(b0*(b0-eta)); 109 const G4double thr = std::atan(eta/( 108 const G4double thr = std::atan(eta/(2.*bInf)); 110 G4double uTemp = (1.-b0/transmission 109 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) + 111 b0/transmissionRadius; << 110 b0/transmissionRadius; 112 if(uTemp>tcos) uTemp=tcos; 111 if(uTemp>tcos) uTemp=tcos; 113 const G4double thd = Math::arcCos(co << 112 const G4double thd = std::acos(cosTheta)-Math::piOverTwo + thr + 114 Math::arcCos(uTemp); << 113 std::acos(uTemp); 115 const G4double c1 = std::sin(thd)*co 114 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd); 116 const G4double c2 = -p*std::sin(thd) 115 const G4double c2 = -p*std::sin(thd)/(r*sinTheta); 117 const ThreeVector newMomentum = mome 116 const ThreeVector newMomentum = momentum*c1 + position*c2; 118 (*particle)->setMomentum(newMomentum 117 (*particle)->setMomentum(newMomentum); 119 } 118 } 120 } 119 } 121 } 120 } 122 } 121 } 123 122 124 G4double CoulombNonRelativistic::maxImpactPa 123 G4double CoulombNonRelativistic::maxImpactParameter(ParticleSpecies const &p, const G4double kinE, 125 << 124 Nucleus const * const n) const { 126 const G4double theMinimumDistance = minimu << 125 G4double theMaxImpactParameter = maxImpactParameterParticle(p, kinE, n); 127 G4double rMax = n->getUniverseRadius(); << 126 if(theMaxImpactParameter <= 0.) >> 127 return 0.; 128 if(p.theType == Composite) 128 if(p.theType == Composite) 129 rMax += 2.*ParticleTable::getLargestNuc << 129 theMaxImpactParameter += 2.*ParticleTable::getNuclearRadius(p.theA, p.theZ); >> 130 return theMaxImpactParameter; >> 131 } >> 132 >> 133 G4double CoulombNonRelativistic::maxImpactParameterParticle(ParticleSpecies const &p, const G4double kinE, >> 134 Nucleus const * const n) const { >> 135 const G4double theMinimumDistance = minimumDistance(p, kinE, n); >> 136 const G4double rMax = n->getCoulombRadius(p); 130 const G4double theMaxImpactParameterSquare 137 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance); 131 if(theMaxImpactParameterSquared<=0.) 138 if(theMaxImpactParameterSquared<=0.) 132 return 0.; 139 return 0.; 133 const G4double theMaxImpactParameter = std << 140 G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared); 134 return theMaxImpactParameter; 141 return theMaxImpactParameter; 135 } 142 } 136 143 137 G4bool CoulombNonRelativistic::coulombDeviat 144 G4bool CoulombNonRelativistic::coulombDeviation(Particle * const p, Nucleus const * const n) const { 138 // Determine the rotation angle and the ne 145 // Determine the rotation angle and the new impact parameter 139 ThreeVector positionTransverse = p->getTra 146 ThreeVector positionTransverse = p->getTransversePosition(); 140 const G4double impactParameterSquared = po << 147 const G4double impactParameter = positionTransverse.mag(); 141 const G4double impactParameter = std::sqrt << 142 148 143 // Some useful variables 149 // Some useful variables 144 const G4double theMinimumDistance = minimu 150 const G4double theMinimumDistance = minimumDistance(p, n); 145 // deltaTheta2 = (pi - Rutherford scatteri 151 // deltaTheta2 = (pi - Rutherford scattering angle)/2 146 G4double deltaTheta2 = std::atan(2.*impact << 152 const G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance); 147 if(deltaTheta2<0.) << 148 deltaTheta2 += Math::pi; << 149 const G4double eccentricity = 1./std::cos( 153 const G4double eccentricity = 1./std::cos(deltaTheta2); 150 154 151 G4double newImpactParameter, alpha; // Par 155 G4double newImpactParameter, alpha; // Parameters that must be determined by the deviation 152 156 153 const G4double radius = getCoulombRadius(p << 157 ParticleSpecies aSpecies = p->getSpecies(); 154 const G4double impactParameterTangentSquar << 158 G4double kineticEnergy = p->getKineticEnergy(); 155 if(impactParameterSquared >= impactParamet << 159 // Note that in the following call to maxImpactParameter we are not 156 // The particle trajectory misses the Co << 160 // interested in the size of the cluster. This is why we call 157 // In this case the new impact parameter << 161 // maxImpactParameterParticle. 158 // approach of the hyperbola << 162 if(impactParameter>maxImpactParameterParticle(aSpecies, kineticEnergy, n)) { 159 // assert(std::abs(1. + 2.*impactParameter*imp << 163 // This should happen only for composite particles, whose trajectory can >> 164 // geometrically miss the nucleus but still trigger a cascade because of >> 165 // the finite extension of the projectile. >> 166 // In this case, the sphere radius is the minimum distance of approach >> 167 // and the kinematics is very simple. 160 newImpactParameter = 0.5 * theMinimumDis 168 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity); // the minimum distance of approach 161 alpha = Math::piOverTwo - deltaTheta2; / 169 alpha = Math::piOverTwo - deltaTheta2; // half the Rutherford scattering angle 162 } else { 170 } else { 163 // The particle trajectory intersects th 171 // The particle trajectory intersects the Coulomb sphere 164 172 165 // Compute the entrance angle 173 // Compute the entrance angle 166 const G4double argument = -(1. + 2.*impa << 174 const G4double radius = n->getCoulombRadius(p->getSpecies()); >> 175 G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance)) 167 / eccentricity; 176 / eccentricity; 168 const G4double thetaIn = Math::twoPi - M << 177 const G4double thetaIn = Math::twoPi - std::acos(argument) - deltaTheta2; 169 178 170 // Velocity angle at the entrance point 179 // Velocity angle at the entrance point 171 alpha = std::atan((1+std::cos(thetaIn)) 180 alpha = std::atan((1+std::cos(thetaIn)) 172 / (std::sqrt(eccentricity*eccentricity << 181 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn))); 173 * Math::sign(theMinimumDistance); << 174 // New impact parameter 182 // New impact parameter 175 newImpactParameter = radius * std::sin(t 183 newImpactParameter = radius * std::sin(thetaIn - alpha); 176 } 184 } 177 185 178 // Modify the impact parameter of the part 186 // Modify the impact parameter of the particle 179 positionTransverse *= newImpactParameter/p 187 positionTransverse *= newImpactParameter/positionTransverse.mag(); 180 const ThreeVector theNewPosition = p->getL 188 const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse; 181 p->setPosition(theNewPosition); 189 p->setPosition(theNewPosition); 182 190 183 // Determine the rotation axis for the inc 191 // Determine the rotation axis for the incoming particle 184 const ThreeVector &momentum = p->getMoment 192 const ThreeVector &momentum = p->getMomentum(); 185 ThreeVector rotationAxis = momentum.vector 193 ThreeVector rotationAxis = momentum.vector(positionTransverse); 186 const G4double axisLength = rotationAxis.m 194 const G4double axisLength = rotationAxis.mag(); 187 // Apply the rotation 195 // Apply the rotation 188 if(axisLength>1E-20) { 196 if(axisLength>1E-20) { 189 rotationAxis /= axisLength; 197 rotationAxis /= axisLength; 190 p->rotatePositionAndMomentum(alpha, rota << 198 p->rotate(alpha, rotationAxis); 191 } 199 } 192 200 193 return true; 201 return true; 194 } << 195 << 196 G4double CoulombNonRelativistic::getCoulombR << 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 << 205 radius = PhysicalConstants::eSquared*Z << 206 } else if(Zp==1 && Ap==3) { // t << 207 barr = 0.5*(0.5009*Math::pow23((G4doub << 208 radius = PhysicalConstants::eSquared*Z << 209 } else if(Zp==2) { // alpha, He3 << 210 barr = 0.5939*Math::pow23((G4double)At << 211 radius = PhysicalConstants::eSquared*Z << 212 } else if(Zp>2) { << 213 // Coulomb radius from the Shen model << 214 const G4double Ap13 = Math::pow13((G4d << 215 const G4double At13 = Math::pow13((G4d << 216 const G4double rp = 1.12*Ap13 - 0.94/A << 217 const G4double rt = 1.12*At13 - 0.94/A << 218 const G4double someRadius = rp+rt+3.2; << 219 const G4double theShenBarrier = Physic << 220 radius = PhysicalConstants::eSquared*Z << 221 } << 222 if(radius<=0.) { << 223 radius = ParticleTable::getLargestNucl << 224 INCL_ERROR("Negative Coulomb radius! U << 225 } << 226 INCL_DEBUG("Coulomb radius for particle << 227 << ParticleTable::getShortName(p) << 228 ", Z=" << Zt << ": " << radius << << 229 return radius; << 230 } else << 231 return n->getUniverseRadius(); << 232 } 202 } 233 203 234 } 204 } 235 205