33 #define INCLXX_IN_GEANT4_MODE 1
88 :propagationModel(0), theA(208), theZ(82),
89 targetInitSuccess(false),
90 maxImpactParameter(0.),
91 maxUniverseRadius(0.),
92 maxInteractionDistance(0.),
93 fixedImpactParameter(0.),
96 forceTransparent(false),
100 #ifdef INCLXX_IN_GEANT4_MODE
102 #else // INCLXX_IN_GEANT4_MODE
105 #endif // INCLXX_IN_GEANT4_MODE
110 #ifdef INCLXX_IN_GEANT4_MODE
112 #else // INCLXX_IN_GEANT4_MODE
114 #endif // INCLXX_IN_GEANT4_MODE
166 cascadeAction->beforeRunAction(
theConfig);
174 #ifndef INCLXX_IN_GEANT4_MODE
191 #ifndef INCLXX_IN_GEANT4_MODE
192 G4INCL::NuclearMassTable::deleteTable();
199 #ifndef INCLXX_IN_GEANT4_MODE
200 G4INCL::Logger::deleteLoggerSlave();
211 if(A < 0 || A > 300 || Z < 1 || Z > 200) {
212 INCL_ERROR(
"Unsupported target: A = " << A <<
" Z = " << Z << std::endl
213 <<
"Target configuration rejected." << std::endl);
217 (projectileSpecies.
theZ==projectileSpecies.
theA || projectileSpecies.
theZ==0)) {
218 INCL_ERROR(
"Unsupported projectile: A = " << projectileSpecies.
theA <<
" Z = " << projectileSpecies.
theZ << std::endl
219 <<
"Projectile configuration rejected." << std::endl);
249 if(projectileSpecies.
theA > 0)
278 INCL_WARN(
"Target initialisation failed for A=" << targetA <<
", Z=" << targetZ << std::endl);
327 INCL_DEBUG(
"Selected impact parameter: " << impactParameter << std::endl);
333 if(effectiveImpactParameter < 0.) {
359 if(avatar == 0)
break;
392 INCL_DEBUG(
"Trying compound nucleus" << std::endl);
396 #ifndef INCLXX_IN_GEANT4_MODE
406 if(projectileRemnant) {
434 INCL_DEBUG(
"Cascade resulted in complete fusion, using realistic fusion kinematics" << std::endl);
440 INCL_WARN(
"Complete-fusion kinematics yields negative excitation energy, returning a transparent!" << std::endl);
457 #ifndef INCLXX_IN_GEANT4_MODE
459 globalConservationChecks(
false);
471 #ifndef INCLXX_IN_GEANT4_MODE
473 globalConservationChecks(
true);
512 std::vector<Particle *> shuffledComponents(initialProjectileComponents.begin(), initialProjectileComponents.end());
517 G4bool atLeastOneNucleonEntering =
false;
518 for(std::vector<Particle*>::const_iterator p=shuffledComponents.begin(), e=shuffledComponents.end(); p!=e; ++p) {
522 (*p)->getPropagationVelocity(),
524 if(!intersectionInteractionDistance.
exists)
528 atLeastOneNucleonEntering =
true;
541 theCNZ += (*p)->getZ();
551 if(!success || !atLeastOneNucleonEntering) {
552 INCL_DEBUG(
"No nucleon entering in forced CN, forcing a transparent" << std::endl);
562 theCNEnergy -= theProjectileRemnant->
getEnergy();
563 theCNMomentum -= theProjectileRemnant->
getMomentum();
571 for(
ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
572 theCNSpin -= (*i)->getAngularMomentum();
577 const G4double theCNInvariantMassSquared = theCNEnergy*theCNEnergy-theCNMomentum.
mag2();
578 if(theCNInvariantMassSquared<0.) {
583 const G4double theCNExcitationEnergy = std::sqrt(theCNInvariantMassSquared) - theCNMass;
584 if(theCNExcitationEnergy<0.) {
586 INCL_DEBUG(
"CN excitation energy is negative, forcing a transparent" << std::endl
587 <<
" theCNA = " << theCNA << std::endl
588 <<
" theCNZ = " << theCNZ << std::endl
589 <<
" theCNEnergy = " << theCNEnergy << std::endl
590 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" << std::endl
591 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy << std::endl
592 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" << std::endl
598 INCL_DEBUG(
"CN excitation energy is positive, forcing a CN" << std::endl
599 <<
" theCNA = " << theCNA << std::endl
600 <<
" theCNZ = " << theCNZ << std::endl
601 <<
" theCNEnergy = " << theCNEnergy << std::endl
602 <<
" theCNMomentum = (" << theCNMomentum.
getX() <<
", "<< theCNMomentum.
getY() <<
", " << theCNMomentum.
getZ() <<
")" << std::endl
603 <<
" theCNExcitationEnergy = " << theCNExcitationEnergy << std::endl
604 <<
" theCNSpin = (" << theCNSpin.
getX() <<
", "<< theCNSpin.
getY() <<
", " << theCNSpin.
getZ() <<
")" << std::endl
631 theRecoilFunctor(theSolution.
x);
633 INCL_WARN(
"Couldn't accommodate remnant recoil while satisfying energy conservation, root-finding algorithm failed." << std::endl);
638 #ifndef INCLXX_IN_GEANT4_MODE
639 void INCL::globalConservationChecks(
G4bool afterRecoil) {
645 if(theBalance.
Z != 0) {
646 INCL_ERROR(
"Violation of charge conservation! ZBalance = " << theBalance.
Z << std::endl);
648 if(theBalance.
A != 0) {
649 INCL_ERROR(
"Violation of baryon-number conservation! ABalance = " << theBalance.
A << std::endl);
651 G4double EThreshold, pLongThreshold, pTransThreshold;
656 pTransThreshold = 1.;
660 pLongThreshold = 0.1;
661 pTransThreshold = 0.1;
663 if(std::abs(theBalance.
energy)>EThreshold) {
664 INCL_WARN(
"Violation of energy conservation > " << EThreshold <<
" MeV. EBalance = " << theBalance.
energy <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber << std::endl);
666 if(std::abs(pLongBalance)>pLongThreshold) {
667 INCL_WARN(
"Violation of longitudinal momentum conservation > " << pLongThreshold <<
" MeV/c. pLongBalance = " << pLongBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber << std::endl);
669 if(std::abs(pTransBalance)>pTransThreshold) {
670 INCL_WARN(
"Violation of transverse momentum conservation > " << pTransThreshold <<
" MeV/c. pTransBalance = " << pTransBalance <<
" afterRecoil = " << afterRecoil <<
" eventNumber=" <<
theEventInfo.
eventNumber << std::endl);
685 <<
"), stopping cascade" << std::endl);
691 INCL_DEBUG(
"No participants in the nucleus and no incoming particles left, stopping cascade" << std::endl);
698 <<
"), stopping cascade" << std::endl);
704 INCL_DEBUG(
"Trying to make a compound nucleus, stopping cascade" << std::endl);
743 G4int nUnmergedSpectators = 0;
746 if(dynSpectators.empty() && geomSpectators.empty()) {
748 }
else if(dynSpectators.size()==1 && geomSpectators.empty()) {
759 nUnmergedSpectators = rejected.size();
767 return nUnmergedSpectators;
781 INCL_DEBUG(
"Initialised interaction distance: r0 = " << r0 << std::endl
782 <<
" theNNDistance = " << theNNDistance << std::endl
792 for(
IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
796 rMax =
std::max(maximumRadius, rMax);
801 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.
G4int shuffleComponentsHelper(G4int range)
Helper function for ProjectileRemnant::shuffleStoredComponents.
Float_t effectiveImpactParameter
Effective (Coulomb-distorted) impact parameter [fm].
G4double maxUniverseRadius
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()
static void distortOut(ParticleList const &pL, Nucleus const *const n)
Modify the momentum of an outgoing particle.
FinalStateValidity getValidity() const
void deleteGenerator()
Delete the generator.
Class for non-relativistic Coulomb distortion.
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.
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.
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
virtual G4INCL::IAvatar * propagate()=0
Propagate the particles and get the next avatar.
ParticleList addAllDynamicalSpectators(ParticleList pL)
Add back all dynamical spectators to the projectile remnant.
G4bool isNucleusNucleusCollision() const
Is it a nucleus-nucleus collision?
Abstract interface to the nuclear potential.
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.
const ThreeVector & getIncomingAngularMomentum() const
Get the incoming angular-momentum vector.
static void setCoulomb(ICoulomb *const coulomb)
Set the Coulomb-distortion algorithm.
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.
Bool_t forcedDeltasOutside
Event involved forced delta decays outside the nucleus.
G4double getEnergy() const
Get the energy of the particle in MeV.
Struct for conservation laws.
Bool_t pionAbsorption
True if the event is a pion absorption.
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.
Abstract interface for the cross-section classes.
Short_t At
Target mass number given as input.
static void deleteBlockers()
Delete blockers.
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.
void finalizeGlobalInfo()
CascadeAction * cascadeAction
Float_t stoppingTime
Cascade stopping time [fm/c].
G4double mag2() const
Get the square of the length.
static void setBlocker(IPauli *const)
Set the Pauli blocker algorithm.
Float_t Ep
Projectile kinetic energy given as input.
void setCrossSections(ICrossSections *c)
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.
static G4ThreadLocal Int_t eventNumber
Number of the event.
void clearCache()
Clear the INuclearPotential cache.
G4int getVerbosity() const
Get the verbosity.
void beforeAvatarAction(IAvatar *a, Nucleus *n)
G4double getImpactParameter() const
G4bool preCascade(ParticleSpecies const projectileSpecies, const G4double kineticEnergy)
Initialise the cascade.
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.
Cross sections used in INCL4.6.
G4bool getCDPP() const
Do we want CDPP?
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.
static G4double maxImpactParameter(ParticleSpecies const &p, const G4double kinE, Nucleus const *const n)
Return the maximum impact parameter for Coulomb-distorted trajectories.
static void deleteCoulomb()
Delete the Coulomb-distortion object.
Simple container for output of calculation-wide results.
Float_t pTransBalance
Transverse momentum-conservation balance [MeV/c].
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
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.
PauliType getPauliType() const
Get the Pauli-blocking algorithm.
void beforePropagationAction(IPropagationModel *pm)
void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z)
Initialize the universe radius.
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.
Placeholder class for no Coulomb distortion.
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 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.
void setGenerator(G4INCL::IRandomGenerator *aGenerator)
Set the random number generator implementation to be used globally by INCL.
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.
void initVerbosityLevelFromEnvvar()
void addParticleEntryAvatar(IAvatar *a)
Add one ParticleEntry avatar.
Float_t errorCompleteFusionCrossSection
Error on the calculated complete-fusion cross section (nParticles==0)
Random::SeedVector getRandomSeeds() const
Get the seeds for the random-number generator.
void setA(const G4int A)
Set the mass number of the cluster.
ClusterAlgorithmType getClusterAlgorithm() const
Get the clustering algorithm.
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)
std::string const & getLogFileName() const
Get the log file name.
Static class for cluster formation.
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.
virtual G4double shoot(ParticleSpecies const projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi)=0
LocalEnergyType getLocalEnergyBBType() const
Get the type of local energy for N-N avatars.
void reset()
Reset the EventInfo members.
CoulombType getCoulombType() const
Get the Coulomb-distortion algorithm.
static void setClusteringModel(IClusteringModel *const model)
Set the clustering model.
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?
G4INCL::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.
Abstract interface for Coulomb distortion.
static void deleteClusteringModel()
Delete clustering model.
Short_t Zt
Target charge number given as input.
static void setCDPP(IPauli *const)
Set the CDPP blocker algorithm.
Intersection-point structure.
Bool_t transparent
True if the event is transparent.
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
Cluster coalescence algorithm used in the IAEA intercomparison.
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.
Cross sections used in INCL4.6.
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
ParticleList::const_iterator ParticleIter
Simple class for computing intersections between a straight line and a sphere.
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.