34 #define INCLXX_IN_GEANT4_MODE 1
79 :propagationModel(0), theA(208), theZ(82),
80 targetInitSuccess(false),
82 maxUniverseRadius(0.),
83 maxInteractionDistance(0.),
84 fixedImpactParameter(0.),
87 forceTransparent(false),
91 #ifdef INCLXX_IN_GEANT4_MODE
93 #else // INCLXX_IN_GEANT4_MODE
95 #endif // INCLXX_IN_GEANT4_MODE
143 #ifndef INCLXX_IN_GEANT4_MODE
158 #ifndef INCLXX_IN_GEANT4_MODE
159 NuclearMassTable::deleteTable();
167 #ifndef INCLXX_IN_GEANT4_MODE
168 Logger::deleteLoggerSlave();
179 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
180 INCL_ERROR(
"Unsupported target: A = " << A <<
" Z = " << Z <<
'\n'
181 <<
"Target configuration rejected." <<
'\n');
185 (projectileSpecies.
theZ==projectileSpecies.
theA || projectileSpecies.
theZ==0)) {
186 INCL_ERROR(
"Unsupported projectile: A = " << projectileSpecies.
theA <<
" Z = " << projectileSpecies.
theZ <<
'\n'
187 <<
"Projectile configuration rejected." <<
'\n');
217 if(projectileSpecies.
theA > 0)
246 INCL_WARN(
"Target initialisation failed for A=" << targetA <<
", Z=" << targetZ <<
'\n');
295 INCL_DEBUG(
"Selected impact parameter: " << impactParameter <<
'\n');
301 if(effectiveImpactParameter < 0.) {
331 if(avatar == 0)
break;
364 INCL_DEBUG(
"Trying compound nucleus" <<
'\n');
368 #ifndef INCLXX_IN_GEANT4_MODE
378 if(projectileRemnant) {
406 INCL_DEBUG(
"Cascade resulted in complete fusion, using realistic fusion kinematics" <<
'\n');
412 INCL_WARN(
"Complete-fusion kinematics yields negative excitation energy, returning a transparent!" <<
'\n');
428 INCL_ERROR(
"Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." <<
'\n');
432 #ifndef INCLXX_IN_GEANT4_MODE
434 globalConservationChecks(
false);
446 #ifndef INCLXX_IN_GEANT4_MODE
448 globalConservationChecks(
true);
487 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
489 std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(),
Random::getAdapter());
492 G4bool atLeastOneNucleonEntering =
false;
493 for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
497 (*p)->getPropagationVelocity(),
499 if(!intersectionInteractionDistance.
exists)
503 atLeastOneNucleonEntering =
true;
516 theCNZ += (*p)->getZ();
526 if(!success || !atLeastOneNucleonEntering) {
527 INCL_DEBUG(
"No nucleon entering in forced CN, forcing a transparent" <<
'\n');
537 theCNEnergy -= theProjectileRemnant->
getEnergy();
538 theCNMomentum -= theProjectileRemnant->
getMomentum();
549 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.
mag2();
550 if(theCNInvariantMassSquared<0.) {
555 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
556 if(theCNExcitationEnergy<0.) {
558 INCL_DEBUG(
"CN excitation energy is negative, forcing a transparent" <<
'\n'
559 <<
" theCNA = " << theCNA <<
'\n'
560 <<
" theCNZ = " << theCNZ <<
'\n'
561 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
562 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
563 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
564 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
570 INCL_DEBUG(
"CN excitation energy is positive, forcing a CN" <<
'\n'
571 <<
" theCNA = " << theCNA <<
'\n'
572 <<
" theCNZ = " << theCNZ <<
'\n'
573 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
574 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
575 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
576 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
603 theRecoilFunctor(theSolution.
x);
605 INCL_WARN(
"Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." <<
'\n');
610 #ifndef INCLXX_IN_GEANT4_MODE
611 void INCL::globalConservationChecks(
G4bool afterRecoil) {
617 if(theBalance.
Z != 0) {
618 INCL_ERROR(
"Violation of charge conservation! ZBalance = " << theBalance.
Z <<
'\n');
620 if(theBalance.
A != 0) {
621 INCL_ERROR(
"Violation of baryon-number conservation! ABalance = " << theBalance.
A <<
'\n');
623 G4double EThreshold, pLongThreshold, pTransThreshold;
628 pTransThreshold = 1.;
632 pLongThreshold = 0.1;
633 pTransThreshold = 0.1;
635 if(std::abs(theBalance.
energy)>EThreshold) {
636 INCL_WARN(
"Violation of energy conservation > " << EThreshold <<
" MeV. EBalance = " << theBalance.
energy <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
638 if(std::abs(pLongBalance)>pLongThreshold) {
639 INCL_WARN(
"Violation of longitudinal momentum conservation > " << pLongThreshold <<
" MeV/c. pLongBalance = " << pLongBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
641 if(std::abs(pTransBalance)>pTransThreshold) {
642 INCL_WARN(
"Violation of transverse momentum conservation > " << pTransThreshold <<
" MeV/c. pTransBalance = " << pTransBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
657 <<
"), stopping cascade" <<
'\n');
663 INCL_DEBUG(
"No participants in the nucleus and no incoming particles left, stopping cascade" <<
'\n');
670 <<
"), stopping cascade" <<
'\n');
676 INCL_DEBUG(
"Trying to make a compound nucleus, stopping cascade" <<
'\n');
720 G4int nUnmergedSpectators = 0;
723 if(dynSpectators.empty() && geomSpectators.empty()) {
725 }
else if(dynSpectators.size()==1 && geomSpectators.empty()) {
736 nUnmergedSpectators = rejected.size();
744 return nUnmergedSpectators;
758 INCL_DEBUG(
"Initialised interaction distance: r0 = " << r0 <<
'\n'
759 <<
" theNNDistance = " << theNNDistance <<
'\n'
769 for(
IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
773 rMax =
std::max(maximumRadius, rMax);
778 rMax =
std::min(pMaximumRadius, nMaximumRadius);
Short_t Ap
Projectile mass number given as input.
Short_t Zp
Charge number of the projectile nucleus.
G4int getA() const
Returns the baryon number.
Float_t effectiveImpactParameter
Effective (Coulomb-distorted) impact parameter [fm].
G4double maxUniverseRadius
void fillFinalState(FinalState *fs)
CascadeActionType getCascadeActionType() const
Get the cascade-action type.
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.
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
G4int minRemnantSize
Remnant size below which cascade stops.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
IsotopeVector::iterator IsotopeIter
The INCL configuration object.
void setMass(G4double mass)
Set the mass of the particle in MeV/c^2.
static void setCutNN(const G4double c)
void deleteCrossSections()
FinalStateValidity getValidity() const
void deleteGenerator()
Delete the generator.
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z)
G4bool continueCascade()
Stopping criterion for the cascade.
void cascade()
The actual cascade loop.
Float_t pLongBalance
Longitudinal momentum-conservation balance [MeV/c].
G4double maxInteractionDistance
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
void postCascade()
Finalise the cascade and clean up.
G4int getTargetZ() const
Get the target charge number.
G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy)
Initialise the cascade.
Int_t nForcedCompoundNucleus
Number of forced compound-nucleus events.
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
std::string cascadeModel
Name of the cascade model.
ParticleList const & getParticles() const
Get the list of particles in the cluster.
void beforeRunAction(Config const *config)
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
Abstract interface to the nuclear potential.
virtual G4INCL::IAvatar * propagate(FinalState const *const fs)=0
Propagate the particles and get the next avatar.
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
void beforeCascadeAction(IPropagationModel *)
Config const *const theConfig
void applyFinalState(FinalState *)
Apply reaction final state information to the nucleus.
void updateGlobalInfo()
Update global counters and other members of theGlobalInfo object.
void setZ(const G4int Z)
Set the charge number of the cluster.
const G4INCL::ThreeVector & getMomentum() const
Get the momentum vector.
void initialize(Config const *const)
Initialize generator according to a Config object.
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
G4double interactionDistanceNN(const ParticleSpecies &aSpecies, const G4double kineticEnergy)
Compute the "interaction distance".
Float_t forcedCNCrossSection
Calculated forced-compound-nucleus cross section.
Float_t pionAbsorptionCrossSection
Pion absorption cross section.
std::vector< Int_t > finalRandomSeeds
Final seeds for the pseudo-random-number generator.
Bool_t forcedDeltasOutside
Event involved forced delta decays outside the nucleus.
G4double getEnergy() const
Get the energy of the particle in MeV.
void deletePhaseSpaceGenerator()
Struct for conservation laws.
Bool_t pionAbsorption
True if the event is a pion absorption.
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
G4bool isNaturalTarget() const
Natural targets.
Bool_t clusterDecay
Event involved cluster decay.
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
void deleteIncoming()
Clear the incoming list and delete the particles.
ParticleList const & getIncomingParticles() const
Return the list of incoming particles (i.e.
Float_t Ep
Projectile kinetic energy given as input.
Short_t At
Target mass number given as input.
Cluster is a particle (inherits from the Particle class) that is actually a collection of elementary ...
Int_t nTransparents
Number of transparent shots.
G4double shoot0()
Return a random number in the ]0,1] interval.
void reset()
Reset the projectile remnant to the state at the beginning of the cascade.
std::vector< Int_t > initialRandomSeeds
Initial seeds for the pseudo-random-number generator.
CascadeAction * cascadeAction
Float_t stoppingTime
Cascade stopping time [fm/c].
G4double mag2() const
Get the square of the length.
Class containing default actions to be performed at intermediate cascade steps.
Float_t Ep
Projectile kinetic energy given as input.
IPropagationModel * propagationModel
Short_t At
Mass number of the target nucleus.
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
void setEnergy(G4double energy)
Set the energy of the particle in MeV.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
static G4ThreadLocal Int_t eventNumber
Number of the event.
SeedVector getSeeds()
Get the seeds of the current generator.
void clearCache()
Clear the INuclearPotential cache.
void beforeAvatarAction(IAvatar *a, Nucleus *n)
G4double getImpactParameter() const
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
Short_t Zt
Charge number of the target nucleus.
Float_t impactParameter
Impact parameter [fm].
G4double fixedImpactParameter
void rescaleOutgoingForRecoil()
Rescale the energies of the outgoing particles.
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
void deleteCoulomb()
Delete the Coulomb-distortion object.
Static class for selecting Coulomb distortion.
INCL(Config const *const config)
void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy)
Initialise the maximum interaction distance.
static std::string const getVersionString()
Get the INCL version string.
Final state of an interaction.
Simple container for output of calculation-wide results.
Float_t pTransBalance
Transverse momentum-conservation balance [MeV/c].
void deleteBlockers()
Delete blockers.
Float_t EBalance
Energy-conservation balance [MeV].
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
ParticleList const & getOutgoingParticles() const
Return the list of outgoing particles (i.e.
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
G4int getCascading() const
virtual G4double shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)=0
Book & getBook()
Return the pointer to the Book object which keeps track of various counters.
Bool_t deltasInside
Event involved deltas in the nucleus at the end of the cascade.
ParticleList extractDynamicalSpectators()
Returns a list of dynamical spectators.
Int_t nForcedTransparents
Number of forced transparents.
Class that stores isotopic abundances for a given element.
std::string deexcitationModel
Name of the de-excitation model.
void beforePropagationAction(IPropagationModel *pm)
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z)
Initialize the universe radius.
void initialize(Config const *const theConfig)
G4int getTargetA() const
Get the target mass number.
Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r)
Compute the first intersection of a straight particle trajectory with a sphere.
const EventInfo & processEvent()
G4int getZ() const
Returns the charge number.
Alternative CascadeAction for dumping avatars.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
Class to adjust remnant recoil in the reaction CM system.
G4double computeExcitationEnergy() const
Compute the current excitation energy.
Float_t completeFusionCrossSection
Calculated complete-fusion cross section (nParticles==0)
Int_t nShots
Number of shots.
Int_t nUnmergedSpectators
Number of dynamical spectators that were merged back into the projectile remnant. ...
void afterCascadeAction(Nucleus *)
static const G4double A[nN]
G4bool decayMe()
Force the phase-space decay of the Nucleus.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
void clearOutgoing()
Clear all outgoing particles from the store.
std::vector< Isotope > IsotopeVector
IsotopeVector const & getIsotopes() const
Get the isotope vector.
Functions that encapsulate a mass table.
Float_t reactionCrossSection
Calculated reaction cross section.
Standard INCL4 particle propagation and avatar prediction.
G4bool getTryCompoundNucleus()
Float_t geometricCrossSection
Geometric cross section.
G4bool initializeTarget(const G4int A, const G4int Z)
Short_t Zp
Projectile charge number given as input.
Adapter const & getAdapter()
void initVerbosityLevelFromEnvvar()
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
void setA(const G4int A)
Set the mass number of the cluster.
ParticleList addAllDynamicalSpectators(ParticleList const &pL)
Add back all dynamical spectators to the projectile remnant.
G4double maxImpactParameter
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Static class for cluster formation.
void deleteClusteringModel()
Delete the clustering model.
std::string getDeExcitationString() const
Get the de-excitation string.
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void setNucleus(G4INCL::Nucleus *nucleus)=0
Set the nucleus for the propagation model.
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
G4double interactionDistancePiN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
void initialize(Config const *const theConfig)
void reset()
Reset the EventInfo members.
void makeCompoundNucleus()
Make a compound nucleus.
virtual G4double getStoppingTime()=0
Get the current stopping time.
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool isEventTransparent() const
Is the event transparent?
FinalState * getFinalState()
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
void fillEventInfo(EventInfo *eventInfo)
Fill the event info which contains INCL output data.
G4int drawRandomNaturalIsotope(const G4int Z)
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
G4double shoot()
Generate flat distribution of random numbers.
void initialize(Config const *const aConfig)
Initialise blockers according to a Config object.
Short_t Zt
Target charge number given as input.
Intersection-point structure.
Bool_t transparent
True if the event is transparent.
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
void afterPropagationAction(IPropagationModel *pm, IAvatar *avatar)
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
virtual G4double getCurrentTime()=0
Returns the current global time of the system.
void afterAvatarAction(IAvatar *a, Nucleus *n, FinalState *fs)
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
G4int makeProjectileRemnant()
Make a projectile pre-fragment out of geometrical spectators.
Int_t nNucleonAbsorptions
Number of nucleon absorptions (no outcoming particles)
Float_t energyViolationInteractionCrossSection
Cross section for attempted collisions/decays for which the energy-conservation algorithm failed to f...
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
Bool_t forcedDeltasInside
Event involved forced delta decays inside the nucleus.
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.
void initialize(Config const *const theConfig=0)
Initialize the particle table.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
Float_t errorForcedCNCrossSection
Error on the calculated forced-compound-nucleus cross section.
G4double getCutNN() const
void distortOut(ParticleList const &pL, Nucleus const *const n)
Modify the momentum of an outgoing particle.
Simple class for computing intersections between a straight line and a sphere.
void initialize(Config const *const theConfig)
Initialize the Coulomb-distortion algorithm.
void clearIncoming()
Clear the incoming list.
Int_t nCompleteFusion
Number of complete-fusion events (nParticles==0)
Int_t projectileType
Projectile particle type.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Set the momentum vector.