33 #define INCLXX_IN_GEANT4_MODE 1
55 :
IAvatar(time), theParticle(aParticle), theNucleus(n),
64 sinRefractionAngle(0.),
65 cosRefractionAngle(0.),
66 refractionIndexRatio(0.),
67 internalReflection(false)
78 INCL_DEBUG(
"Particle " << theParticle->
getID() <<
" is a spectator, reflection" << std::endl);
88 INCL_DEBUG(
"Particle " << theParticle->
getID() <<
" is a resonance below Tf, reflection" << std::endl
89 <<
" Tf=" << theFermiEnergy <<
", EKin=" << theParticle->
getKineticEnergy() << std::endl);
98 const G4double cosR = cosRefractionAngle;
100 INCL_DEBUG(
"Transmission probability for particle " << theParticle->
getID() <<
" = " << transmissionProbability << std::endl);
111 && transmissionProbability>1.E-4) {
115 if(candidateCluster != 0 &&
118 INCL_DEBUG(
"Cluster algorithm succeded. Candidate cluster:" << std::endl << candidateCluster->
print() << std::endl);
124 INCL_DEBUG(
"Transmission probability for cluster " << candidateCluster->
getID() <<
" = " << clusterTransmissionProbability << std::endl);
126 if (x <= clusterTransmissionProbability) {
128 INCL_DEBUG(
"Cluster " << candidateCluster->
getID() <<
" passes the Coulomb barrier, transmitting." << std::endl);
131 INCL_DEBUG(
"Cluster " << candidateCluster->
getID() <<
" does not pass the Coulomb barrier. Falling back to transmission of the leading particle." << std::endl);
132 delete candidateCluster;
135 delete candidateCluster;
145 INCL_DEBUG(
"Particle " << theParticle->
getID() <<
" is a projectile spectator, transmission" << std::endl);
152 if(x <= transmissionProbability) {
153 INCL_DEBUG(
"Particle " << theParticle->
getID() <<
" passes the Coulomb barrier, transmitting." << std::endl);
160 INCL_DEBUG(
"Particle " << theParticle->
getID() <<
" does not pass the Coulomb barrier, reflection." << std::endl);
173 if(!outgoing.empty()) {
180 for(
ParticleIter i=components.begin(), e=components.end(); i!=e; ++i) {
181 if(!(*i)->isTargetSpectator())
193 std::stringstream ss;
194 ss <<
"(avatar " <<
theTime <<
" 'reflection" << std::endl
195 <<
"(list " << std::endl
196 << theParticle->
dump()
197 <<
"))" << std::endl;
203 particleMass = particle->
getMass();
212 if (particleTOut <= V)
215 TMinusV = particleTOut-V;
216 TMinusV2 = TMinusV*TMinusV;
220 const G4double particlePOut2 = 2.*particleMass*TMinusV+TMinusV2;
221 particlePIn = std::sqrt(particlePIn2);
222 particlePOut = std::sqrt(particlePOut2);
225 G4double theTransmissionProbability;
228 initializeRefractionVariables(particle);
230 if(internalReflection)
234 const G4double x = refractionIndexRatio*cosIncidentAngle;
235 const G4double y = (x - cosRefractionAngle) / (x + cosRefractionAngle);
237 theTransmissionProbability = 1. - y*
y;
242 const G4double y = particlePIn+particlePOut;
245 theTransmissionProbability = 4.*particlePIn*particlePOut/(y*
y);
250 const G4int particleZ = particle->
getZ();
251 if (particleZ <= 0 || particleZ >= theZ)
252 return theTransmissionProbability;
256 if (TMinusV >= theTransmissionBarrier)
257 return theTransmissionProbability;
260 const G4double px = std::sqrt(TMinusV/theTransmissionBarrier);
261 const G4double logCoulombTransmission =
262 particleZ*(theZ-particleZ)/137.03*std::sqrt(2.*particleMass/TMinusV/(1.+TMinusV/2./particleMass))
263 *(std::acos(px)-px*std::sqrt(1.-px*px));
264 INCL_DEBUG(
"Coulomb barrier, logCoulombTransmission=" << logCoulombTransmission << std::endl);
265 if (logCoulombTransmission > 35.)
267 theTransmissionProbability *= std::exp(-2.*logCoulombTransmission);
269 return theTransmissionProbability;
272 void SurfaceAvatar::initializeRefractionVariables(
Particle const *
const particle) {
274 if(cosIncidentAngle>1.)
276 sinIncidentAngle = std::sqrt(1. - cosIncidentAngle*cosIncidentAngle);
277 refractionIndexRatio = particlePIn/particlePOut;
278 const G4double sinCandidate = refractionIndexRatio*sinIncidentAngle;
279 internalReflection = (std::fabs(sinCandidate)>1.);
280 if(internalReflection) {
281 sinRefractionAngle = 1.;
282 cosRefractionAngle = 0.;
284 sinRefractionAngle = sinCandidate;
285 cosRefractionAngle = std::sqrt(1. - sinRefractionAngle*sinRefractionAngle);
287 INCL_DEBUG(
"Refraction parameters initialised as follows:\n"
288 <<
" cosIncidentAngle=" << cosIncidentAngle <<
'\n'
289 <<
" sinIncidentAngle=" << sinIncidentAngle <<
'\n'
290 <<
" cosRefractionAngle=" << cosRefractionAngle <<
'\n'
291 <<
" sinRefractionAngle=" << sinRefractionAngle <<
'\n'
292 <<
" refractionIndexRatio=" << refractionIndexRatio <<
'\n'
293 <<
" internalReflection=" << internalReflection << std::endl);
G4int getA() const
Returns the baryon number.
G4double getCosRPAngle() const
Get the cosine of the angle between position and momentum.
G4bool isResonance() const
Is it a resonance?
G4INCL::IChannel * getChannel()
G4double getMass() const
Get the cached particle mass.
virtual FinalState * postInteraction(FinalState *)
ParticleList const & getParticles() const
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
G4INCL::FinalState * getFinalState()
G4bool isTargetSpectator() const
const G4INCL::ThreeVector & getMomentum() const
Config const * getConfig()
G4double getTransmissionProbability(Particle const *const particle)
Calculate the transmission probability for the particle.
G4double getPotentialEnergy() const
Get the particle potential energy.
static G4bool clusterCanEscape(Nucleus const *const n, Cluster const *const c)
G4double getFermiEnergy(const Particle *const p) const
Return the Fermi energy for a particle.
void rpCorrelate()
Make the particle follow a strict r-p correlation.
ParticleList const & getOutgoingParticles() const
void setType(AvatarType t)
SurfaceAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *aNucleus)
std::string print() const
G4int getZ() const
Returns the charge number.
void incrementEmittedClusters()
NuclearPotential::INuclearPotential const * getPotential() const
Getter for thePotential.
void decrementCascading()
Static class for cluster formation.
G4double getTransmissionBarrier(Particle const *const p)
Get the transmission barrier.
static Cluster * getCluster(Nucleus *n, Particle *p)
G4bool isProjectileSpectator() const
virtual G4INCL::FinalState * getFinalState()=0
virtual void preInteraction()
G4double getKineticEnergy() const
Get the particle kinetic energy.
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
ParticleList::const_iterator ParticleIter
G4bool getRefraction() const
True if we should use refraction.