33 #define INCLXX_IN_GEANT4_MODE 1
44 #ifndef G4INCLNucleus_hh
45 #define G4INCLNucleus_hh 1
66 theInitialZ(charge), theInitialA(mass),
67 theNpInitial(0), theNnInitial(0),
68 initialInternalEnergy(0.),
69 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
70 initialCenterOfMass(0.,0.,0.),
77 theUniverseRadius(universeRadius),
78 isNucleusNucleus(false),
79 theProjectileRemnant(NULL),
148 for(
ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
150 if(!(*iter)->isOutOfWell()) {
151 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
157 for(
ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
162 for(
ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
164 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
169 for(
ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
170 if((*iter)->isCluster()) {
173 #ifdef INCLXX_IN_GEANT4_MODE
178 for(
ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
183 totalEnergy += (*iter)->getEnergy();
184 theA -= (*iter)->getA();
185 theZ -= (*iter)->getZ();
191 for(
ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
193 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
199 INCL_DEBUG(
"A Particle is entering below the Fermi sea:" << std::endl << finalstate->
print() << std::endl);
202 for(
ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
209 INCL_ERROR(
"Energy nonconservation! Energy at the beginning of the event = "
211 <<
" and after interaction = "
212 << totalEnergy << std::endl
213 << finalstate->
print());
218 INCL_WARN(
"Useless Nucleus::propagateParticles -method called." << std::endl);
224 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
225 if((*p)->isNucleon())
226 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
227 else if((*p)->isResonance())
230 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
250 for(
ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p) {
252 theSpin -= (*p)->getAngularMomentum();
272 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
273 const G4double mass = (*p)->getMass();
274 cm += (*p)->getPosition() * mass;
290 std::stringstream ss;
291 ss <<
"Particles in the nucleus:" << std::endl
292 <<
"Inside:" << std::endl;
295 for(
ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
296 ss <<
"index = " << counter << std::endl
300 ss <<
"Outgoing:" << std::endl;
302 for(
ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p)
311 for(
ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
312 if((*i)->isDelta()) deltas.push_back((*i));
314 if(deltas.empty())
return false;
316 for(
ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
317 INCL_DEBUG(
"Decay outgoing delta particle:" << std::endl
318 << (*i)->print() << std::endl);
320 const G4double deltaMass = (*i)->getMass();
325 (*i)->setEnergy((*i)->getMass());
338 newMomentum *= decayMomentum / newMomentum.
mag();
350 nucleon->
boost(beta);
369 const G4bool unphysicalRemnant = (theZ<0 || theZ>
theA);
376 for(
ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
377 if((*i)->isDelta()) deltas.push_back((*i));
380 for(
ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
381 INCL_DEBUG(
"Decay inside delta particle:" << std::endl
382 << (*i)->print() << std::endl);
388 if(unphysicalRemnant) {
389 INCL_WARN(
"Forcing delta decay inside an unphysical remnant (A=" <<
theA
390 <<
", Z=" <<
theZ <<
"). Might lead to energy-violation warnings."
409 if(unphysicalRemnant) {
410 INCL_DEBUG(
"Remnant is unphysical: Z=" <<
theZ <<
", A=" <<
theA <<
", emitting all the pions" << std::endl);
420 for(
ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
421 if((*i)->isCluster()) clusters.push_back((*i));
423 if(clusters.empty())
return false;
425 for(
ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
428 #ifdef INCLXX_IN_GEANT4_MODE
434 for(
ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j)
446 for(
ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j)
457 INCL_WARN(
"Forcing emissions of all pions in the nucleus." << std::endl);
460 const G4double tinyPionEnergy = 0.1;
465 for(
ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
468 INCL_DEBUG(
"Forcing emission of the following particle: "
469 << thePion->
print() << std::endl);
475 if(kineticEnergyOutside > 0.0)
482 toEject.push_back(thePion);
485 for(
ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
497 if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
508 INCL_ERROR(
"Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << std::endl);
527 for(
ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
535 if(outgoing.size() == 2) {
537 INCL_DEBUG(
"Two particles in the outgoing channel, applying exact two-body kinematics" << std::endl);
541 Particle *p1 = outgoing.front(), *p2 = outgoing.back();
544 p1->
boost(aBoostVector);
552 p2->adjustEnergyFromMomentum();
554 p1->
boost(-aBoostVector);
555 p2->boost(-aBoostVector);
559 INCL_DEBUG(
"Trying to adjust final-state momenta to achieve energy and momentum conservation" << std::endl);
561 const G4int maxIterations=8;
563 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
565 std::vector<ThreeVector> minMomenta;
568 minMomenta.reserve(outgoing.size());
571 totalMomentum.
setX(0.0);
572 totalMomentum.
setY(0.0);
573 totalMomentum.
setZ(0.0);
574 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
575 totalMomentum += (*i)->getMomentum();
579 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
580 totalEnergy += (*i)->getEnergy();
583 for(
G4int iterations=0; iterations < maxIterations; ++iterations) {
586 oldOldOldVal = oldOldVal;
590 if(iterations%2 == 0) {
595 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
596 pOldTot += (*i)->getMomentum().
mag();
597 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
599 (*i)->setMomentum(mom + deltaP*mom.
mag()/pOldTot);
600 (*i)->adjustEnergyFromMomentum();
606 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
608 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.
mag2();
610 (*i)->setEnergy((*i)->getEnergy()*energyScale);
611 (*i)->adjustMomentumFromEnergy();
617 totalMomentum.
setX(0.0);
618 totalMomentum.
setY(0.0);
619 totalMomentum.
setZ(0.0);
621 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
622 totalMomentum += (*i)->getMomentum();
623 totalEnergy += (*i)->getEnergy();
629 INCL_DEBUG(
"Merit function: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal << std::endl);
633 INCL_DEBUG(
"New minimum found, storing the particle momenta" << std::endl);
635 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
636 minMomenta.push_back((*i)->getMomentum());
640 if(val > oldOldVal && oldVal > oldOldOldVal) {
641 INCL_DEBUG(
"Search is diverging, breaking out of the iteration loop: val=" << val <<
", oldVal=" << oldVal <<
", oldOldVal=" << oldOldVal <<
", oldOldOldVal=" << oldOldOldVal << std::endl);
650 INCL_DEBUG(
"Applying the solution" << std::endl);
651 std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
652 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
653 (*i)->setMomentum(*v);
654 (*i)->adjustEnergyFromMomentum();
664 G4bool isNucleonAbsorption =
false;
666 G4bool isPionAbsorption =
false;
672 isPionAbsorption =
true;
683 if(outgoingParticles.size() == 0 &&
686 isNucleonAbsorption =
true;
693 for(
ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
696 if(isPionAbsorption) {
698 isPionAbsorption =
false;
702 eventInfo->
A[eventInfo->
nParticles] = (*i)->getA();
703 eventInfo->
Z[eventInfo->
nParticles] = (*i)->getZ();
705 eventInfo->
EKin[eventInfo->
nParticles] = (*i)->getKineticEnergy();
713 eventInfo->
history.push_back(
"");
725 if(std::abs(eStar)<1E-10)
729 INCL_WARN(
"Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->
EStarRem[eventInfo->
nRemnants] << std::endl);
756 INCL_WARN(
"Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->
EStarRem[eventInfo->
nRemnants] << std::endl);
797 theBalance.
Z = theEventInfo.
Zp + theEventInfo.
Zt;
798 theBalance.
A = theEventInfo.
Ap + theEventInfo.
At;
805 for(
ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
806 theBalance.
Z -= (*i)->getZ();
807 theBalance.
A -= (*i)->getA();
810 theBalance.
energy -= (*i)->getEnergy();
811 theBalance.
momentum -= (*i)->getMomentum();
826 theBalance.
Z -=
getZ();
827 theBalance.
A -=
getA();
855 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.
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.
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
Particle * getBlockedDelta()
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
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 getCurrentTime()
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
Int_t nRemnants
Number of remnants.
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 that a particle has been updated.
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?
G4INCL::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.
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.