| 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 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 /** \file G4INCLCoulombNonRelativistic.cc 38 /** \file G4INCLCoulombNonRelativistic.cc
39 * \brief Class for non-relativistic Coulomb d 39 * \brief Class for non-relativistic Coulomb distortion.
40 * 40 *
41 * \date 14 February 2011 41 * \date 14 February 2011
42 * \author Davide Mancusi 42 * \author Davide Mancusi
43 */ 43 */
44 44
45 #include "G4INCLCoulombNonRelativistic.hh" 45 #include "G4INCLCoulombNonRelativistic.hh"
46 #include "G4INCLGlobals.hh" 46 #include "G4INCLGlobals.hh"
47 47
48 namespace G4INCL { 48 namespace G4INCL {
49 49
50 ParticleEntryAvatar *CoulombNonRelativistic: 50 ParticleEntryAvatar *CoulombNonRelativistic::bringToSurface(Particle * const p, Nucleus * const n) const {
51 // No distortion for neutral particles 51 // No distortion for neutral particles
52 if(p->getZ()!=0) { 52 if(p->getZ()!=0) {
53 const G4bool success = coulombDeviation( 53 const G4bool success = coulombDeviation(p, n);
54 if(!success) // transparent 54 if(!success) // transparent
55 return NULL; 55 return NULL;
56 } 56 }
57 57
58 // Rely on the CoulombNone slave to comput 58 // Rely on the CoulombNone slave to compute the straight-line intersection
59 // and actually bring the particle to the 59 // and actually bring the particle to the surface of the nucleus
60 return theCoulombNoneSlave.bringToSurface( 60 return theCoulombNoneSlave.bringToSurface(p,n);
61 } 61 }
62 62
63 IAvatarList CoulombNonRelativistic::bringToS 63 IAvatarList CoulombNonRelativistic::bringToSurface(Cluster * const c, Nucleus * const n) const {
64 // Neutral clusters?! 64 // Neutral clusters?!
65 // assert(c->getZ()>0); 65 // assert(c->getZ()>0);
66 66
67 // Perform the actual Coulomb deviation 67 // Perform the actual Coulomb deviation
68 const G4bool success = coulombDeviation(c, 68 const G4bool success = coulombDeviation(c, n);
69 if(!success) { 69 if(!success) {
70 return IAvatarList(); 70 return IAvatarList();
71 } 71 }
72 72
73 // Rely on the CoulombNone slave to comput 73 // Rely on the CoulombNone slave to compute the straight-line intersection
74 // and actually bring the particle to the 74 // and actually bring the particle to the surface of the nucleus
75 return theCoulombNoneSlave.bringToSurface( 75 return theCoulombNoneSlave.bringToSurface(c,n);
76 } 76 }
77 77
78 void CoulombNonRelativistic::distortOut(Part 78 void CoulombNonRelativistic::distortOut(ParticleList const &pL,
79 Nucleus const * const nucleus) const { 79 Nucleus const * const nucleus) const {
80 80
81 for(ParticleIter particle=pL.begin(), e=pL 81 for(ParticleIter particle=pL.begin(), e=pL.end(); particle!=e; ++particle) {
82 82
83 const G4int Z = (*particle)->getZ(); 83 const G4int Z = (*particle)->getZ();
84 if(Z == 0) continue; 84 if(Z == 0) continue;
85 85
86 const G4double tcos=1.-0.000001; 86 const G4double tcos=1.-0.000001;
87 87
88 const G4double et1 = PhysicalConstants:: 88 const G4double et1 = PhysicalConstants::eSquared * nucleus->getZ();
89 const G4double transmissionRadius = 89 const G4double transmissionRadius =
90 nucleus->getDensity()->getTransmission 90 nucleus->getDensity()->getTransmissionRadius(*particle);
91 91
92 const ThreeVector position = (*particle) 92 const ThreeVector position = (*particle)->getPosition();
93 ThreeVector momentum = (*particle)->getM 93 ThreeVector momentum = (*particle)->getMomentum();
94 const G4double r = position.mag(); 94 const G4double r = position.mag();
95 const G4double p = momentum.mag(); 95 const G4double p = momentum.mag();
96 const G4double cosTheta = position.dot(m 96 const G4double cosTheta = position.dot(momentum)/(r*p);
97 if(cosTheta < 0.999999) { 97 if(cosTheta < 0.999999) {
98 const G4double sinTheta = std::sqrt(1. 98 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
99 const G4double eta = et1 * Z / (*parti 99 const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
100 if(eta > transmissionRadius-0.0001) { 100 if(eta > transmissionRadius-0.0001) {
101 // If below the Coulomb barrier, rad 101 // If below the Coulomb barrier, radial emission:
102 momentum = position * (p/r); 102 momentum = position * (p/r);
103 (*particle)->setMomentum(momentum); 103 (*particle)->setMomentum(momentum);
104 } else { 104 } else {
105 const G4double b0 = 0.5 * (eta + std 105 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
106 4. * std::pow(transmissionRadi 106 4. * std::pow(transmissionRadius*sinTheta,2)
107 * (1.-eta/transmissionRadius)) 107 * (1.-eta/transmissionRadius)));
108 const G4double bInf = std::sqrt(b0*( 108 const G4double bInf = std::sqrt(b0*(b0-eta));
109 const G4double thr = std::atan(eta/( 109 const G4double thr = std::atan(eta/(2.*bInf));
110 G4double uTemp = (1.-b0/transmission 110 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
111 b0/transmissionRadius; 111 b0/transmissionRadius;
112 if(uTemp>tcos) uTemp=tcos; 112 if(uTemp>tcos) uTemp=tcos;
113 const G4double thd = Math::arcCos(co 113 const G4double thd = Math::arcCos(cosTheta)-Math::piOverTwo + thr +
114 Math::arcCos(uTemp); 114 Math::arcCos(uTemp);
115 const G4double c1 = std::sin(thd)*co 115 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
116 const G4double c2 = -p*std::sin(thd) 116 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
117 const ThreeVector newMomentum = mome 117 const ThreeVector newMomentum = momentum*c1 + position*c2;
118 (*particle)->setMomentum(newMomentum 118 (*particle)->setMomentum(newMomentum);
119 } 119 }
120 } 120 }
121 } 121 }
122 } 122 }
123 123
124 G4double CoulombNonRelativistic::maxImpactPa 124 G4double CoulombNonRelativistic::maxImpactParameter(ParticleSpecies const &p, const G4double kinE,
125 125 Nucleus const * const n) const {
126 const G4double theMinimumDistance = minimu 126 const G4double theMinimumDistance = minimumDistance(p, kinE, n);
127 G4double rMax = n->getUniverseRadius(); 127 G4double rMax = n->getUniverseRadius();
128 if(p.theType == Composite) 128 if(p.theType == Composite)
129 rMax += 2.*ParticleTable::getLargestNuc 129 rMax += 2.*ParticleTable::getLargestNuclearRadius(p.theA, p.theZ);
130 const G4double theMaxImpactParameterSquare 130 const G4double theMaxImpactParameterSquared = rMax*(rMax-theMinimumDistance);
131 if(theMaxImpactParameterSquared<=0.) 131 if(theMaxImpactParameterSquared<=0.)
132 return 0.; 132 return 0.;
133 const G4double theMaxImpactParameter = std 133 const G4double theMaxImpactParameter = std::sqrt(theMaxImpactParameterSquared);
134 return theMaxImpactParameter; 134 return theMaxImpactParameter;
135 } 135 }
136 136
137 G4bool CoulombNonRelativistic::coulombDeviat 137 G4bool CoulombNonRelativistic::coulombDeviation(Particle * const p, Nucleus const * const n) const {
138 // Determine the rotation angle and the ne 138 // Determine the rotation angle and the new impact parameter
139 ThreeVector positionTransverse = p->getTra 139 ThreeVector positionTransverse = p->getTransversePosition();
140 const G4double impactParameterSquared = po 140 const G4double impactParameterSquared = positionTransverse.mag2();
141 const G4double impactParameter = std::sqrt 141 const G4double impactParameter = std::sqrt(impactParameterSquared);
142 142
143 // Some useful variables 143 // Some useful variables
144 const G4double theMinimumDistance = minimu 144 const G4double theMinimumDistance = minimumDistance(p, n);
145 // deltaTheta2 = (pi - Rutherford scatteri 145 // deltaTheta2 = (pi - Rutherford scattering angle)/2
146 G4double deltaTheta2 = std::atan(2.*impact 146 G4double deltaTheta2 = std::atan(2.*impactParameter/theMinimumDistance);
147 if(deltaTheta2<0.) 147 if(deltaTheta2<0.)
148 deltaTheta2 += Math::pi; 148 deltaTheta2 += Math::pi;
149 const G4double eccentricity = 1./std::cos( 149 const G4double eccentricity = 1./std::cos(deltaTheta2);
150 150
151 G4double newImpactParameter, alpha; // Par 151 G4double newImpactParameter, alpha; // Parameters that must be determined by the deviation
152 152
153 const G4double radius = getCoulombRadius(p 153 const G4double radius = getCoulombRadius(p->getSpecies(), n);
154 const G4double impactParameterTangentSquar 154 const G4double impactParameterTangentSquared = radius*(radius-theMinimumDistance);
155 if(impactParameterSquared >= impactParamet 155 if(impactParameterSquared >= impactParameterTangentSquared) {
156 // The particle trajectory misses the Co 156 // The particle trajectory misses the Coulomb sphere
157 // In this case the new impact parameter 157 // In this case the new impact parameter is the minimum distance of
158 // approach of the hyperbola 158 // approach of the hyperbola
159 // assert(std::abs(1. + 2.*impactParameter*imp 159 // assert(std::abs(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))>=eccentricity);
160 newImpactParameter = 0.5 * theMinimumDis 160 newImpactParameter = 0.5 * theMinimumDistance * (1.+eccentricity); // the minimum distance of approach
161 alpha = Math::piOverTwo - deltaTheta2; / 161 alpha = Math::piOverTwo - deltaTheta2; // half the Rutherford scattering angle
162 } else { 162 } else {
163 // The particle trajectory intersects th 163 // The particle trajectory intersects the Coulomb sphere
164 164
165 // Compute the entrance angle 165 // Compute the entrance angle
166 const G4double argument = -(1. + 2.*impa 166 const G4double argument = -(1. + 2.*impactParameter*impactParameter/(radius*theMinimumDistance))
167 / eccentricity; 167 / eccentricity;
168 const G4double thetaIn = Math::twoPi - M 168 const G4double thetaIn = Math::twoPi - Math::arcCos(argument) - deltaTheta2;
169 169
170 // Velocity angle at the entrance point 170 // Velocity angle at the entrance point
171 alpha = std::atan((1+std::cos(thetaIn)) 171 alpha = std::atan((1+std::cos(thetaIn))
172 / (std::sqrt(eccentricity*eccentricity 172 / (std::sqrt(eccentricity*eccentricity-1.) - std::sin(thetaIn)))
173 * Math::sign(theMinimumDistance); 173 * Math::sign(theMinimumDistance);
174 // New impact parameter 174 // New impact parameter
175 newImpactParameter = radius * std::sin(t 175 newImpactParameter = radius * std::sin(thetaIn - alpha);
176 } 176 }
177 177
178 // Modify the impact parameter of the part 178 // Modify the impact parameter of the particle
179 positionTransverse *= newImpactParameter/p 179 positionTransverse *= newImpactParameter/positionTransverse.mag();
180 const ThreeVector theNewPosition = p->getL 180 const ThreeVector theNewPosition = p->getLongitudinalPosition() + positionTransverse;
181 p->setPosition(theNewPosition); 181 p->setPosition(theNewPosition);
182 182
183 // Determine the rotation axis for the inc 183 // Determine the rotation axis for the incoming particle
184 const ThreeVector &momentum = p->getMoment 184 const ThreeVector &momentum = p->getMomentum();
185 ThreeVector rotationAxis = momentum.vector 185 ThreeVector rotationAxis = momentum.vector(positionTransverse);
186 const G4double axisLength = rotationAxis.m 186 const G4double axisLength = rotationAxis.mag();
187 // Apply the rotation 187 // Apply the rotation
188 if(axisLength>1E-20) { 188 if(axisLength>1E-20) {
189 rotationAxis /= axisLength; 189 rotationAxis /= axisLength;
190 p->rotatePositionAndMomentum(alpha, rota 190 p->rotatePositionAndMomentum(alpha, rotationAxis);
191 } 191 }
192 192
193 return true; 193 return true;
194 } 194 }
195 195
196 G4double CoulombNonRelativistic::getCoulombR 196 G4double CoulombNonRelativistic::getCoulombRadius(ParticleSpecies const &p, Nucleus const * const n) const {
197 if(p.theType == Composite) { 197 if(p.theType == Composite) {
198 const G4int Zp = p.theZ; 198 const G4int Zp = p.theZ;
199 const G4int Ap = p.theA; 199 const G4int Ap = p.theA;
200 const G4int Zt = n->getZ(); 200 const G4int Zt = n->getZ();
201 const G4int At = n->getA(); 201 const G4int At = n->getA();
202 G4double barr, radius = 0.; 202 G4double barr, radius = 0.;
203 if(Zp==1 && Ap==2) { // d 203 if(Zp==1 && Ap==2) { // d
204 barr = 0.2565*Math::pow23((G4double)At 204 barr = 0.2565*Math::pow23((G4double)At)-0.78;
205 radius = PhysicalConstants::eSquared*Z 205 radius = PhysicalConstants::eSquared*Zp*Zt/barr - 2.5;
206 } else if(Zp==1 && Ap==3) { // t 206 } else if(Zp==1 && Ap==3) { // t
207 barr = 0.5*(0.5009*Math::pow23((G4doub 207 barr = 0.5*(0.5009*Math::pow23((G4double)At)-1.16);
208 radius = PhysicalConstants::eSquared*Z 208 radius = PhysicalConstants::eSquared*Zt/barr - 0.5;
209 } else if(Zp==2) { // alpha, He3 209 } else if(Zp==2) { // alpha, He3
210 barr = 0.5939*Math::pow23((G4double)At 210 barr = 0.5939*Math::pow23((G4double)At)-1.64;
211 radius = PhysicalConstants::eSquared*Z 211 radius = PhysicalConstants::eSquared*Zp*Zt/barr - 0.5;
212 } else if(Zp>2) { 212 } else if(Zp>2) {
213 // Coulomb radius from the Shen model 213 // Coulomb radius from the Shen model
214 const G4double Ap13 = Math::pow13((G4d 214 const G4double Ap13 = Math::pow13((G4double)Ap);
215 const G4double At13 = Math::pow13((G4d 215 const G4double At13 = Math::pow13((G4double)At);
216 const G4double rp = 1.12*Ap13 - 0.94/A 216 const G4double rp = 1.12*Ap13 - 0.94/Ap13;
217 const G4double rt = 1.12*At13 - 0.94/A 217 const G4double rt = 1.12*At13 - 0.94/At13;
218 const G4double someRadius = rp+rt+3.2; 218 const G4double someRadius = rp+rt+3.2;
219 const G4double theShenBarrier = Physic 219 const G4double theShenBarrier = PhysicalConstants::eSquared*Zp*Zt/someRadius - rt*rp/(rt+rp);
220 radius = PhysicalConstants::eSquared*Z 220 radius = PhysicalConstants::eSquared*Zp*Zt/theShenBarrier;
221 } 221 }
222 if(radius<=0.) { 222 if(radius<=0.) {
223 radius = ParticleTable::getLargestNucl 223 radius = ParticleTable::getLargestNuclearRadius(Ap,Zp) + ParticleTable::getLargestNuclearRadius(At, Zt);
224 INCL_ERROR("Negative Coulomb radius! U 224 INCL_ERROR("Negative Coulomb radius! Using the sum of nuclear radii = " << radius << '\n');
225 } 225 }
226 INCL_DEBUG("Coulomb radius for particle 226 INCL_DEBUG("Coulomb radius for particle "
227 << ParticleTable::getShortName(p) 227 << ParticleTable::getShortName(p) << " in nucleus A=" << At <<
228 ", Z=" << Zt << ": " << radius << 228 ", Z=" << Zt << ": " << radius << '\n');
229 return radius; 229 return radius;
230 } else 230 } else
231 return n->getUniverseRadius(); 231 return n->getUniverseRadius();
232 } 232 }
233 233
234 } 234 }
235 235