34 #define INCLXX_IN_GEANT4_MODE 1
80 :propagationModel(0), theA(208), theZ(82), theS(0),
81 targetInitSuccess(false),
83 maxUniverseRadius(0.),
84 maxInteractionDistance(0.),
85 fixedImpactParameter(0.),
88 forceTransparent(false),
92 #ifdef INCLXX_IN_GEANT4_MODE
94 #else // INCLXX_IN_GEANT4_MODE
96 #endif // INCLXX_IN_GEANT4_MODE
147 #ifndef INCLXX_IN_GEANT4_MODE
162 #ifndef INCLXX_IN_GEANT4_MODE
163 NuclearMassTable::deleteTable();
171 #ifndef INCLXX_IN_GEANT4_MODE
172 Logger::deleteLoggerSlave();
183 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
184 INCL_ERROR(
"Unsupported target: A = " << A <<
" Z = " << Z <<
" S = " << S <<
'\n'
185 <<
"Target configuration rejected." <<
'\n');
189 (projectileSpecies.
theZ==projectileSpecies.
theA || projectileSpecies.
theZ==0)) {
190 INCL_ERROR(
"Unsupported projectile: A = " << projectileSpecies.
theA <<
" Z = " << projectileSpecies.
theZ <<
" S = " << projectileSpecies.
theS <<
'\n'
191 <<
"Projectile configuration rejected." <<
'\n');
222 if(projectileSpecies.
theA > 0)
257 INCL_WARN(
"Target initialisation failed for A=" << targetA <<
", Z=" << targetZ <<
", S=" << targetS <<
'\n');
308 INCL_DEBUG(
"Selected impact parameter: " << impactParameter <<
'\n');
314 if(effectiveImpactParameter < 0.) {
331 unsigned long loopCounter = 0;
332 const unsigned long maxLoopCounter = 10000000;
346 if(avatar == 0)
break;
380 INCL_DEBUG(
"Trying compound nucleus" <<
'\n');
384 #ifndef INCLXX_IN_GEANT4_MODE
394 if(projectileRemnant) {
416 #ifdef INCLXX_IN_GEANT4_MODE
418 #endif // INCLXX_IN_GEANT4_MODE
449 INCL_DEBUG(
"Cascade resulted in complete fusion, using realistic fusion kinematics" <<
'\n');
455 INCL_WARN(
"Complete-fusion kinematics yields negative excitation energy, returning a transparent!" <<
'\n');
471 INCL_ERROR(
"Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." <<
'\n');
475 #ifndef INCLXX_IN_GEANT4_MODE
477 globalConservationChecks(
false);
489 #ifndef INCLXX_IN_GEANT4_MODE
491 globalConservationChecks(
true);
530 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
532 std::random_shuffle(shuffledComponents.begin(), shuffledComponents.end(),
Random::getAdapter());
535 G4bool atLeastOneNucleonEntering =
false;
536 for(std::vector<Particle*>::const_iterator
p=shuffledComponents.begin(),
e=shuffledComponents.end();
p!=
e; ++
p) {
540 (*p)->getPropagationVelocity(),
542 if(!intersectionInteractionDistance.
exists)
546 atLeastOneNucleonEntering =
true;
559 theCNZ += (*p)->getZ();
569 if(!success || !atLeastOneNucleonEntering) {
570 INCL_DEBUG(
"No nucleon entering in forced CN, forcing a transparent" <<
'\n');
580 theCNEnergy -= theProjectileRemnant->
getEnergy();
581 theCNMomentum -= theProjectileRemnant->
getMomentum();
592 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.
mag2();
593 if(theCNInvariantMassSquared<0.) {
598 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
599 if(theCNExcitationEnergy<0.) {
601 INCL_DEBUG(
"CN excitation energy is negative, forcing a transparent" <<
'\n'
602 <<
" theCNA = " << theCNA <<
'\n'
603 <<
" theCNZ = " << theCNZ <<
'\n'
604 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
605 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
606 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
607 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
613 INCL_DEBUG(
"CN excitation energy is positive, forcing a CN" <<
'\n'
614 <<
" theCNA = " << theCNA <<
'\n'
615 <<
" theCNZ = " << theCNZ <<
'\n'
616 <<
" theCNEnergy = " << theCNEnergy <<
'\n'
617 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" <<
'\n'
618 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy <<
'\n'
619 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" <<
'\n'
650 theRecoilFunctor(theSolution.
x);
652 INCL_WARN(
"Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." <<
'\n');
657 #ifndef INCLXX_IN_GEANT4_MODE
658 void INCL::globalConservationChecks(
G4bool afterRecoil) {
664 if(theBalance.
Z != 0) {
665 INCL_ERROR(
"Violation of charge conservation! ZBalance = " << theBalance.
Z <<
'\n');
667 if(theBalance.
A != 0) {
670 G4double EThreshold, pLongThreshold, pTransThreshold;
675 pTransThreshold = 1.;
679 pLongThreshold = 0.1;
680 pTransThreshold = 0.1;
682 if(std::abs(theBalance.
energy)>EThreshold) {
685 if(std::abs(pLongBalance)>pLongThreshold) {
686 INCL_WARN(
"Violation of longitudinal momentum conservation > " << pLongThreshold <<
" MeV/c. pLongBalance = " << pLongBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
688 if(std::abs(pTransBalance)>pTransThreshold) {
689 INCL_WARN(
"Violation of transverse momentum conservation > " << pTransThreshold <<
" MeV/c. pTransBalance = " << pTransBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber <<
'\n');
704 <<
"), stopping cascade" <<
'\n');
710 INCL_DEBUG(
"No participants in the nucleus and no incoming particles left, stopping cascade" <<
'\n');
717 <<
"), stopping cascade" <<
'\n');
723 INCL_DEBUG(
"Trying to make a compound nucleus, stopping cascade" <<
'\n');
767 G4int nUnmergedSpectators = 0;
770 if(dynSpectators.empty() && geomSpectators.empty()) {
772 }
else if(dynSpectators.size()==1 && geomSpectators.empty()) {
783 nUnmergedSpectators = rejected.size();
791 return nUnmergedSpectators;
805 INCL_DEBUG(
"Initialised interaction distance: r0 = " << r0 <<
'\n'
806 <<
" theNNDistance = " << theNNDistance <<
'\n'
816 for(
IsotopeIter i=theIsotopes.begin(),
e=theIsotopes.end(); i!=
e; ++i) {
820 rMax =
std::max(maximumRadius, rMax);
826 rMax =
std::max(maximumRadius, rMax);
void afterCascadeAction(Nucleus *)
void finalizeGlobalInfo(Random::SeedVector const &initialSeeds)
virtual G4INCL::IAvatar * propagate(FinalState const *const fs)=0
Float_t effectiveImpactParameter
Effective (Coulomb-distorted) impact parameter [fm].
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual void setNucleus(G4INCL::Nucleus *nucleus)=0
void setZ(const G4int Z)
Set the charge number of the cluster.
Intersection-point structure.
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
void beforeAvatarAction(IAvatar *a, Nucleus *n)
G4double interactionDistanceKN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
G4double getProjectileKineticEnergy() const
Get the projectile kinetic energy.
G4double computeExcitationEnergy() const
Compute the current excitation energy.
G4int emitInsideLambda()
Force emission of all Lambda (desexitation code with strangeness not implanted yet) ...
G4int minRemnantSize
Remnant size below which cascade stops.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
G4bool preCascade(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy)
Initialise the cascade.
G4double getEnergy() const
void deleteClusteringModel()
Delete the clustering model.
CascadeActionType getCascadeActionType() const
Get the cascade-action type.
static std::vector< G4double > INCLBiasVector
Time ordered vector of all bias applied.
static void setCutNN(const G4double c)
Int_t nTransparents
Number of transparent shots.
G4bool containsKaon()
Returns true if the nucleus contains any Kaons.
Bool_t kaonsInside
Event involved kaons in the nucleus at the end of the cascade.
virtual G4double shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)=0
Short_t St
Strangeness number of the target nucleus.
void initialize(Config const *const)
Initialize generator according to a Config object.
CascadeAction * cascadeAction
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
void setMass(G4double mass)
INCL(Config const *const config)
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z)
Initialize the universe radius.
static std::string const getVersionString()
Get the INCL version string.
const G4INCL::ThreeVector & getMomentum() const
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
G4int getS() const
Returns the strangeness number.
Short_t At
Target mass number given as input.
G4bool decayInsideStrangeParticles()
Force the transformation of strange particles into a Lambda;.
void makeCompoundNucleus()
Make a compound nucleus.
ParticleList addAllDynamicalSpectators(ParticleList const &pL)
Add back all dynamical spectators to the projectile remnant.
G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S)
IsotopeVector const & getIsotopes() const
Get the isotope vector.
FinalStateValidity getValidity() const
ProjectileRemnant * getProjectileRemnant() const
Get the projectile remnant.
Float_t errorReactionCrossSection
Error on the calculated reaction cross section.
G4double getDecayTimeThreshold() const
Get the decay time threshold time.
Static class for cluster formation.
void reset()
Reset the EventInfo members.
static void deleteBackupParticles()
Release the memory allocated for the backup particles.
Short_t Zt
Target charge number given as input.
G4double fixedImpactParameter
G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
Float_t reactionCrossSection
Calculated reaction cross section.
G4double interactionDistanceYN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
void initialize(Config const *const aConfig)
Initialise blockers according to a Config object.
Bool_t sigmasInside
Event involved sigmas in the nucleus at the end of the cascade.
Functions that encapsulate a mass table.
void applyFinalState(FinalState *)
G4int getTargetZ() const
Get the target charge number.
Float_t Ep
Projectile kinetic energy given as input.
Float_t completeFusionCrossSection
Calculated complete-fusion cross section (nParticles==0)
void rescaleOutgoingForRecoil()
Rescale the energies of the outgoing particles.
G4double getImpactParameter() const
G4double maxInteractionDistance
G4int drawRandomNaturalIsotope(const G4int Z)
FinalState * getFinalState()
std::string getDeExcitationString() const
Get the de-excitation string.
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
G4double getHadronizationTime() const
Get the hadronization time.
void deleteCrossSections()
Int_t nForcedTransparents
Number of forced transparents.
ParticleList const & getIncomingParticles() const
void initVerbosityLevelFromEnvvar()
Alternative CascadeAction for dumping avatars.
void fillFinalState(FinalState *fs)
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
void initializeParticles()
G4double interactionDistanceNN(const ParticleSpecies &aSpecies, const G4double kineticEnergy)
Compute the "interaction distance".
Float_t pLongBalance
Longitudinal momentum-conservation balance [MeV/c].
G4int getTargetA() const
Get the target mass number.
static G4double getTotalBias()
General bias vector function.
void beforePropagationAction(IPropagationModel *pm)
Int_t nCompleteFusion
Number of complete-fusion events (nParticles==0)
void fillEventInfo(EventInfo *eventInfo)
Config const *const theConfig
Int_t nUnmergedSpectators
Number of dynamical spectators that were merged back into the projectile remnant. ...
Adapter const & getAdapter()
ParticleSpecies getProjectileSpecies() const
Get the projectile species.
ParticleList const & getParticles() const
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
void deleteCoulomb()
Delete the Coulomb-distortion object.
Int_t nShots
Number of shots.
Bool_t pionAbsorption
True if the event is a pion absorption.
Class that stores isotopic abundances for a given element.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
void cascade()
The actual cascade loop.
Float_t impactParameter
Impact parameter [fm].
Solution solve(RootFunctor const *const f, const G4double x0)
Numerically solve a one-dimensional equation.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Float_t EBalance
Energy-conservation balance [MeV].
Short_t Sp
Strangeness number of the projectile nucleus.
Int_t nPionAbsorptions
Number of nucleon absorptions (no outcoming pions)
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
G4int makeProjectileRemnant()
Make a projectile pre-fragment out of geometrical spectators.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
Short_t Zt
Charge number of the target nucleus.
Class to adjust remnant recoil in the reaction CM system.
Bool_t forcedDeltasOutside
Event involved forced delta decays outside the nucleus.
void afterPropagationAction(IPropagationModel *pm, IAvatar *avatar)
double A(double temperature)
void deleteBlockers()
Delete blockers.
G4int getZ() const
Returns the charge number.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
Bool_t antikaonsInside
Event involved antikaons in the nucleus at the end of the cascade.
void clearIncoming()
Clear the incoming list.
void beforeRunAction(Config const *config)
G4int getA() const
Returns the baryon number.
LocalEnergyType getLocalEnergyPiType() const
Get the type of local energy for pi-N and decay avatars.
G4bool decayMe()
Force the phase-space decay of the Nucleus.
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
void emitInsideStrangeParticles()
Force emission of all strange particles inside the nucleus.
Float_t eventBias
Event bias.
G4bool isEventTransparent() const
Is the event transparent?
G4bool containsSigma()
Returns true if the nucleus contains any Sigma.
virtual G4double getStoppingTime()=0
void reset()
Reset the projectile remnant to the state at the beginning of the cascade.
IsotopeVector::iterator IsotopeIter
std::vector< Int_t > initialRandomSeeds
Initial seeds for the pseudo-random-number generator.
G4bool emitInsideKaon()
Force emission of all Kaon inside the nucleus.
Simple container for output of calculation-wide results.
G4bool isNaturalTarget() const
Natural targets.
Float_t pTransBalance
Transverse momentum-conservation balance [MeV/c].
Struct for conservation laws.
G4double getCutNN() const
static G4ThreadLocal Int_t eventNumber
Number of the event.
const EventInfo & processEvent()
Float_t Ep
Projectile kinetic energy given as input.
void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy)
Initialise the maximum interaction distance.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Bool_t emitKaon
Event involved forced Kaon emission.
IPropagationModel * propagationModel
void updateGlobalInfo()
Update global counters and other members of theGlobalInfo object.
void deleteIncoming()
Clear the incoming list and delete the particles.
G4bool decayOutgoingSigmaZero(G4double timeThreshold)
Force the decay of outgoing Neutral Sigma.
Bool_t forcedPionResonancesOutside
Event involved forced eta/omega decays outside the nucleus.
Simple class for computing intersections between a straight line and a sphere.
Int_t nForcedCompoundNucleus
Number of forced compound-nucleus events.
G4double interactionDistancePiN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
Short_t nCascadeParticles
Number of cascade particles.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
void initialize(Config const *const theConfig)
Initialize the Coulomb-distortion algorithm.
Int_t emitLambda
Number of forced Lambda emit out of the nucleus.
Float_t errorForcedCNCrossSection
Error on the calculated forced-compound-nucleus cross section.
G4bool containsDeltas()
Returns true if the nucleus contains any deltas.
Bool_t deltasInside
Event involved deltas in the nucleus at the end of the cascade.
Short_t Ap
Mass number of the projectile nucleus.
void postCascade()
Finalise the cascade and clean up.
Float_t nucleonAbsorptionCrossSection
Nucleon absorption cross section.
G4double interactionDistanceKbarN(const G4double projectileKineticEnergy)
Compute the "interaction distance".
void setEnergy(G4double energy)
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
std::string cascadeModel
Name of the cascade model.
Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r)
Compute the first intersection of a straight particle trajectory with a sphere.
std::string deexcitationModel
Name of the de-excitation model.
Bool_t forcedDeltasInside
Event involved forced delta decays inside the nucleus.
void clearCache()
Clear the INuclearPotential cache.
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
Class containing default actions to be performed at intermediate cascade steps.
virtual G4double getCurrentTime()=0
void distortOut(ParticleList const &pL, Nucleus const *const n)
Modify the momentum of an outgoing particle.
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
static G4ThreadLocal G4int nextBiasedCollisionID
Float_t stoppingTime
Cascade stopping time [fm/c].
void initialize(Config const *const theConfig)
Short_t Zp
Charge number of the projectile nucleus.
std::vector< Int_t > finalRandomSeeds
Final seeds for the pseudo-random-number generator.
G4bool decayOutgoingNeutralKaon()
Force the transformation of outgoing Neutral Kaon into propation eigenstate.
Int_t projectileType
Projectile particle type.
G4bool containsAntiKaon()
Returns true if the nucleus contains any anti Kaons.
Float_t energyViolationInteractionCrossSection
Cross section for attempted collisions/decays for which the energy-conservation algorithm failed to f...
Static class for selecting Coulomb distortion.
Short_t Ap
Projectile mass number given as input.
Short_t Zp
Projectile charge number given as input.
G4bool containsLambda()
Returns true if the nucleus contains any Lambda.
Bool_t lambdasInside
Event involved lambdas in the nucleus at the end of the cascade.
Bool_t absorbedStrangeParticle
Event involved forced strange absorbtion inside the nucleus.
Float_t pionAbsorptionCrossSection
Pion absorption cross section.
void initialize(Config const *const theConfig)
ParticleList const & getOutgoingParticles() const
G4double getBias() const
Get the bias.
void setA(const G4int A)
Set the mass number of the cluster.
ParticleList extractDynamicalSpectators()
Returns a list of dynamical spectators.
Abstract interface to the nuclear potential.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
void beforeCascadeAction(IPropagationModel *)
std::vector< Isotope > IsotopeVector
G4double maxImpactParameter
Short_t At
Mass number of the target nucleus.
void initialize(Config const *const theConfig=0)
Initialize the particle table.
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
Float_t geometricCrossSection
Geometric cross section.
G4double maxUniverseRadius
G4int getCascading() const
G4bool getTryCompoundNucleus()
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
G4bool continueCascade()
Stopping criterion for the cascade.
Int_t nNucleonAbsorptions
Number of nucleon absorptions (no outcoming particles)
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void deletePhaseSpaceGenerator()
static void setBias(const G4double b)
Set the bias.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
Bool_t clusterDecay
Event involved cluster decay.
G4bool initializeTarget(const G4int A, const G4int Z, const G4int S)
void afterAvatarAction(IAvatar *a, Nucleus *n, FinalState *fs)
Bool_t transparent
True if the event is transparent.
Float_t forcedCNCrossSection
Calculated forced-compound-nucleus cross section.