33 #define INCLXX_IN_GEANT4_MODE 1
63 :theNucleus(0), maximumTime(70.0), currentTime(0.0), firstAvatar(true),
64 theLocalEnergyType(localEnergyType),
65 theLocalEnergyDeltaType(localEnergyDeltaType)
81 return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
93 G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
99 temfin = 29.8 * std::pow(theNucleus->
getA(), 0.16);
102 temfin = 30.18 * std::pow(theNucleus->
getA(), 0.17*(1.0 - 5.7E-5*tlab));
105 maximumTime = temfin;
111 const G4double traversalTime = distance / projectileVelocity;
112 if(maximumTime < traversalTime)
113 maximumTime = traversalTime;
114 INCL_DEBUG(
"Cascade stopping time is " << maximumTime << std::endl);
120 INCL_DEBUG(
"impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
126 impactParameter * std::sin(phi),
168 maximumTime = 29.8 * std::pow(theNucleus->
getA(), 0.16);
173 const G4double distance = 2.*rMax + 2.725*rms;
175 const G4double traversalTime = distance / projectileVelocity;
176 if(maximumTime < traversalTime)
177 maximumTime = traversalTime;
178 INCL_DEBUG(
"Cascade stopping time is " << maximumTime << std::endl);
184 INCL_DEBUG(
"impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
191 impactParameter * std::sin(phi),
208 if(theAvatarList.empty()) {
209 INCL_DEBUG(
"No ParticleEntryAvatar found, transparent event" << std::endl);
253 theNucleus = nucleus;
258 if(anAvatar) theNucleus->
getStore()->
add(anAvatar);
270 G4double minDistOfApproachSquared = 0.0;
272 if(t>maximumTime || t<currentTime)
return NULL;
285 const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->
isPion());
286 const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->
isPion());
288 if(p1HasLocalEnergy) {
289 backupParticle1 = *p1;
292 *p1 = backupParticle1;
297 if(p2HasLocalEnergy) {
298 backupParticle2 = *p2;
301 *p2 = backupParticle2;
302 if(p1HasLocalEnergy) {
303 *p1 = backupParticle1;
315 if(p1HasLocalEnergy) {
316 *p1 = backupParticle1;
318 if(p2HasLocalEnergy) {
319 *p2 = backupParticle2;
328 if(
Math::tenPi*minDistOfApproachSquared > totalCrossSection)
return NULL;
344 if(theIntersection.
exists) {
345 time = currentTime + theIntersection.
time;
347 INCL_ERROR(
"Imaginary reflection time for particle: " << std::endl
348 << aParticle->
print());
365 (*minDistOfApproach) = 100000.0;
366 return currentTime + 100000.0;
370 (*minDistOfApproach) = distance.
mag2() + time * t7;
371 return currentTime + time;
377 for(
ParticleIter updated=updatedParticles.begin(), e=updatedParticles.end(); updated!=e; ++updated)
380 for(
ParticleIter particle=particles.begin(), end=particles.end(); particle!=end; ++particle)
386 if((*particle)->isInList(updatedParticles))
continue;
395 for(
ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1) {
397 for(
ParticleIter p2 = p1 + 1; p2 != particles.end(); ++p2) {
405 const G4bool haveExcept = !except.empty();
408 for(
ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1)
412 for(++p2; p2 != particles.end(); ++p2)
415 if(haveExcept && (*p1)->isInList(except) && (*p2)->isInList(except))
continue;
425 for(
ParticleIter iter=particles.begin(), e=particles.end(); iter!=e; ++iter) {
436 for(
ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
444 #ifdef INCL_REGENERATE_AVATARS
445 void StandardPropagationModel::generateAllAvatarsExceptUpdated() {
448 for(
ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
459 for(
ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
460 if((*i)->isDelta()) {
462 G4double time = currentTime + decayTime;
463 if(time <= maximumTime) {
474 #ifdef INCL_REGENERATE_AVATARS
475 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
481 generateAllAvatarsExceptUpdated();
489 updatedParticles.push_back(blockedDelta);
494 needNewAvatars.insert(needNewAvatars.end(), created.begin(), created.end());
499 if(theAvatar == 0)
return 0;
502 if(theAvatar->
getTime() < currentTime) {
503 INCL_ERROR(
"Avatar time = " << theAvatar->
getTime() <<
", currentTime = " << currentTime << std::endl);
505 }
else if(theAvatar->
getTime() > currentTime) {
508 currentTime = theAvatar->
getTime();
515 void StandardPropagationModel::putSpectatorsOnShell(
IAvatarList const &entryAvatars,
ParticleList const &spectators) {
517 for(
ParticleIter p=spectators.begin(), e=spectators.end();
p!=e; ++
p) {
519 const G4double oldEnergy = (*p)->getEnergy();
520 (*p)->setTableMass();
521 (*p)->adjustEnergyFromMomentum();
522 deltaE += (*p)->getEnergy() - oldEnergy;
525 deltaE /= entryAvatars.size();
527 for(
IAvatarIter a=entryAvatars.begin(), e=entryAvatars.end();
a!=e; ++
a) {
529 Particle *
p = (*a)->getParticles().front();
534 p->setEnergy(energy);
535 const G4double newMass = std::sqrt(energy*energy - p->getMomentum().mag2());
G4int getA() const
Returns the baryon number.
void registerAvatar(G4INCL::IAvatar *anAvatar)
G4bool isResonance() const
Is it a resonance?
G4double shootComposite(ParticleSpecies const s, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
void generateAllAvatars()
(Re)Generate all possible avatars.
G4double getMass() const
Get the cached particle mass.
void setIncomingAngularMomentum(const ThreeVector &j)
Set the incoming angular-momentum vector.
void setIncomingMomentum(const ThreeVector &p)
Set the incoming momentum vector.
G4ThreadLocal ParticleMassFn getTableParticleMass
Static pointer to the mass function for particles.
G4double dot(const ThreeVector &v) const
ParticleList const & getParticles() const
void transformToLocalEnergyFrame(Nucleus const *const n, Particle *const p)
G4int getAcceptedCollisions() const
ThreeVector getPropagationVelocity() const
Get the propagation velocity of the particle.
void setInitialEnergy(const G4double e)
Set the initial energy.
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType)
void generateCollisions(const ParticleList &particles)
Generate and register collisions among particles in a list, except between those in another list...
UnorderedVector< IAvatar * >::const_iterator IAvatarIter
IAvatar * findSmallestTime()
const G4INCL::ThreeVector & getMomentum() const
void setParticleNucleusCollision()
Set a particle-nucleus collision.
std::string print() const
G4double getReflectionTime(G4INCL::Particle const *const aParticle)
Get the reflection time.
void updateAvatars(const ParticleList &particles)
void storeComponents()
Store the projectile components.
G4double getEnergy() const
void setINCLMass()
Set the mass of the Particle to its table mass.
virtual void makeProjectileSpectator()
void propagate(G4double step)
void setNucleus(G4INCL::Nucleus *nucleus)
G4double getSurfaceRadius(Particle const *const particle) const
Get the maximum allowed radius for a given particle.
G4bool isParticipant() const
ThreeVector boostVector() const
void setStoppingTime(G4double)
static ParticleEntryAvatar * bringToSurface(Particle *p, Nucleus *const n)
Modify the momentum of an incoming particle and position it on the surface of the nucleus...
G4double getTime(G4INCL::Particle const *const particleA, G4INCL::Particle const *const particleB, G4double *minDistOfApproach) const
void setEnergy(G4double energy)
Intersection getLaterTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r)
Compute the second intersection of a straight particle trajectory with a sphere.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
static G4double computeDecayTime(Particle *p)
void addParticleEntryAvatars(IAvatarList const &al)
Add one ParticleEntry avatar.
G4double shoot(ParticleSpecies const projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
double precision function energy(A, Z)
Static class for selecting Coulomb distortion.
UnorderedVector< Particle * > ParticleList
static G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
void initialiseParticleAvatarConnections()
Initialise the particleAvatarConnections map.
ParticipantType getParticipantType() const
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
static G4double getCutNNSquared()
virtual G4INCL::ThreeVector getAngularMomentum() const
virtual void setPosition(const G4INCL::ThreeVector &position)
G4double getCurrentTime()
G4int getZ() const
Returns the charge number.
ParticleList const & getUpdatedParticles() const
ParticleList const & getCreatedParticles() const
ThreeVector getTransversePosition() const
Transverse component of the position w.r.t. the momentum.
IAvatar * generateBinaryCollisionAvatar(Particle *const p1, Particle *const p2)
Generate a two-particle avatar.
void generateDecays(const ParticleList &particles)
Generate decays for particles that can decay.
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
void setProjectileChargeNumber(G4int n)
Set the charge number of the projectile.
G4INCL::Nucleus * getNucleus()
G4double total(Particle const *const p1, Particle const *const p2)
G4double getINCLMass(const G4int A, const G4int Z)
Get INCL nuclear mass (in MeV/c^2)
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
void setProjectileRemnant(ProjectileRemnant *const c)
Set the projectile remnant.
const G4INCL::ThreeVector & getPosition() const
void setProjectileMassNumber(G4int n)
Set the mass number of the projectile.
void setCurrentTime(G4double t)
virtual G4INCL::ParticleSpecies getSpecies() const
Get the particle species.
G4double getUniverseRadius() const
Getter for theUniverseRadius.
Intersection-point structure.
virtual ~StandardPropagationModel()
G4double shootParticle(ParticleType const t, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)
void generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles)
Generate and register collisions between a list of updated particles and all the other particles...
Particle * getBlockedDelta() const
Get the delta that could not decay.
G4bool isPion() const
Is this a pion?
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
void setNucleusNucleusCollision()
Set a nucleus-nucleus collision.
G4double getStoppingTime()
void timeStep(G4double step)
void setPosition(const ThreeVector &position)
Set the position of the cluster.
ParticleList::const_iterator ParticleIter
Simple class for computing intersections between a straight line and a sphere.
G4INCL::IAvatar * propagate()