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