| 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 /* 38 /*
39 * G4INCLReflectionAvatar.cc 39 * G4INCLReflectionAvatar.cc
40 * 40 *
41 * \date Jun 8, 2009 41 * \date Jun 8, 2009
42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi
43 */ 43 */
44 44
45 #include "G4INCLSurfaceAvatar.hh" 45 #include "G4INCLSurfaceAvatar.hh"
46 #include "G4INCLRandom.hh" 46 #include "G4INCLRandom.hh"
47 #include "G4INCLReflectionChannel.hh" 47 #include "G4INCLReflectionChannel.hh"
48 #include "G4INCLTransmissionChannel.hh" 48 #include "G4INCLTransmissionChannel.hh"
49 #include "G4INCLClustering.hh" 49 #include "G4INCLClustering.hh"
50 #include <sstream> 50 #include <sstream>
51 #include <string> 51 #include <string>
52 52
53 namespace G4INCL { 53 namespace G4INCL {
54 54
55 SurfaceAvatar::SurfaceAvatar(G4INCL::Particl 55 SurfaceAvatar::SurfaceAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *n)
56 :IAvatar(time), theParticle(aParticle), th 56 :IAvatar(time), theParticle(aParticle), theNucleus(n),
57 particlePIn(0.), 57 particlePIn(0.),
58 particlePOut(0.), 58 particlePOut(0.),
59 particleTOut(0.), 59 particleTOut(0.),
60 TMinusV(0.), 60 TMinusV(0.),
61 TMinusV2(0.), 61 TMinusV2(0.),
62 particleMass(0.), 62 particleMass(0.),
63 sinIncidentAngle(0.), 63 sinIncidentAngle(0.),
64 cosIncidentAngle(0.), 64 cosIncidentAngle(0.),
65 sinRefractionAngle(0.), 65 sinRefractionAngle(0.),
66 cosRefractionAngle(0.), 66 cosRefractionAngle(0.),
67 refractionIndexRatio(0.), 67 refractionIndexRatio(0.),
68 internalReflection(false) 68 internalReflection(false)
69 { 69 {
70 setType(SurfaceAvatarType); 70 setType(SurfaceAvatarType);
71 } 71 }
72 72
73 SurfaceAvatar::~SurfaceAvatar() 73 SurfaceAvatar::~SurfaceAvatar()
74 { 74 {
75 } 75 }
76 76
77 G4INCL::IChannel* SurfaceAvatar::getChannel( 77 G4INCL::IChannel* SurfaceAvatar::getChannel() {
78 if(theParticle->isTargetSpectator()) { 78 if(theParticle->isTargetSpectator()) {
79 INCL_DEBUG("Particle " << theParticle->g 79 INCL_DEBUG("Particle " << theParticle->getID() << " is a spectator, reflection" << '\n');
80 return new ReflectionChannel(theNucleus, 80 return new ReflectionChannel(theNucleus, theParticle);
81 } 81 }
82 82
83 // We forbid transmission of resonances be 83 // We forbid transmission of resonances below the Fermi energy. Emitting a
84 // delta particle below Tf can lead to neg 84 // delta particle below Tf can lead to negative excitation energies, since
85 // CDPP assumes that particles stay in the 85 // CDPP assumes that particles stay in the Fermi sea.
86 if(theParticle->isResonance()) { 86 if(theParticle->isResonance()) {
87 const G4double theFermiEnergy = theNucle 87 const G4double theFermiEnergy = theNucleus->getPotential()->getFermiEnergy(theParticle);
88 if(theParticle->getKineticEnergy()<theFe 88 if(theParticle->getKineticEnergy()<theFermiEnergy) {
89 INCL_DEBUG("Particle " << theParticle- 89 INCL_DEBUG("Particle " << theParticle->getID() << " is a resonance below Tf, reflection" << '\n'
90 << " Tf=" << theFermiEnergy << ", 90 << " Tf=" << theFermiEnergy << ", EKin=" << theParticle->getKineticEnergy() << '\n');
91 return new ReflectionChannel(theNucleu 91 return new ReflectionChannel(theNucleus, theParticle);
92 } 92 }
93 } 93 }
94 94
95 // Don't try to make a cluster if the lead 95 // Don't try to make a cluster if the leading particle is too slow
96 const G4double transmissionProbability = g 96 const G4double transmissionProbability = getTransmissionProbability(theParticle);
97 const G4double TOut = TMinusV; 97 const G4double TOut = TMinusV;
98 const G4double kOut = particlePOut; 98 const G4double kOut = particlePOut;
99 const G4double cosR = cosRefractionAngle; 99 const G4double cosR = cosRefractionAngle;
100 100
101 INCL_DEBUG("Transmission probability for p 101 INCL_DEBUG("Transmission probability for particle " << theParticle->getID() << " = " << transmissionProbability << '\n');
102 /* Don't attempt to construct clusters whe 102 /* Don't attempt to construct clusters when a projectile spectator is
103 * trying to escape during a nucleus-nucle 103 * trying to escape during a nucleus-nucleus collision. The idea behind
104 * this is that projectile spectators will 104 * this is that projectile spectators will later be collected in the
105 * projectile remnant, and trying to clust 105 * projectile remnant, and trying to clusterise them somewhat feels like
106 * G4double counting. Moreover, applying t 106 * G4double counting. Moreover, applying the clustering algorithm on escaping
107 * projectile spectators makes the code *r 107 * projectile spectators makes the code *really* slow if the projectile is
108 * large. 108 * large.
109 */ 109 */
110 if(theParticle->isNucleonorLambda() 110 if(theParticle->isNucleonorLambda()
111 && (!theParticle->isProjectileSpectato 111 && (!theParticle->isProjectileSpectator() || !theNucleus->isNucleusNucleusCollision())
112 && transmissionProbability>1.E-4) { 112 && transmissionProbability>1.E-4) {
113 Cluster *candidateCluster = 0; 113 Cluster *candidateCluster = 0;
114 114
115 candidateCluster = Clustering::getCluste 115 candidateCluster = Clustering::getCluster(theNucleus, theParticle);
116 if(candidateCluster != 0 && 116 if(candidateCluster != 0 &&
117 Clustering::clusterCanEscape(theNucl 117 Clustering::clusterCanEscape(theNucleus, candidateCluster)) {
118 118
119 INCL_DEBUG("Cluster algorithm succeede 119 INCL_DEBUG("Cluster algorithm succeeded. Candidate cluster:" << '\n' << candidateCluster->print() << '\n');
120 120
121 // Check if the cluster can penetrate 121 // Check if the cluster can penetrate the Coulomb barrier
122 const G4double clusterTransmissionProb 122 const G4double clusterTransmissionProbability = getTransmissionProbability(candidateCluster);
123 const G4double x = Random::shoot(); 123 const G4double x = Random::shoot();
124 124
125 INCL_DEBUG("Transmission probability f 125 INCL_DEBUG("Transmission probability for cluster " << candidateCluster->getID() << " = " << clusterTransmissionProbability << '\n');
126 126
127 if (x <= clusterTransmissionProbabilit 127 if (x <= clusterTransmissionProbability) {
128 theNucleus->getStore()->getBook().in 128 theNucleus->getStore()->getBook().incrementEmittedClusters();
129 INCL_DEBUG("Cluster " << candidateCl 129 INCL_DEBUG("Cluster " << candidateCluster->getID() << " passes the Coulomb barrier, transmitting." << '\n');
130 return new TransmissionChannel(theNu 130 return new TransmissionChannel(theNucleus, candidateCluster);
131 } else { 131 } else {
132 INCL_DEBUG("Cluster " << candidateCl 132 INCL_DEBUG("Cluster " << candidateCluster->getID() << " does not pass the Coulomb barrier. Falling back to transmission of the leading particle." << '\n');
133 delete candidateCluster; 133 delete candidateCluster;
134 } 134 }
135 } else { 135 } else {
136 delete candidateCluster; 136 delete candidateCluster;
137 } 137 }
138 } 138 }
139 139
140 // If we haven't transmitted a cluster (ma 140 // If we haven't transmitted a cluster (maybe cluster feature was
141 // disabled or maybe we just can't produce 141 // disabled or maybe we just can't produce an acceptable cluster):
142 142
143 // Always transmit projectile spectators i 143 // Always transmit projectile spectators if no cluster was formed and if
144 // transmission is energetically allowed 144 // transmission is energetically allowed
145 if(theParticle->isProjectileSpectator() && 145 if(theParticle->isProjectileSpectator() && transmissionProbability>0.) {
146 INCL_DEBUG("Particle " << theParticle->g 146 INCL_DEBUG("Particle " << theParticle->getID() << " is a projectile spectator, transmission" << '\n');
147 return new TransmissionChannel(theNucleu 147 return new TransmissionChannel(theNucleus, theParticle, TOut);
148 } 148 }
149 149
150 // Transmit or reflect depending on the tr 150 // Transmit or reflect depending on the transmission probability
151 const G4double x = Random::shoot(); 151 const G4double x = Random::shoot();
152 152
153 if(x <= transmissionProbability) { // Tran 153 if(x <= transmissionProbability) { // Transmission
154 INCL_DEBUG("Particle " << theParticle->g 154 INCL_DEBUG("Particle " << theParticle->getID() << " passes the Coulomb barrier, transmitting." << '\n');
155 if(theParticle->isKaon()) theNucleus->se 155 if(theParticle->isKaon()) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()-1);
156 if(theNucleus->getStore()->getConfig()-> 156 if(theNucleus->getStore()->getConfig()->getRefraction()) {
157 return new TransmissionChannel(theNucl 157 return new TransmissionChannel(theNucleus, theParticle, kOut, cosR);
158 } else { 158 } else {
159 return new TransmissionChannel(theNucl 159 return new TransmissionChannel(theNucleus, theParticle, TOut);
160 } 160 }
161 } else { // Reflection 161 } else { // Reflection
162 INCL_DEBUG("Particle " << theParticle->g 162 INCL_DEBUG("Particle " << theParticle->getID() << " does not pass the Coulomb barrier, reflection." << '\n');
163 return new ReflectionChannel(theNucleus, 163 return new ReflectionChannel(theNucleus, theParticle);
164 } 164 }
165 } 165 }
166 166
167 void SurfaceAvatar::fillFinalState(FinalStat 167 void SurfaceAvatar::fillFinalState(FinalState *fs) {
168 getChannel()->fillFinalState(fs); 168 getChannel()->fillFinalState(fs);
169 } 169 }
170 170
171 void SurfaceAvatar::preInteraction() {} 171 void SurfaceAvatar::preInteraction() {}
172 172
173 void SurfaceAvatar::postInteraction(FinalSta 173 void SurfaceAvatar::postInteraction(FinalState *fs) {
174 ParticleList const &outgoing = fs->getOutg 174 ParticleList const &outgoing = fs->getOutgoingParticles();
175 if(!outgoing.empty()) { // Transmission 175 if(!outgoing.empty()) { // Transmission
176 // assert(outgoing.size()==1); 176 // assert(outgoing.size()==1);
177 Particle *out = outgoing.front(); 177 Particle *out = outgoing.front();
178 out->rpCorrelate(); 178 out->rpCorrelate();
179 if(out->isCluster()) { 179 if(out->isCluster()) {
180 Cluster *clusterOut = dynamic_cast<Clu 180 Cluster *clusterOut = dynamic_cast<Cluster*>(out);
181 ParticleList const &components = clust 181 ParticleList const &components = clusterOut->getParticles();
182 for(ParticleIter i=components.begin(), 182 for(ParticleIter i=components.begin(), e=components.end(); i!=e; ++i) {
183 if(!(*i)->isTargetSpectator()) 183 if(!(*i)->isTargetSpectator())
184 theNucleus->getStore()->getBook(). 184 theNucleus->getStore()->getBook().decrementCascading();
185 } 185 }
186 out->setBiasCollisionVector(components 186 out->setBiasCollisionVector(components.getParticleListBiasVector());
187 } else if(!theParticle->isTargetSpectato 187 } else if(!theParticle->isTargetSpectator()) {
188 // assert(out==theParticle); 188 // assert(out==theParticle);
189 theNucleus->getStore()->getBook().decr 189 theNucleus->getStore()->getBook().decrementCascading();
190 } 190 }
191 } 191 }
192 } 192 }
193 193
194 std::string SurfaceAvatar::dump() const { 194 std::string SurfaceAvatar::dump() const {
195 std::stringstream ss; 195 std::stringstream ss;
196 ss << "(avatar " << theTime << " 'reflecti 196 ss << "(avatar " << theTime << " 'reflection" << '\n'
197 << "(list " << '\n' 197 << "(list " << '\n'
198 << theParticle->dump() 198 << theParticle->dump()
199 << "))" << '\n'; 199 << "))" << '\n';
200 return ss.str(); 200 return ss.str();
201 } 201 }
202 202
203 G4double SurfaceAvatar::getTransmissionProba 203 G4double SurfaceAvatar::getTransmissionProbability(Particle const * const particle) {
204 204
205 particleMass = particle->getMass(); 205 particleMass = particle->getMass();
206 const G4double V = particle->getPotentialE 206 const G4double V = particle->getPotentialEnergy();
207 207
208 // Correction to the particle kinetic ener 208 // Correction to the particle kinetic energy if using real masses
209 const G4int theA = theNucleus->getA(); 209 const G4int theA = theNucleus->getA();
210 const G4int theZ = theNucleus->getZ(); 210 const G4int theZ = theNucleus->getZ();
211 const G4int theS = theNucleus->getS(); 211 const G4int theS = theNucleus->getS();
212 const G4double correction = particle->getE 212 const G4double correction = particle->getEmissionQValueCorrection(theA, theZ, theS);
213 particleTOut = particle->getKineticEnergy( 213 particleTOut = particle->getKineticEnergy() + correction;
214 214
215 if (particleTOut <= V) // No transmission 215 if (particleTOut <= V) // No transmission if total energy < 0
216 return 0.0; 216 return 0.0;
217 217
218 TMinusV = particleTOut-V; 218 TMinusV = particleTOut-V;
219 TMinusV2 = TMinusV*TMinusV; 219 TMinusV2 = TMinusV*TMinusV;
220 220
221 // Momenta in and out 221 // Momenta in and out
222 const G4double particlePIn2 = particle->g 222 const G4double particlePIn2 = particle->getMomentum().mag2();
223 const G4double particlePOut2 = 2.*particle 223 const G4double particlePOut2 = 2.*particleMass*TMinusV+TMinusV2;
224 particlePIn = std::sqrt(particlePIn2); 224 particlePIn = std::sqrt(particlePIn2);
225 particlePOut = std::sqrt(particlePOut2); 225 particlePOut = std::sqrt(particlePOut2);
226 226
227 if (0. > V) // Automatic transmission for 227 if (0. > V) // Automatic transmission for repulsive potential
228 return 1.0; 228 return 1.0;
229 229
230 // Compute the transmission probability 230 // Compute the transmission probability
231 G4double theTransmissionProbability; 231 G4double theTransmissionProbability;
232 if(theNucleus->getStore()->getConfig()->ge 232 if(theNucleus->getStore()->getConfig()->getRefraction()) {
233 // Use the formula with refraction 233 // Use the formula with refraction
234 initializeRefractionVariables(particle); 234 initializeRefractionVariables(particle);
235 235
236 if(internalReflection) 236 if(internalReflection)
237 return 0.; // total internal reflectio 237 return 0.; // total internal reflection
238 238
239 // Intermediate variables for calculatio 239 // Intermediate variables for calculation
240 const G4double x = refractionIndexRatio* 240 const G4double x = refractionIndexRatio*cosIncidentAngle;
241 const G4double y = (x - cosRefractionAng 241 const G4double y = (x - cosRefractionAngle) / (x + cosRefractionAngle);
242 242
243 theTransmissionProbability = 1. - y*y; 243 theTransmissionProbability = 1. - y*y;
244 } else { 244 } else {
245 // Use the formula without refraction 245 // Use the formula without refraction
246 // Intermediate variable for calculation 246 // Intermediate variable for calculation
247 const G4double y = particlePIn+particleP 247 const G4double y = particlePIn+particlePOut;
248 248
249 // The transmission probability for a po 249 // The transmission probability for a potential step
250 theTransmissionProbability = 4.*particle 250 theTransmissionProbability = 4.*particlePIn*particlePOut/(y*y);
251 } 251 }
252 252
253 // For neutral and negative particles, no 253 // For neutral and negative particles, no Coulomb transmission
254 // Also, no Coulomb if the particle takes 254 // Also, no Coulomb if the particle takes away all of the nuclear charge
255 const G4int particleZ = particle->getZ(); 255 const G4int particleZ = particle->getZ();
256 if (particleZ <= 0 || particleZ >= theZ) 256 if (particleZ <= 0 || particleZ >= theZ)
257 return theTransmissionProbability; 257 return theTransmissionProbability;
258 258
259 // Nominal Coulomb barrier 259 // Nominal Coulomb barrier
260 const G4double theTransmissionBarrier = th 260 const G4double theTransmissionBarrier = theNucleus->getTransmissionBarrier(particle);
261 if (TMinusV >= theTransmissionBarrier) // 261 if (TMinusV >= theTransmissionBarrier) // Above the Coulomb barrier
262 return theTransmissionProbability; 262 return theTransmissionProbability;
263 263
264 // Coulomb-penetration factor 264 // Coulomb-penetration factor
265 const G4double px = std::sqrt(TMinusV/theT 265 const G4double px = std::sqrt(TMinusV/theTransmissionBarrier);
266 const G4double logCoulombTransmission = 266 const G4double logCoulombTransmission =
267 particleZ*(theZ-particleZ)/137.03*std::s 267 particleZ*(theZ-particleZ)/137.03*std::sqrt(2.*particleMass/TMinusV/(1.+TMinusV/2./particleMass))
268 *(Math::arcCos(px)-px*std::sqrt(1.-px*px 268 *(Math::arcCos(px)-px*std::sqrt(1.-px*px));
269 INCL_DEBUG("Coulomb barrier, logCoulombTra 269 INCL_DEBUG("Coulomb barrier, logCoulombTransmission=" << logCoulombTransmission << '\n');
270 if (logCoulombTransmission > 35.) // Trans 270 if (logCoulombTransmission > 35.) // Transmission is forbidden by Coulomb
271 return 0.; 271 return 0.;
272 theTransmissionProbability *= std::exp(-2. 272 theTransmissionProbability *= std::exp(-2.*logCoulombTransmission);
273 273
274 return theTransmissionProbability; 274 return theTransmissionProbability;
275 } 275 }
276 276
277 void SurfaceAvatar::initializeRefractionVari 277 void SurfaceAvatar::initializeRefractionVariables(Particle const * const particle) {
278 cosIncidentAngle = particle->getCosRPAngle 278 cosIncidentAngle = particle->getCosRPAngle();
279 if(cosIncidentAngle>1.) 279 if(cosIncidentAngle>1.)
280 cosIncidentAngle=1.; 280 cosIncidentAngle=1.;
281 sinIncidentAngle = std::sqrt(1. - cosIncid 281 sinIncidentAngle = std::sqrt(1. - cosIncidentAngle*cosIncidentAngle);
282 refractionIndexRatio = particlePIn/particl 282 refractionIndexRatio = particlePIn/particlePOut;
283 const G4double sinCandidate = refractionIn 283 const G4double sinCandidate = refractionIndexRatio*sinIncidentAngle;
284 internalReflection = (std::fabs(sinCandida 284 internalReflection = (std::fabs(sinCandidate)>1.);
285 if(internalReflection) { 285 if(internalReflection) {
286 sinRefractionAngle = 1.; 286 sinRefractionAngle = 1.;
287 cosRefractionAngle = 0.; 287 cosRefractionAngle = 0.;
288 } else { 288 } else {
289 sinRefractionAngle = sinCandidate; 289 sinRefractionAngle = sinCandidate;
290 cosRefractionAngle = std::sqrt(1. - sinR 290 cosRefractionAngle = std::sqrt(1. - sinRefractionAngle*sinRefractionAngle);
291 } 291 }
292 INCL_DEBUG("Refraction parameters initiali 292 INCL_DEBUG("Refraction parameters initialised as follows:\n"
293 << " cosIncidentAngle=" << cosIncid 293 << " cosIncidentAngle=" << cosIncidentAngle << '\n'
294 << " sinIncidentAngle=" << sinIncid 294 << " sinIncidentAngle=" << sinIncidentAngle << '\n'
295 << " cosRefractionAngle=" << cosRef 295 << " cosRefractionAngle=" << cosRefractionAngle << '\n'
296 << " sinRefractionAngle=" << sinRef 296 << " sinRefractionAngle=" << sinRefractionAngle << '\n'
297 << " refractionIndexRatio=" << refr 297 << " refractionIndexRatio=" << refractionIndexRatio << '\n'
298 << " internalReflection=" << intern 298 << " internalReflection=" << internalReflection << '\n');
299 } 299 }
300 } 300 }
301 301