| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /** \file G4INCLCoulombNonRelativistic.cc
39 * \brief Class for non-relativistic Coulomb d
40 *
41 * \date 14 February 2011
42 * \author Davide Mancusi
43 */
44
45 #include "G4INCLCoulombNonRelativistic.hh"
46 #include "G4INCLGlobals.hh"
47
48 namespace G4INCL {
49
50 ParticleEntryAvatar *CoulombNonRelativistic:
51 // No distortion for neutral particles
52 if(p->getZ()!=0) {
53 const G4bool success = coulombDeviation(
54 if(!success) // transparent
55 return NULL;
56 }
57
58 // Rely on the CoulombNone slave to comput
59 // and actually bring the particle to the
60 return theCoulombNoneSlave.bringToSurface(
61 }
62
63 IAvatarList CoulombNonRelativistic::bringToS
64 // Neutral clusters?!
65 // assert(c->getZ()>0);
66
67 // Perform the actual Coulomb deviation
68 const G4bool success = coulombDeviation(c,
69 if(!success) {
70 return IAvatarList();
71 }
72
73 // Rely on the CoulombNone slave to comput
74 // and actually bring the particle to the
75 return theCoulombNoneSlave.bringToSurface(
76 }
77
78 void CoulombNonRelativistic::distortOut(Part
79 Nucleus const * const nucleus) const {
80
81 for(ParticleIter particle=pL.begin(), e=pL
82
83 const G4int Z = (*particle)->getZ();
84 if(Z == 0) continue;
85
86 const G4double tcos=1.-0.000001;
87
88 const G4double et1 = PhysicalConstants::
89 const G4double transmissionRadius =
90 nucleus->getDensity()->getTransmission
91
92 const ThreeVector position = (*particle)
93 ThreeVector momentum = (*particle)->getM
94 const G4double r = position.mag();
95 const G4double p = momentum.mag();
96 const G4double cosTheta = position.dot(m
97 if(cosTheta < 0.999999) {
98 const G4double sinTheta = std::sqrt(1.
99 const G4double eta = et1 * Z / (*parti
100 if(eta > transmissionRadius-0.0001) {
101 // If below the Coulomb barrier, rad
102 momentum = position * (p/r);
103 (*particle)->setMomentum(momentum);
104 } else {
105 const G4double b0 = 0.5 * (eta + std
106 4. * std::pow(transmissionRadi
107 * (1.-eta/transmissionRadius))
108 const G4double bInf = std::sqrt(b0*(
109 const G4double thr = std::atan(eta/(
110 G4double uTemp = (1.-b0/transmission
111 b0/transmissionRadius;
112 if(uTemp>tcos) uTemp=tcos;
113 const G4double thd = Math::arcCos(co
114 Math::arcCos(uTemp);
115 const G4double c1 = std::sin(thd)*co
116 const G4double c2 = -p*std::sin(thd)
117 const ThreeVector newMomentum = mome
118 (*particle)->setMomentum(newMomentum
119 }
120 }
121 }
122 }
123
124 G4double CoulombNonRelativistic::maxImpactPa
125
126 const G4double theMinimumDistance = minimu
127 G4double rMax = n->getUniverseRadius();
128 if(p.theType == Composite)
129 rMax += 2.*ParticleTable::getLargestNuc
130 const G4double theMaxImpactParameterSquare
131 if(theMaxImpactParameterSquared<=0.)
132 return 0.;
133 const G4double theMaxImpactParameter = std
134 return theMaxImpactParameter;
135 }
136
137 G4bool CoulombNonRelativistic::coulombDeviat
138 // Determine the rotation angle and the ne
139 ThreeVector positionTransverse = p->getTra
140 const G4double impactParameterSquared = po
141 const G4double impactParameter = std::sqrt
142
143 // Some useful variables
144 const G4double theMinimumDistance = minimu
145 // deltaTheta2 = (pi - Rutherford scatteri
146 G4double deltaTheta2 = std::atan(2.*impact
147 if(deltaTheta2<0.)
148 deltaTheta2 += Math::pi;
149 const G4double eccentricity = 1./std::cos(
150
151 G4double newImpactParameter, alpha; // Par
152
153 const G4double radius = getCoulombRadius(p
154 const G4double impactParameterTangentSquar
155 if(impactParameterSquared >= impactParamet
156 // The particle trajectory misses the Co
157 // In this case the new impact parameter
158 // approach of the hyperbola
159 // assert(std::abs(1. + 2.*impactParameter*imp
160 newImpactParameter = 0.5 * theMinimumDis
161 alpha = Math::piOverTwo - deltaTheta2; /
162 } else {
163 // The particle trajectory intersects th
164
165 // Compute the entrance angle
166 const G4double argument = -(1. + 2.*impa
167 / eccentricity;
168 const G4double thetaIn = Math::twoPi - M
169
170 // Velocity angle at the entrance point
171 alpha = std::atan((1+std::cos(thetaIn))
172 / (std::sqrt(eccentricity*eccentricity
173 * Math::sign(theMinimumDistance);
174 // New impact parameter
175 newImpactParameter = radius * std::sin(t
176 }
177
178 // Modify the impact parameter of the part
179 positionTransverse *= newImpactParameter/p
180 const ThreeVector theNewPosition = p->getL
181 p->setPosition(theNewPosition);
182
183 // Determine the rotation axis for the inc
184 const ThreeVector &momentum = p->getMoment
185 ThreeVector rotationAxis = momentum.vector
186 const G4double axisLength = rotationAxis.m
187 // Apply the rotation
188 if(axisLength>1E-20) {
189 rotationAxis /= axisLength;
190 p->rotatePositionAndMomentum(alpha, rota
191 }
192
193 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 }
233
234 }
235