| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
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 Helsinki Institute of Physics, Finland
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 /*
39 * G4INCLReflectionAvatar.cc
40 *
41 * \date Jun 8, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #include "G4INCLSurfaceAvatar.hh"
46 #include "G4INCLRandom.hh"
47 #include "G4INCLReflectionChannel.hh"
48 #include "G4INCLTransmissionChannel.hh"
49 #include "G4INCLClustering.hh"
50 #include <sstream>
51 #include <string>
52
53 namespace G4INCL {
54
55 SurfaceAvatar::SurfaceAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *n)
56 :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 {
70 setType(SurfaceAvatarType);
71 }
72
73 SurfaceAvatar::~SurfaceAvatar()
74 {
75 }
76
77 G4INCL::IChannel* SurfaceAvatar::getChannel() {
78 if(theParticle->isTargetSpectator()) {
79 INCL_DEBUG("Particle " << theParticle->getID() << " is a spectator, reflection" << '\n');
80 return new ReflectionChannel(theNucleus, theParticle);
81 }
82
83 // We forbid transmission of resonances below the Fermi energy. Emitting a
84 // delta particle below Tf can lead to negative excitation energies, since
85 // CDPP assumes that particles stay in the Fermi sea.
86 if(theParticle->isResonance()) {
87 const G4double theFermiEnergy = theNucleus->getPotential()->getFermiEnergy(theParticle);
88 if(theParticle->getKineticEnergy()<theFermiEnergy) {
89 INCL_DEBUG("Particle " << theParticle->getID() << " is a resonance below Tf, reflection" << '\n'
90 << " Tf=" << theFermiEnergy << ", EKin=" << theParticle->getKineticEnergy() << '\n');
91 return new ReflectionChannel(theNucleus, theParticle);
92 }
93 }
94
95 // Don't try to make a cluster if the leading particle is too slow
96 const G4double transmissionProbability = getTransmissionProbability(theParticle);
97 const G4double TOut = TMinusV;
98 const G4double kOut = particlePOut;
99 const G4double cosR = cosRefractionAngle;
100
101 INCL_DEBUG("Transmission probability for particle " << theParticle->getID() << " = " << transmissionProbability << '\n');
102 /* Don't attempt to construct clusters when a projectile spectator is
103 * trying to escape during a nucleus-nucleus collision. The idea behind
104 * this is that projectile spectators will later be collected in the
105 * projectile remnant, and trying to clusterise them somewhat feels like
106 * G4double counting. Moreover, applying the clustering algorithm on escaping
107 * projectile spectators makes the code *really* slow if the projectile is
108 * large.
109 */
110 if(theParticle->isNucleonorLambda()
111 && (!theParticle->isProjectileSpectator() || !theNucleus->isNucleusNucleusCollision())
112 && transmissionProbability>1.E-4) {
113 Cluster *candidateCluster = 0;
114
115 candidateCluster = Clustering::getCluster(theNucleus, theParticle);
116 if(candidateCluster != 0 &&
117 Clustering::clusterCanEscape(theNucleus, candidateCluster)) {
118
119 INCL_DEBUG("Cluster algorithm succeeded. Candidate cluster:" << '\n' << candidateCluster->print() << '\n');
120
121 // Check if the cluster can penetrate the Coulomb barrier
122 const G4double clusterTransmissionProbability = getTransmissionProbability(candidateCluster);
123 const G4double x = Random::shoot();
124
125 INCL_DEBUG("Transmission probability for cluster " << candidateCluster->getID() << " = " << clusterTransmissionProbability << '\n');
126
127 if (x <= clusterTransmissionProbability) {
128 theNucleus->getStore()->getBook().incrementEmittedClusters();
129 INCL_DEBUG("Cluster " << candidateCluster->getID() << " passes the Coulomb barrier, transmitting." << '\n');
130 return new TransmissionChannel(theNucleus, candidateCluster);
131 } else {
132 INCL_DEBUG("Cluster " << candidateCluster->getID() << " does not pass the Coulomb barrier. Falling back to transmission of the leading particle." << '\n');
133 delete candidateCluster;
134 }
135 } else {
136 delete candidateCluster;
137 }
138 }
139
140 // If we haven't transmitted a cluster (maybe cluster feature was
141 // disabled or maybe we just can't produce an acceptable cluster):
142
143 // Always transmit projectile spectators if no cluster was formed and if
144 // transmission is energetically allowed
145 if(theParticle->isProjectileSpectator() && transmissionProbability>0.) {
146 INCL_DEBUG("Particle " << theParticle->getID() << " is a projectile spectator, transmission" << '\n');
147 return new TransmissionChannel(theNucleus, theParticle, TOut);
148 }
149
150 // Transmit or reflect depending on the transmission probability
151 const G4double x = Random::shoot();
152
153 if(x <= transmissionProbability) { // Transmission
154 INCL_DEBUG("Particle " << theParticle->getID() << " passes the Coulomb barrier, transmitting." << '\n');
155 if(theParticle->isKaon()) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()-1);
156 if(theNucleus->getStore()->getConfig()->getRefraction()) {
157 return new TransmissionChannel(theNucleus, theParticle, kOut, cosR);
158 } else {
159 return new TransmissionChannel(theNucleus, theParticle, TOut);
160 }
161 } else { // Reflection
162 INCL_DEBUG("Particle " << theParticle->getID() << " does not pass the Coulomb barrier, reflection." << '\n');
163 return new ReflectionChannel(theNucleus, theParticle);
164 }
165 }
166
167 void SurfaceAvatar::fillFinalState(FinalState *fs) {
168 getChannel()->fillFinalState(fs);
169 }
170
171 void SurfaceAvatar::preInteraction() {}
172
173 void SurfaceAvatar::postInteraction(FinalState *fs) {
174 ParticleList const &outgoing = fs->getOutgoingParticles();
175 if(!outgoing.empty()) { // Transmission
176 // assert(outgoing.size()==1);
177 Particle *out = outgoing.front();
178 out->rpCorrelate();
179 if(out->isCluster()) {
180 Cluster *clusterOut = dynamic_cast<Cluster*>(out);
181 ParticleList const &components = clusterOut->getParticles();
182 for(ParticleIter i=components.begin(), e=components.end(); i!=e; ++i) {
183 if(!(*i)->isTargetSpectator())
184 theNucleus->getStore()->getBook().decrementCascading();
185 }
186 out->setBiasCollisionVector(components.getParticleListBiasVector());
187 } else if(!theParticle->isTargetSpectator()) {
188 // assert(out==theParticle);
189 theNucleus->getStore()->getBook().decrementCascading();
190 }
191 }
192 }
193
194 std::string SurfaceAvatar::dump() const {
195 std::stringstream ss;
196 ss << "(avatar " << theTime << " 'reflection" << '\n'
197 << "(list " << '\n'
198 << theParticle->dump()
199 << "))" << '\n';
200 return ss.str();
201 }
202
203 G4double SurfaceAvatar::getTransmissionProbability(Particle const * const particle) {
204
205 particleMass = particle->getMass();
206 const G4double V = particle->getPotentialEnergy();
207
208 // Correction to the particle kinetic energy if using real masses
209 const G4int theA = theNucleus->getA();
210 const G4int theZ = theNucleus->getZ();
211 const G4int theS = theNucleus->getS();
212 const G4double correction = particle->getEmissionQValueCorrection(theA, theZ, theS);
213 particleTOut = particle->getKineticEnergy() + correction;
214
215 if (particleTOut <= V) // No transmission if total energy < 0
216 return 0.0;
217
218 TMinusV = particleTOut-V;
219 TMinusV2 = TMinusV*TMinusV;
220
221 // Momenta in and out
222 const G4double particlePIn2 = particle->getMomentum().mag2();
223 const G4double particlePOut2 = 2.*particleMass*TMinusV+TMinusV2;
224 particlePIn = std::sqrt(particlePIn2);
225 particlePOut = std::sqrt(particlePOut2);
226
227 if (0. > V) // Automatic transmission for repulsive potential
228 return 1.0;
229
230 // Compute the transmission probability
231 G4double theTransmissionProbability;
232 if(theNucleus->getStore()->getConfig()->getRefraction()) {
233 // Use the formula with refraction
234 initializeRefractionVariables(particle);
235
236 if(internalReflection)
237 return 0.; // total internal reflection
238
239 // Intermediate variables for calculation
240 const G4double x = refractionIndexRatio*cosIncidentAngle;
241 const G4double y = (x - cosRefractionAngle) / (x + cosRefractionAngle);
242
243 theTransmissionProbability = 1. - y*y;
244 } else {
245 // Use the formula without refraction
246 // Intermediate variable for calculation
247 const G4double y = particlePIn+particlePOut;
248
249 // The transmission probability for a potential step
250 theTransmissionProbability = 4.*particlePIn*particlePOut/(y*y);
251 }
252
253 // For neutral and negative particles, no Coulomb transmission
254 // Also, no Coulomb if the particle takes away all of the nuclear charge
255 const G4int particleZ = particle->getZ();
256 if (particleZ <= 0 || particleZ >= theZ)
257 return theTransmissionProbability;
258
259 // Nominal Coulomb barrier
260 const G4double theTransmissionBarrier = theNucleus->getTransmissionBarrier(particle);
261 if (TMinusV >= theTransmissionBarrier) // Above the Coulomb barrier
262 return theTransmissionProbability;
263
264 // Coulomb-penetration factor
265 const G4double px = std::sqrt(TMinusV/theTransmissionBarrier);
266 const G4double logCoulombTransmission =
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));
269 INCL_DEBUG("Coulomb barrier, logCoulombTransmission=" << logCoulombTransmission << '\n');
270 if (logCoulombTransmission > 35.) // Transmission is forbidden by Coulomb
271 return 0.;
272 theTransmissionProbability *= std::exp(-2.*logCoulombTransmission);
273
274 return theTransmissionProbability;
275 }
276
277 void SurfaceAvatar::initializeRefractionVariables(Particle const * const particle) {
278 cosIncidentAngle = particle->getCosRPAngle();
279 if(cosIncidentAngle>1.)
280 cosIncidentAngle=1.;
281 sinIncidentAngle = std::sqrt(1. - cosIncidentAngle*cosIncidentAngle);
282 refractionIndexRatio = particlePIn/particlePOut;
283 const G4double sinCandidate = refractionIndexRatio*sinIncidentAngle;
284 internalReflection = (std::fabs(sinCandidate)>1.);
285 if(internalReflection) {
286 sinRefractionAngle = 1.;
287 cosRefractionAngle = 0.;
288 } else {
289 sinRefractionAngle = sinCandidate;
290 cosRefractionAngle = std::sqrt(1. - sinRefractionAngle*sinRefractionAngle);
291 }
292 INCL_DEBUG("Refraction parameters initialised as follows:\n"
293 << " cosIncidentAngle=" << cosIncidentAngle << '\n'
294 << " sinIncidentAngle=" << sinIncidentAngle << '\n'
295 << " cosRefractionAngle=" << cosRefractionAngle << '\n'
296 << " sinRefractionAngle=" << sinRefractionAngle << '\n'
297 << " refractionIndexRatio=" << refractionIndexRatio << '\n'
298 << " internalReflection=" << internalReflection << '\n');
299 }
300 }
301