| 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 << 32 //
>> 33 // INCL++ revision: v5.0_rc3
33 // 34 //
34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1
35 36
36 #include "globals.hh" 37 #include "globals.hh"
37 38
38 /** \file G4INCLCoulombNonRelativistic.cc 39 /** \file G4INCLCoulombNonRelativistic.cc
39 * \brief Class for non-relativistic Coulomb d 40 * \brief Class for non-relativistic Coulomb distortion.
40 * 41 *
41 * \date 14 February 2011 << 42 * Created on: 14 February 2011
42 * \author Davide Mancusi << 43 * Author: Davide Mancusi
43 */ 44 */
44 45
45 #include "G4INCLCoulombNonRelativistic.hh" 46 #include "G4INCLCoulombNonRelativistic.hh"
46 #include "G4INCLGlobals.hh" 47 #include "G4INCLGlobals.hh"
47 48
48 namespace G4INCL { 49 namespace G4INCL {
49 50
50 ParticleEntryAvatar *CoulombNonRelativistic: << 51 void CoulombNonRelativistic::bringToSurface(Particle * const p, Nucleus const * const n) const {
>> 52 ThreeVector momentumUnitVector = p->getMomentum();
>> 53 momentumUnitVector /= momentumUnitVector.mag();
>> 54
>> 55 ThreeVector positionTransverse = p->getTransversePosition();
>> 56 const G4double impactParameter = positionTransverse.mag();
>> 57
>> 58 const G4double radius = n->getSurfaceRadius(p);
>> 59
51 // No distortion for neutral particles 60 // No distortion for neutral particles
52 if(p->getZ()!=0) { << 61 G4double newImpactParameter;
53 const G4bool success = coulombDeviation( << 62 G4double alpha;
54 if(!success) // transparent << 63 if(p->getZ()==0) {
55 return NULL; << 64 newImpactParameter = impactParameter;
>> 65 alpha = 0.;
>> 66 } else {
>> 67 const G4double theCoulombFactor = coulombFactor(p, n);
>> 68 const G4double thrs2 = std::atan(theCoulombFactor/(2*impactParameter));
>> 69 const G4double eccentricity = -1./std::sin(thrs2);
>> 70 const G4double bMin = 0.5 * (
>> 71 theCoulombFactor +
>> 72 std::sqrt(theCoulombFactor*theCoulombFactor +
>> 73 4*impactParameter*impactParameter)
>> 74 );
>> 75 const G4double phyp = (1.+eccentricity) * bMin;
>> 76 const G4double thetaMax = std::acos((phyp/radius - 1.)/eccentricity);
>> 77 newImpactParameter = radius * std::cos(thrs2 + thetaMax);
>> 78 const G4double psi = std::atan( (1.+eccentricity*std::cos(thetaMax)) /
>> 79 (eccentricity*std::sin(thetaMax)));
>> 80 alpha = psi - Math::piOverTwo + thrs2 + thetaMax;
56 } 81 }
>> 82 const G4double distanceZ2 = radius*radius - newImpactParameter*newImpactParameter;
>> 83 const G4double distanceZ = (distanceZ2>0. ? std::sqrt(distanceZ2) : 0.);
57 84
58 // Rely on the CoulombNone slave to comput << 85 positionTransverse *= newImpactParameter/impactParameter;
59 // and actually bring the particle to the <<
60 return theCoulombNoneSlave.bringToSurface( <<
61 } <<
62 86
63 IAvatarList CoulombNonRelativistic::bringToS << 87 const ThreeVector position = positionTransverse - momentumUnitVector *
64 // Neutral clusters?! << 88 distanceZ;
65 // assert(c->getZ()>0); << 89 p->setPosition(position);
66 << 90
67 // Perform the actual Coulomb deviation << 91 positionTransverse /= positionTransverse.mag();
68 const G4bool success = coulombDeviation(c, << 92 const G4double momentum = p->getMomentum().mag();
69 if(!success) { << 93 const ThreeVector newMomentum = p->getMomentum() * std::cos(alpha) +
70 return IAvatarList(); << 94 positionTransverse * (std::sin(alpha) * momentum);
71 } << 95
>> 96 p->setMomentum(newMomentum);
72 97
73 // Rely on the CoulombNone slave to comput <<
74 // and actually bring the particle to the <<
75 return theCoulombNoneSlave.bringToSurface( <<
76 } 98 }
77 99
78 void CoulombNonRelativistic::distortOut(Part 100 void CoulombNonRelativistic::distortOut(ParticleList const &pL,
79 Nucleus const * const nucleus) const { 101 Nucleus const * const nucleus) const {
80 102
81 for(ParticleIter particle=pL.begin(), e=pL << 103 for(ParticleIter particle=pL.begin(); particle!=pL.end(); ++particle) {
82 104
83 const G4int Z = (*particle)->getZ(); 105 const G4int Z = (*particle)->getZ();
84 if(Z == 0) continue; 106 if(Z == 0) continue;
85 107
86 const G4double tcos=1.-0.000001; 108 const G4double tcos=1.-0.000001;
87 109
88 const G4double et1 = PhysicalConstants:: << 110 const G4double et1 = eSquared * nucleus->getZ();
89 const G4double transmissionRadius = 111 const G4double transmissionRadius =
90 nucleus->getDensity()->getTransmission 112 nucleus->getDensity()->getTransmissionRadius(*particle);
91 113
92 const ThreeVector position = (*particle) 114 const ThreeVector position = (*particle)->getPosition();
93 ThreeVector momentum = (*particle)->getM 115 ThreeVector momentum = (*particle)->getMomentum();
94 const G4double r = position.mag(); 116 const G4double r = position.mag();
95 const G4double p = momentum.mag(); 117 const G4double p = momentum.mag();
96 const G4double cosTheta = position.dot(m 118 const G4double cosTheta = position.dot(momentum)/(r*p);
97 if(cosTheta < 0.999999) { 119 if(cosTheta < 0.999999) {
98 const G4double sinTheta = std::sqrt(1. 120 const G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
99 const G4double eta = et1 * Z / (*parti 121 const G4double eta = et1 * Z / (*particle)->getKineticEnergy();
100 if(eta > transmissionRadius-0.0001) { 122 if(eta > transmissionRadius-0.0001) {
101 // If below the Coulomb barrier, rad 123 // If below the Coulomb barrier, radial emission:
102 momentum = position * (p/r); 124 momentum = position * (p/r);
103 (*particle)->setMomentum(momentum); 125 (*particle)->setMomentum(momentum);
104 } else { 126 } else {
105 const G4double b0 = 0.5 * (eta + std 127 const G4double b0 = 0.5 * (eta + std::sqrt(eta*eta +
106 4. * std::pow(transmissionRadi 128 4. * std::pow(transmissionRadius*sinTheta,2)
107 * (1.-eta/transmissionRadius)) 129 * (1.-eta/transmissionRadius)));
108 const G4double bInf = std::sqrt(b0*( 130 const G4double bInf = std::sqrt(b0*(b0-eta));
109 const G4double thr = std::atan(eta/( 131 const G4double thr = std::atan(eta/(2.*bInf));
110 G4double uTemp = (1.-b0/transmission 132 G4double uTemp = (1.-b0/transmissionRadius) * std::sin(thr) +
111 b0/transmissionRadius; << 133 b0/transmissionRadius;
112 if(uTemp>tcos) uTemp=tcos; 134 if(uTemp>tcos) uTemp=tcos;
113 const G4double thd = Math::arcCos(co << 135 const G4double thd = std::acos(cosTheta)-Math::piOverTwo + thr +
114 Math::arcCos(uTemp); << 136 std::acos(uTemp);
115 const G4double c1 = std::sin(thd)*co 137 const G4double c1 = std::sin(thd)*cosTheta/sinTheta + std::cos(thd);
116 const G4double c2 = -p*std::sin(thd) 138 const G4double c2 = -p*std::sin(thd)/(r*sinTheta);
117 const ThreeVector newMomentum = mome 139 const ThreeVector newMomentum = momentum*c1 + position*c2;
118 (*particle)->setMomentum(newMomentum 140 (*particle)->setMomentum(newMomentum);
119 } 141 }
120 } 142 }
121 } 143 }
122 } 144 }
123 145
124 G4double CoulombNonRelativistic::maxImpactPa << 146 G4double CoulombNonRelativistic::maxImpactParameter(Particle const * const p,
125 << 147 Nucleus const * const n) const {
126 const G4double theMinimumDistance = minimu << 148 const G4double theCoulombFactor = coulombFactor(p, n);
127 G4double rMax = n->getUniverseRadius(); << 149 const G4double rMax = n->getSurfaceRadius(p);
128 if(p.theType == Composite) << 150 const G4double theMaxImpactParameterSquared = rMax*(rMax-theCoulombFactor);
129 rMax += 2.*ParticleTable::getLargestNuc << 151 return (theMaxImpactParameterSquared>0. ?
130 const G4double theMaxImpactParameterSquare << 152 std::sqrt(theMaxImpactParameterSquared) : 0.);
131 if(theMaxImpactParameterSquared<=0.) <<
132 return 0.; <<
133 const G4double theMaxImpactParameter = std <<
134 return theMaxImpactParameter; <<
135 } 153 }
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 } 154 }
235 155