34 #define INCLXX_IN_GEANT4_MODE 1
45 #ifndef G4INCLNucleus_hh
46 #define G4INCLNucleus_hh 1
67 theInitialZ(charge), theInitialA(mass),
68 theNpInitial(0), theNnInitial(0),
69 initialInternalEnergy(0.),
70 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
71 initialCenterOfMass(0.,0.,0.),
75 theUniverseRadius(universeRadius),
76 isNucleusNucleus(false),
77 theProjectileRemnant(NULL),
141 for(
ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
143 if(!(*iter)->isOutOfWell()) {
144 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
149 for(
ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
154 for(
ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
156 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
160 for(
ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
161 if((*iter)->isCluster()) {
164 #ifdef INCLXX_IN_GEANT4_MODE
169 for(
ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
174 totalEnergy += (*iter)->getEnergy();
175 theA -= (*iter)->getA();
176 theZ -= (*iter)->getZ();
182 for(
ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
184 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
190 INCL_DEBUG(
"A Particle is entering below the Fermi sea:" <<
'\n' << finalstate->
print() <<
'\n');
193 for(
ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
200 INCL_ERROR(
"Energy nonconservation! Energy at the beginning of the event = "
202 <<
" and after interaction = "
203 << totalEnergy <<
'\n'
204 << finalstate->
print());
209 INCL_WARN(
"Useless Nucleus::propagateParticles -method called." <<
'\n');
215 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
216 if((*p)->isNucleon())
217 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
218 else if((*p)->isResonance())
221 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
241 for(
ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p) {
243 theSpin -= (*p)->getAngularMomentum();
263 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
264 const G4double mass = (*p)->getMass();
265 cm += (*p)->getPosition() * mass;
281 std::stringstream ss;
282 ss <<
"Particles in the nucleus:" <<
'\n'
283 <<
"Inside:" <<
'\n';
286 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
287 ss <<
"index = " << counter <<
'\n'
291 ss <<
"Outgoing:" <<
'\n';
293 for(
ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p)
302 for(
ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
303 if((*i)->isDelta()) deltas.push_back((*i));
305 if(deltas.empty())
return false;
307 for(
ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
308 INCL_DEBUG(
"Decay outgoing delta particle:" <<
'\n'
309 << (*i)->print() <<
'\n');
311 const G4double deltaMass = (*i)->getMass();
316 (*i)->setEnergy((*i)->getMass());
329 newMomentum *= decayMomentum / newMomentum.
mag();
340 nucleon->
boost(beta);
359 const G4bool unphysicalRemnant = (theZ<0 || theZ>
theA);
366 for(
ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
367 if((*i)->isDelta()) deltas.push_back((*i));
370 for(
ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
371 INCL_DEBUG(
"Decay inside delta particle:" <<
'\n'
372 << (*i)->print() <<
'\n');
378 if(unphysicalRemnant) {
379 INCL_WARN(
"Forcing delta decay inside an unphysical remnant (A=" <<
theA
380 <<
", Z=" <<
theZ <<
"). Might lead to energy-violation warnings."
399 if(unphysicalRemnant) {
400 INCL_DEBUG(
"Remnant is unphysical: Z=" <<
theZ <<
", A=" <<
theA <<
", emitting all the pions" <<
'\n');
410 for(
ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
411 if((*i)->isCluster()) clusters.push_back((*i));
413 if(clusters.empty())
return false;
415 for(
ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
418 #ifdef INCLXX_IN_GEANT4_MODE
424 for(
ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j)
436 for(
ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j)
447 INCL_WARN(
"Forcing emissions of all pions in the nucleus." <<
'\n');
450 const G4double tinyPionEnergy = 0.1;
455 for(
ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
458 INCL_DEBUG(
"Forcing emission of the following particle: "
459 << thePion->
print() <<
'\n');
465 if(kineticEnergyOutside > 0.0)
472 toEject.push_back(thePion);
475 for(
ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
487 if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
515 for(
ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
523 if(outgoing.size() == 2) {
525 INCL_DEBUG(
"Two particles in the outgoing channel, applying exact two-body kinematics" <<
'\n');
532 p1->
boost(aBoostVector);
540 p2->adjustEnergyFromMomentum();
542 p1->
boost(-aBoostVector);
543 p2->boost(-aBoostVector);
547 INCL_DEBUG(
"Trying to adjust final-state momenta to achieve energy and momentum conservation" <<
'\n');
549 const G4int maxIterations=8;
551 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
553 std::vector<ThreeVector> minMomenta;
556 minMomenta.reserve(outgoing.size());
559 totalMomentum.
setX(0.0);
560 totalMomentum.
setY(0.0);
561 totalMomentum.
setZ(0.0);
562 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
563 totalMomentum += (*i)->getMomentum();
567 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
568 totalEnergy += (*i)->getEnergy();
571 for(
G4int iterations=0; iterations < maxIterations; ++iterations) {
574 oldOldOldVal = oldOldVal;
578 if(iterations%2 == 0) {
583 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
584 pOldTot += (*i)->getMomentum().
mag();
585 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
587 (*i)->setMomentum(mom + deltaP*mom.
mag()/pOldTot);
588 (*i)->adjustEnergyFromMomentum();
594 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
596 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.
mag2();
598 (*i)->setEnergy((*i)->getEnergy()*energyScale);
599 (*i)->adjustMomentumFromEnergy();
605 totalMomentum.
setX(0.0);
606 totalMomentum.
setY(0.0);
607 totalMomentum.
setZ(0.0);
609 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
610 totalMomentum += (*i)->getMomentum();
611 totalEnergy += (*i)->getEnergy();
617 INCL_DEBUG(
"Merit function: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal <<
'\n');
621 INCL_DEBUG(
"New minimum found, storing the particle momenta" <<
'\n');
623 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
624 minMomenta.push_back((*i)->getMomentum());
628 if(val > oldOldVal && oldVal > oldOldOldVal) {
629 INCL_DEBUG(
"Search is diverging, breaking out of the iteration loop: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal <<
'\n');
639 std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
640 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
641 (*i)->setMomentum(*v);
642 (*i)->adjustEnergyFromMomentum();
652 G4bool isNucleonAbsorption =
false;
654 G4bool isPionAbsorption =
false;
660 isPionAbsorption =
true;
671 if(outgoingParticles.size() == 0 &&
674 isNucleonAbsorption =
true;
681 for(
ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
684 if(isPionAbsorption) {
686 isPionAbsorption =
false;
690 eventInfo->
A[eventInfo->
nParticles] = (*i)->getA();
691 eventInfo->
Z[eventInfo->
nParticles] = (*i)->getZ();
693 eventInfo->
EKin[eventInfo->
nParticles] = (*i)->getKineticEnergy();
701 eventInfo->
history.push_back(
"");
713 if(std::abs(eStar)<1E-10)
717 INCL_WARN(
"Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->
EStarRem[eventInfo->
nRemnants] <<
'\n');
744 INCL_WARN(
"Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->
EStarRem[eventInfo->
nRemnants] <<
'\n');
785 theBalance.
Z = theEventInfo.
Zp + theEventInfo.
Zt;
786 theBalance.
A = theEventInfo.
Ap + theEventInfo.
At;
793 for(
ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
794 theBalance.
Z -= (*i)->getZ();
795 theBalance.
A -= (*i)->getA();
798 theBalance.
energy -= (*i)->getEnergy();
799 theBalance.
momentum -= (*i)->getMomentum();
814 theBalance.
Z -=
getZ();
815 theBalance.
A -=
getA();
843 const G4double anExcitationEnergy = aMass
Short_t Zp
Charge number of the projectile nucleus.
G4int getA() const
Returns the baryon number.
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
void initializeParticles()
Call the Cluster method to generate the initial distribution of particles.
Short_t nRemnants
Number of remnants.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
The INCL configuration object.
void setMass(G4double mass)
Set the mass of the particle in MeV/c^2.
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
G4double getFirstCollisionSpectatorMomentum() const
ThreeVector incomingMomentum
FinalStateValidity getValidity() const
virtual G4double getTableMass() const
Get the real particle mass.
G4double getFirstCollisionSpectatorPosition() const
Simple class implementing De Jong's spin model for nucleus-nucleus collisions.
G4double getMass() const
Get the cached particle mass.
ThreeVector incomingAngularMomentum
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
ParticleList const & getParticles() const
Return the list of "active" particles (i.e.
Float_t firstCollisionTime
Time of the first collision [fm/c].
std::string print()
Print the nucleus info.
G4int getEmittedClusters() const
G4int getAcceptedCollisions() const
Float_t emissionTime[maxSizeParticles]
Emission time [fm/c].
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
ParticleList const & getParticles() const
Get the list of particles in the cluster.
G4double getEmissionTime()
Int_t nCollisions
Number of accepted two-body collisions.
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
Float_t firstCollisionXSec
Cross section of the first collision (mb)
ParticleList const & getModifiedParticles() const
void setEmissionTime(G4double t)
Static class for carrying out cluster decays.
void computeOneNucleonRecoilKinematics()
Compute the recoil kinematics for a 1-nucleon remnant.
void applyFinalState(FinalState *)
Apply reaction final state information to the nucleus.
const G4double hc
[MeV*fm]
void boost(const ThreeVector &aBoostVector)
Boost the particle using a boost vector.
const G4INCL::ThreeVector & getMomentum() const
Get the momentum vector.
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
Bool_t forcedCompoundNucleus
True if the event is a forced CN.
G4double initialInternalEnergy
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
G4bool isDelta() const
Is it a Delta?
std::string print() const
void add(Particle *p)
Add one particle to the store.
Int_t nBlockedCollisions
Number of two-body collisions blocked by Pauli or CDPP.
Struct for conservation laws.
Bool_t pionAbsorption
True if the event is a pion absorption.
G4double toDegrees(G4double radians)
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
Float_t JRem[maxSizeRemnants]
Remnant spin [ ].
Int_t nReflectionAvatars
Number of reflection avatars.
Cluster is a particle (inherits from the Particle class) that is actually a collection of elementary ...
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
G4int getBlockedCollisions() const
G4double mag2() const
Get the square of the length.
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
G4double getFirstCollisionTime() const
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
G4double getPotentialEnergy() const
Get the particle potential energy.
Int_t nCollisionAvatars
Number of collision avatars.
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
G4double getInvariantMass() const
Get the the particle invariant mass.
Short_t At
Mass number of the target nucleus.
void setY(G4double ay)
Set the y coordinate.
void setEnergy(G4double energy)
Set the energy of the particle in MeV.
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
ParticleList const & getOutgoingParticles() const
Short_t ARem[maxSizeRemnants]
Remnant mass number.
ParticleList const & getCreatedParticles() const
Short_t Zt
Charge number of the target nucleus.
G4int getAvatars(AvatarType type) const
void removeScheduledAvatars()
Remove avatars that have been scheduled.
Float_t theta[maxSizeParticles]
Particle momentum polar angle [radians].
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
G4double getInitialEnergy() const
Get the initial energy.
Final state of an interaction.
G4double getSeparationEnergy(const Particle *const p) const
Return the separation energy for a particle.
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
ParticleList const & getOutgoingParticles() const
Return the list of outgoing particles (i.e.
Short_t Z[maxSizeParticles]
Particle charge number.
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
Book & getBook()
Return the pointer to the Book object which keeps track of various counters.
Short_t nParticles
Number of particles in the final state.
G4bool nucleon(G4int ityp)
G4int getBlockedDecays() const
void particleHasBeenDestroyed(Particle *const)
Remove the particle from the system.
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
virtual void initializeParticles()
Initialise the NuclearDensity pointer and sample the particles.
G4double phi() const
Phi angle.
void particleHasBeenEjected(Particle *const)
Mark the particle as ejected.
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
G4int getZ() const
Returns the charge number.
G4bool getPionPotential() const
Do we want the pion potential?
NuclearDensity const * theDensity
Pointer to the NuclearDensity object.
void setPotentialEnergy(G4double v)
Set the particle potential energy.
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
G4double computeSeparationEnergyBalance() const
Outgoing - incoming separation energies.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
NuclearDensity const * createDensity(const G4int A, const G4int Z)
G4double computeExcitationEnergy() const
Compute the current excitation energy.
Int_t nDecayAvatars
Number of decay avatars.
void deleteProjectileRemnant()
Delete the projectile remnant.
G4INCL::ThreeVector thePosition
G4bool decayMe()
Force the phase-space decay of the Nucleus.
G4int getAcceptedDecays() const
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
ParticleSampler * theParticleSampler
Float_t firstCollisionSpectatorPosition
Position of the spectator on the first collision (fm)
Float_t phiRem[maxSizeRemnants]
Remnant momentum azimuthal angle [radians].
G4bool getFirstCollisionIsElastic() const
Int_t nBlockedDecays
Number of decays blocked by Pauli or CDPP.
G4bool hasPionPotential() const
Do we have a pion potential?
ThreeVector computeCenterOfMass() const
Compute the current center-of-mass position.
G4double getTotalEnergyBeforeInteraction() const
G4double getFirstCollisionXSec() const
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Short_t A[maxSizeParticles]
Particle mass number.
Int_t nDecays
Number of accepted Delta decays.
Bool_t firstCollisionIsElastic
True if the first collision was elastic.
G4double theExcitationEnergy
void particleHasBeenUpdated(Particle *const)
Notify the Store about a particle update.
void propagateParticles(G4double step)
Propagate the particles one time step.
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
PotentialType getPotentialType() const
Get the type of the potential for nucleons.
NuclearPotential::INuclearPotential const * thePotential
Pointer to the NuclearPotential object.
std::string print() const
void emitInsidePions()
Force emission of all pions inside the nucleus.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
void setTableMass()
Set the mass of the Particle to its table mass.
ParticleList const & getEnteringParticles() const
ParticleList const & getDestroyedParticles() const
void setPotential(NuclearPotential::INuclearPotential const *const p)
Setter for thePotential.
Float_t thetaRem[maxSizeRemnants]
Remnant momentum polar angle [radians].
G4bool isEventTransparent() const
Is the event transparent?
FinalState * getFinalState()
void fillEventInfo(EventInfo *eventInfo)
Fill the event info which contains INCL output data.
G4double computeTotalEnergy() const
Compute the current total energy.
G4INCL::ThreeVector theMomentum
ThreeVector const & getSpin() const
Get the spin of the nucleus.
virtual G4double getTableMass() const
Get the tabulated particle mass.
void setX(G4double ax)
Set the x coordinate.
G4double getExcitationEnergy() const
Get the excitation energy of the cluster.
G4double getKineticEnergy() const
Get the particle kinetic energy.
The purpose of the Store object is to act as a "particle manager" that keeps track ofall the particle...
void insertParticle(Particle *p)
Insert a new particle (e.g. a projectile) in the nucleus.
G4double theta() const
Theta angle.
Float_t phi[maxSizeParticles]
Particle momentum azimuthal angle [radians].
G4double mag() const
Get the length of the vector.
G4double getCurrentTime() const
Nucleus(G4int mass, G4int charge, Config const *const conf, const G4double universeRadius=-1.)
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
void setZ(G4double az)
Set the z coordinate.
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
INuclearPotential const * createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential)
Create an INuclearPotential object.
Float_t firstCollisionSpectatorMomentum
Momentum of the spectator on the first collision (fm)
ThreeVector initialCenterOfMass
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
G4int getEnergyViolationInteraction() const
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
Short_t Ap
Mass number of the projectile nucleus.
Short_t nCascadeParticles
Number of cascade particles.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
const G4double effectiveNucleonMass
ProjectileRemnant * theProjectileRemnant
Pointer to the quasi-projectile.
Short_t origin[maxSizeParticles]
Origin of the particle.
#define INCL_DATABLOCK(x)
void setNeutronSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
ParticleList::const_iterator ParticleIter
Int_t projectileType
Projectile particle type.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Set the momentum vector.
G4double theUniverseRadius
The radius of the universe.
std::vector< std::string > history
History of the particle.