| 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 /*
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::Particl
56 :IAvatar(time), theParticle(aParticle), th
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->g
80 return new ReflectionChannel(theNucleus,
81 }
82
83 // We forbid transmission of resonances be
84 // delta particle below Tf can lead to neg
85 // CDPP assumes that particles stay in the
86 if(theParticle->isResonance()) {
87 const G4double theFermiEnergy = theNucle
88 if(theParticle->getKineticEnergy()<theFe
89 INCL_DEBUG("Particle " << theParticle-
90 << " Tf=" << theFermiEnergy << ",
91 return new ReflectionChannel(theNucleu
92 }
93 }
94
95 // Don't try to make a cluster if the lead
96 const G4double transmissionProbability = g
97 const G4double TOut = TMinusV;
98 const G4double kOut = particlePOut;
99 const G4double cosR = cosRefractionAngle;
100
101 INCL_DEBUG("Transmission probability for p
102 /* Don't attempt to construct clusters whe
103 * trying to escape during a nucleus-nucle
104 * this is that projectile spectators will
105 * projectile remnant, and trying to clust
106 * G4double counting. Moreover, applying t
107 * projectile spectators makes the code *r
108 * large.
109 */
110 if(theParticle->isNucleonorLambda()
111 && (!theParticle->isProjectileSpectato
112 && transmissionProbability>1.E-4) {
113 Cluster *candidateCluster = 0;
114
115 candidateCluster = Clustering::getCluste
116 if(candidateCluster != 0 &&
117 Clustering::clusterCanEscape(theNucl
118
119 INCL_DEBUG("Cluster algorithm succeede
120
121 // Check if the cluster can penetrate
122 const G4double clusterTransmissionProb
123 const G4double x = Random::shoot();
124
125 INCL_DEBUG("Transmission probability f
126
127 if (x <= clusterTransmissionProbabilit
128 theNucleus->getStore()->getBook().in
129 INCL_DEBUG("Cluster " << candidateCl
130 return new TransmissionChannel(theNu
131 } else {
132 INCL_DEBUG("Cluster " << candidateCl
133 delete candidateCluster;
134 }
135 } else {
136 delete candidateCluster;
137 }
138 }
139
140 // If we haven't transmitted a cluster (ma
141 // disabled or maybe we just can't produce
142
143 // Always transmit projectile spectators i
144 // transmission is energetically allowed
145 if(theParticle->isProjectileSpectator() &&
146 INCL_DEBUG("Particle " << theParticle->g
147 return new TransmissionChannel(theNucleu
148 }
149
150 // Transmit or reflect depending on the tr
151 const G4double x = Random::shoot();
152
153 if(x <= transmissionProbability) { // Tran
154 INCL_DEBUG("Particle " << theParticle->g
155 if(theParticle->isKaon()) theNucleus->se
156 if(theNucleus->getStore()->getConfig()->
157 return new TransmissionChannel(theNucl
158 } else {
159 return new TransmissionChannel(theNucl
160 }
161 } else { // Reflection
162 INCL_DEBUG("Particle " << theParticle->g
163 return new ReflectionChannel(theNucleus,
164 }
165 }
166
167 void SurfaceAvatar::fillFinalState(FinalStat
168 getChannel()->fillFinalState(fs);
169 }
170
171 void SurfaceAvatar::preInteraction() {}
172
173 void SurfaceAvatar::postInteraction(FinalSta
174 ParticleList const &outgoing = fs->getOutg
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<Clu
181 ParticleList const &components = clust
182 for(ParticleIter i=components.begin(),
183 if(!(*i)->isTargetSpectator())
184 theNucleus->getStore()->getBook().
185 }
186 out->setBiasCollisionVector(components
187 } else if(!theParticle->isTargetSpectato
188 // assert(out==theParticle);
189 theNucleus->getStore()->getBook().decr
190 }
191 }
192 }
193
194 std::string SurfaceAvatar::dump() const {
195 std::stringstream ss;
196 ss << "(avatar " << theTime << " 'reflecti
197 << "(list " << '\n'
198 << theParticle->dump()
199 << "))" << '\n';
200 return ss.str();
201 }
202
203 G4double SurfaceAvatar::getTransmissionProba
204
205 particleMass = particle->getMass();
206 const G4double V = particle->getPotentialE
207
208 // Correction to the particle kinetic ener
209 const G4int theA = theNucleus->getA();
210 const G4int theZ = theNucleus->getZ();
211 const G4int theS = theNucleus->getS();
212 const G4double correction = particle->getE
213 particleTOut = particle->getKineticEnergy(
214
215 if (particleTOut <= V) // No transmission
216 return 0.0;
217
218 TMinusV = particleTOut-V;
219 TMinusV2 = TMinusV*TMinusV;
220
221 // Momenta in and out
222 const G4double particlePIn2 = particle->g
223 const G4double particlePOut2 = 2.*particle
224 particlePIn = std::sqrt(particlePIn2);
225 particlePOut = std::sqrt(particlePOut2);
226
227 if (0. > V) // Automatic transmission for
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
253 // For neutral and negative particles, no
254 // Also, no Coulomb if the particle takes
255 const G4int particleZ = particle->getZ();
256 if (particleZ <= 0 || particleZ >= theZ)
257 return theTransmissionProbability;
258
259 // Nominal Coulomb barrier
260 const G4double theTransmissionBarrier = th
261 if (TMinusV >= theTransmissionBarrier) //
262 return theTransmissionProbability;
263
264 // Coulomb-penetration factor
265 const G4double px = std::sqrt(TMinusV/theT
266 const G4double logCoulombTransmission =
267 particleZ*(theZ-particleZ)/137.03*std::s
268 *(Math::arcCos(px)-px*std::sqrt(1.-px*px
269 INCL_DEBUG("Coulomb barrier, logCoulombTra
270 if (logCoulombTransmission > 35.) // Trans
271 return 0.;
272 theTransmissionProbability *= std::exp(-2.
273
274 return theTransmissionProbability;
275 }
276
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 }
301