34#define INCLXX_IN_GEANT4_MODE 1
71 thePreCompoundModel(aPreCompound),
74 complainedAboutBackupModel(false),
75 complainedAboutPreCompound(false),
77 theINCLXXLevelDensity(NULL),
78 theINCLXXFissionProbability(NULL),
89 if(std::getenv(
"G4INCLXX_NO_DE_EXCITATION")) {
90 G4String message =
"de-excitation is completely disabled!";
103 if(theFissionChannelCast) {
117 if(std::getenv(
"G4INCLXX_DUMP_REMNANT"))
142 std::stringstream ss;
143 ss <<
"the model does not know how to handle a collision between a "
165 if(pA > theMaxProjMassINCL)
167 else if(tA > theMaxProjMassINCL)
185 if((isIonTrack && ((trackZ<=0 && trackL==0) || trackA<=trackZ)) ||
186 (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
195 if(trackA<=1 && nucleusA<=1) {
202 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
205 std::stringstream ss;
206 ss <<
"INCL++ refuses to handle reactions between nuclei with A>"
207 << theMaxProjMassINCL
208 <<
". A backup model ("
210 <<
") will be used instead.";
220 && trackKinE < cascadeMinEnergyPerNucleon) {
223 std::stringstream ss;
224 ss <<
"INCL++ refuses to handle nucleon-induced reactions below "
225 << cascadeMinEnergyPerNucleon /
MeV
226 <<
" MeV. A PreCoumpound model ("
228 <<
") will be used instead.";
237 const G4double theTrackEnergy = trackKinE + theTrackMass;
238 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
239 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
244 fourMomentumIn.
setE(theTrackEnergy + theNucleusMass);
245 fourMomentumIn.
setVect(theTrackMomentum);
254 G4Nucleus *theTargetNucleus = &theNucleus;
255 if(inverseKinematics) {
259 if(oldProjectileDef != 0 && oldTargetDef != 0) {
263 if(newTargetA > 0 && newTargetZ > 0) {
265 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ, newTargetL);
268 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
271 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
276 G4String message =
"oldProjectileDef or oldTargetDef was null";
312 std::list<G4Fragment> remnants;
314 const G4int maxTries = 200;
330 -theTargetNucleus->
GetL());
337 if(inverseKinematics) {
358 momentum *= toLabFrame;
360 if(inverseKinematics) {
361 momentum *= *toDirectKinematics;
367 fourMomentumOut += momentum;
371 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
401 if(std::abs(scaling - 1.0) > 0.01) {
402 std::stringstream ss;
403 ss <<
"momentum scaling = " << scaling
404 <<
"\n Lorentz vector = " << fourMomentum
405 <<
")\n A = " <<
A <<
", Z = " <<
Z <<
", S = " <<
S
406 <<
"\n E* = " << excitationE <<
", nuclearMass = " << nuclearMass
407 <<
"\n remnant i=" << i <<
", nRemnants=" << eventInfo.
nRemnants
411 <<
", in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
416 fourMomentum *= toLabFrame;
419 if(inverseKinematics) {
420 fourMomentum *= *toDirectKinematics;
424 fourMomentumOut += fourMomentum;
429 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
431 remnants.push_back(remnant);
435 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
436 const G4double energyViolation = std::abs(violation4Momentum.
e());
437 const G4double momentumViolation = violation4Momentum.
rho();
439 std::stringstream ss;
440 ss <<
"energy conservation violated by " << energyViolation/
MeV <<
" MeV in "
443 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
447 for(
G4int j=0; j<nSecondaries; ++j)
453 std::stringstream ss;
454 ss <<
"momentum conservation violated by " << momentumViolation/
MeV <<
" MeV in "
457 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
461 for(
G4int j=0; j<nSecondaries; ++j)
469 }
while(!eventIsOK && nTries < maxTries);
472 if(inverseKinematics) {
473 delete aProjectileTrack;
474 delete theTargetNucleus;
475 delete toInverseKinematics;
476 delete toDirectKinematics;
480 std::stringstream ss;
481 ss <<
"maximum number of tries exceeded for the proposed "
484 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
495 for(std::list<G4Fragment>::iterator i = remnants.begin();
496 i != remnants.end(); ++i) {
499 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
500 fragment != deExcitationResult->end(); ++fragment) {
508 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
509 fragment != deExcitationResult->end(); ++fragment) {
512 deExcitationResult->clear();
513 delete deExcitationResult;
575 }
else if(PDGCode == 221) {
return G4Eta::Eta();
589 }
else if(PDGCode == 2003) {
return G4He3::He3();
598 }
else if(
A > 0 &&
Z > 0 &&
A >
Z) {
622 const G4double p2 = px*px + py*py + pz*pz;
624 const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
625 return std::sqrt(pnew2)/std::sqrt(p2);
633 <<
"The Li�ge Intranuclear Cascade (INCL++) is a model for reactions induced\n"
634 <<
"by nucleons, pions and light ion on any nucleus. The reaction is\n"
635 <<
"described as an avalanche of binary nucleon-nucleon collisions, which can\n"
636 <<
"lead to the emission of energetic particles and to the formation of an\n"
637 <<
"excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
638 <<
"outside the scope of INCL++ and is typically described by another model.\n\n"
639 <<
"INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
640 <<
"pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
641 <<
"Most tests involved target nuclei close to the stability valley, with\n"
642 <<
"numbers between 4 and 250.\n\n"
643 <<
"Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
G4double S(G4double temp)
Header file for the G4INCLXXInterfaceStore class.
CLHEP::HepLorentzRotation G4LorentzRotation
std::vector< G4ReactionProduct * > G4ReactionProductVector
static constexpr double MeV
G4GLOB_DLL std::ostream G4cerr
HepLorentzRotation inverse() const
Hep3Vector boostVector() const
void setVect(const Hep3Vector &)
HepRotation inverse() const
HepRotation & rotateZ(double delta)
HepRotation & rotateY(double delta)
static G4AntiKaonZero * AntiKaonZero()
void SetEmissionStrategy(G4VEmissionProbability *aFissionProb)
void SetLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)
static G4Deuteron * Deuteron()
static G4DoubleHyperDoubleNeutron * Definition()
static G4DoubleHyperH4 * Definition()
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4VEvaporation * GetEvaporation()
Revised level-density parameter for fission after INCL++.
void SetFissionLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)
void SetCreatorModelID(G4int value)
void SetAngularMomentum(const G4ThreeVector &)
static G4GenericIon * GenericIon()
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
std::size_t GetNumberOfSecondaries() const
void SetEnergyChange(G4double anEnergy)
G4HadSecondary * GetSecondary(size_t i)
void SetMomentumChange(const G4ThreeVector &aV)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4DynamicParticle * GetParticle()
G4HadronicInteraction * FindModel(const G4String &name)
static G4HadronicInteractionRegistry * Instance()
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
const G4String & GetModelName() const
static G4HyperAlpha * Definition()
static G4HyperH4 * Definition()
static G4HyperHe5 * Definition()
static G4double GetNuclearMass(G4int A, G4int Z, G4int L)
static G4HyperTriton * Definition()
Singleton class for configuring the INCL++ Geant4 interface.
void EmitWarning(const G4String &message)
Emit a warning to G4cout.
static G4INCLXXInterfaceStore * GetInstance()
Get the singleton instance.
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4INCL::INCL * GetINCLModel()
Get the cached INCL model engine.
G4double GetConservationTolerance() const
Getter for conservationTolerance.
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
G4INCLXXVInterfaceTally * GetTally() const
Getter for the interface tally.
G4bool GetAccurateProjectile() const
Getter for accurateProjectile.
G4bool complainedAboutPreCompound
virtual void ModelDescription(std::ostream &outFile) const
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
G4INCLXXInterfaceStore *const theInterfaceStore
G4HadFinalState theResult
G4INCL::ParticleType toINCLParticleType(G4ParticleDefinition const *const) const
Convert G4ParticleDefinition to corresponding INCL particle type.
G4ParticleDefinition * toG4ParticleDefinition(G4int A, G4int Z, G4int S, G4int PDGCode) const
Convert A, Z and S to a G4ParticleDefinition.
G4double remnant4MomentumScaling(G4double mass, G4double kineticE, G4double px, G4double py, G4double pz) const
Rescale remnant momentum if necessary.
G4INCL::ParticleSpecies toINCLParticleSpecies(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle species.
G4double toINCLKineticEnergy(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle kinetic energy.
G4INCLXXInterface(G4VPreCompoundModel *const aPreCompound=0)
G4bool complainedAboutBackupModel
G4FissionProbability * theINCLXXFissionProbability
G4DynamicParticle * toG4Particle(G4int A, G4int Z, G4int S, G4int PDGCode, G4double kinE, G4double px, G4double py, G4double pz) const
Convert an INCL particle to a G4DynamicParticle.
G4HadronicInteraction * theBackupModelNucleon
G4VLevelDensityParameter * theINCLXXLevelDensity
G4String const & GetDeExcitationModelName() const
G4HadronicInteraction * theBackupModel
G4bool AccurateProjectile(const G4HadProjectile &aTrack, const G4Nucleus &theTargetNucleus) const
G4INCLXXVInterfaceTally * theTally
G4IonTable *const theIonTable
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
G4VPreCompoundModel * thePreCompoundModel
G4INCL::INCL * theINCLModel
virtual void Tally(G4HadProjectile const &aTrack, G4Nucleus const &theNucleus, G4HadFinalState const &result)=0
const EventInfo & processEvent()
const G4String & GetIonName(G4int Z, G4int A, G4int lvl=0) const
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4double GetIonMass(G4int Z, G4int A, G4int nL=0, G4int lvl=0) const
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4KaonZero * KaonZero()
static G4Lambda * Lambda()
static G4Neutron * NeutronDefinition()
static G4Neutron * Neutron()
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetAtomicNumber() const
const G4String & GetParticleType() const
G4double GetPDGMass() const
G4int GetAtomicMass() const
G4double GetPDGCharge() const
G4int GetNumberOfLambdasInHypernucleus() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
static G4int GetModelID(const G4int modelIndex)
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4PionZero * PionZero()
static G4Proton * ProtonDefinition()
static G4Proton * Proton()
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4SigmaZero * SigmaZero()
static G4Triton * Triton()
G4VEvaporationChannel * GetFissionChannel()
G4VPreCompoundModel * theDeExcitation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
G4ExcitationHandler * GetExcitationHandler() const
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
G4double energy(const ThreeVector &p, const G4double m)
Short_t S[maxSizeParticles]
Particle strangeness number.
Short_t nRemnants
Number of remnants.
Bool_t transparent
True if the event is transparent.
Short_t A[maxSizeParticles]
Particle mass number.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
Short_t Z[maxSizeParticles]
Particle charge number.
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
Short_t nParticles
Number of particles in the final state.
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
Short_t SRem[maxSizeRemnants]
Remnant strangeness number.
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
Int_t PDGCode[maxSizeParticles]
PDG numbering of the particles.
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
Short_t ARem[maxSizeRemnants]
Remnant mass number.
Short_t ZRem[maxSizeRemnants]
Remnant charge number.