34 #define INCLXX_IN_GEANT4_MODE 1
49 #include "G4String.hh"
62 thePreCompoundModel(aPreCompound),
65 complainedAboutBackupModel(false),
66 complainedAboutPreCompound(false),
68 theINCLXXLevelDensity(NULL),
69 theINCLXXFissionProbability(NULL)
79 if(getenv(
"G4INCLXX_NO_DE_EXCITATION")) {
80 G4String message =
"de-excitation is completely disabled!";
93 if(theFissionChannelCast) {
107 if(getenv(
"G4INCLXX_DUMP_REMNANT"))
135 std::stringstream ss;
136 ss <<
"the model does not know how to handle a collision between a "
158 if(pA > theMaxProjMassINCL)
160 else if(tA > theMaxProjMassINCL)
177 if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
186 if(trackA<=1 && nucleusA<=1) {
193 if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
196 std::stringstream ss;
197 ss <<
"INCL++ refuses to handle reactions between nuclei with A>"
198 << theMaxProjMassINCL
199 <<
". A backup model ("
201 <<
") will be used instead.";
211 && trackKinE < cascadeMinEnergyPerNucleon) {
214 std::stringstream ss;
215 ss <<
"INCL++ refuses to handle nucleon-induced reactions below "
216 << cascadeMinEnergyPerNucleon /
MeV
217 <<
" MeV. A PreCoumpound model ("
219 <<
") will be used instead.";
228 const G4double theTrackEnergy = trackKinE + theTrackMass;
229 const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
230 const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
235 fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
236 fourMomentumIn.setVect(theTrackMomentum);
245 G4Nucleus *theTargetNucleus = &theNucleus;
246 if(inverseKinematics) {
250 if(oldProjectileDef != 0 && oldTargetDef != 0) {
254 if(newTargetA > 0 && newTargetZ > 0) {
256 theTargetNucleus =
new G4Nucleus(newTargetA, newTargetZ);
259 G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
262 G4String message =
"badly defined target after swapping. Falling back to normal (non-swapped) mode.";
267 G4String message =
"oldProjectileDef or oldTargetDef was null";
288 toZ.rotateZ(-projectileMomentum.phi());
289 toZ.rotateY(-projectileMomentum.theta());
303 std::list<G4Fragment> remnants;
305 const G4int maxTries = 200;
325 if(inverseKinematics) {
344 momentum *= toLabFrame;
346 if(inverseKinematics) {
347 momentum *= *toDirectKinematics;
348 momentum.setVect(-momentum.vect());
353 fourMomentumOut += momentum;
357 G4String message =
"the model produced a particle that couldn't be converted to Geant4 particle.";
371 eventInfo.
jxRem[i]*hbar_Planck,
372 eventInfo.
jyRem[i]*hbar_Planck,
373 eventInfo.
jzRem[i]*hbar_Planck
382 if(std::abs(scaling - 1.0) > 0.01) {
383 std::stringstream ss;
384 ss <<
"momentum scaling = " << scaling
385 <<
"\n Lorentz vector = " << fourMomentum
386 <<
")\n A = " << A <<
", Z = " << Z
387 <<
"\n E* = " << excitationE <<
", nuclearMass = " << nuclearMass
388 <<
"\n remnant i=" << i <<
", nRemnants=" << eventInfo.
nRemnants
392 <<
", in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
397 fourMomentum *= toLabFrame;
400 if(inverseKinematics) {
401 fourMomentum *= *toDirectKinematics;
402 fourMomentum.setVect(-fourMomentum.vect());
405 fourMomentumOut += fourMomentum;
409 G4cerr <<
"G4INCLXX_DUMP_REMNANT: " << remnant <<
" spin: " << spin <<
G4endl;
411 remnants.push_back(remnant);
415 const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
416 const G4double energyViolation = std::abs(violation4Momentum.e());
417 const G4double momentumViolation = violation4Momentum.rho();
419 std::stringstream ss;
420 ss <<
"energy conservation violated by " << energyViolation/
MeV <<
" MeV in "
423 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
427 for(
G4int j=0; j<nSecondaries; ++j)
433 std::stringstream ss;
434 ss <<
"momentum conservation violated by " << momentumViolation/
MeV <<
" MeV in "
437 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics. Will resample.";
441 for(
G4int j=0; j<nSecondaries; ++j)
449 }
while(!eventIsOK && nTries < maxTries);
452 if(inverseKinematics) {
453 delete aProjectileTrack;
454 delete theTargetNucleus;
455 delete toInverseKinematics;
456 delete toDirectKinematics;
460 std::stringstream ss;
461 ss <<
"maximum number of tries exceeded for the proposed "
464 <<
" inelastic reaction, in " << (inverseKinematics ?
"inverse" :
"direct") <<
" kinematics.";
475 for(std::list<G4Fragment>::iterator i = remnants.begin();
476 i != remnants.end(); ++i) {
479 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
480 fragment != deExcitationResult->end(); ++fragment) {
488 for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
489 fragment != deExcitationResult->end(); ++fragment) {
492 deExcitationResult->clear();
493 delete deExcitationResult;
549 else if(A == 3 && Z == 2)
return G4He3::He3();
551 else if(A > 0 && Z > 0 && A > Z) {
577 const G4double p2 = px*px + py*py + pz*pz;
579 const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
580 return std::sqrt(pnew2)/std::sqrt(p2);
588 <<
"The Liège Intranuclear Cascade (INCL++) is a model for reactions induced\n"
589 <<
"by nucleons, pions and light ion on any nucleus. The reaction is\n"
590 <<
"described as an avalanche of binary nucleon-nucleon collisions, which can\n"
591 <<
"lead to the emission of energetic particles and to the formation of an\n"
592 <<
"excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
593 <<
"outside the scope of INCL++ and is typically described by another model.\n\n"
594 <<
"INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
595 <<
"pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
596 <<
"Most tests involved target nuclei close to the stability valley, with\n"
597 <<
"numbers between 4 and 250.\n\n"
598 <<
"Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
Short_t nRemnants
Number of remnants.
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetMaxProjMassINCL() const
Getter for theMaxProjMassINCL.
G4double GetCascadeMinEnergyPerNucleon() const
Getter for cascadeMinEnergyPerNucleon.
G4HadronicInteraction * theBackupModel
G4HadSecondary * GetSecondary(size_t i)
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepRotation G4RotationMatrix
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
G4VEvaporationChannel * GetFissionChannel()
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Proton * ProtonDefinition()
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
G4DynamicParticle * toG4Particle(G4int A, G4int Z, G4double kinE, G4double px, G4double py, G4double pz) const
Convert an INCL particle to a G4DynamicParticle.
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
G4INCLXXInterface(G4VPreCompoundModel *const aPreCompound=0)
const G4String & GetModelName() const
G4VEvaporation * GetEvaporation()
void EmitBigWarning(const G4String &message) const
Emit a BIG warning to G4cout.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
G4INCL::ParticleType toINCLParticleType(G4ParticleDefinition const *const) const
Convert G4ParticleDefinition to corresponding INCL particle type.
G4bool complainedAboutPreCompound
void SetEmissionStrategy(G4VEmissionProbability *aFissionProb)
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
Header file for the G4INCLXXInterfaceStore class.
const G4String & GetParticleName() const
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
const G4String & GetIonName(G4int Z, G4int A, G4int lvl=0) const
G4HadronicInteraction * theBackupModelNucleon
void SetFissionLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)
G4int GetAtomicNumber() const
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
G4VPreCompoundModel * thePreCompoundModel
G4bool complainedAboutBackupModel
void SetStatusChange(G4HadFinalStateStatus aS)
Short_t ARem[maxSizeRemnants]
Remnant mass number.
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4double remnant4MomentumScaling(G4double mass, G4double kineticE, G4double px, G4double py, G4double pz) const
Rescale remnant momentum if necessary.
static G4INCLXXInterfaceStore * GetInstance()
Get the singleton instance.
G4ExcitationHandler * GetExcitationHandler() const
Singleton class for configuring the INCL++ Geant4 interface.
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
Main method to apply the INCL physics model.
double A(double temperature)
Short_t Z[maxSizeParticles]
Particle charge number.
const G4ParticleDefinition * GetDefinition() const
Short_t nParticles
Number of particles in the final state.
G4HadFinalState theResult
G4INCL::INCL * GetINCLModel()
Get the cached INCL model engine.
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
const EventInfo & processEvent()
G4double GetKineticEnergy() const
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
const G4String & GetParticleType() const
G4bool AccurateProjectile(const G4HadProjectile &aTrack, const G4Nucleus &theTargetNucleus) const
static G4Neutron * Neutron()
Revised level-density parameter for fission after INCL++.
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
static G4Deuteron * Deuteron()
G4int GetAtomicMass() const
G4ParticleDefinition * toG4ParticleDefinition(G4int A, G4int Z) const
Convert A and Z to a G4ParticleDefinition.
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
G4LorentzVector Get4Momentum() const
G4String const & GetDeExcitationModelName() const
G4INCL::ParticleSpecies toINCLParticleSpecies(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle species.
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
G4HadronicInteraction * FindModel(const G4String &name)
Short_t A[maxSizeParticles]
Particle mass number.
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)=0
static G4GenericIon * GenericIon()
G4double GetPDGMass() const
virtual void ModelDescription(std::ostream &outFile) const
G4double toINCLKineticEnergy(G4HadProjectile const &) const
Convert G4HadProjectile to corresponding INCL particle kinetic energy.
G4double energy(const ThreeVector &p, const G4double m)
G4INCL::INCL * theINCLModel
static G4HadronicInteractionRegistry * Instance()
G4DynamicParticle * GetParticle()
static G4PionMinus * PionMinus()
G4bool GetAccurateProjectile() const
Getter for accurateProjectile.
G4INCLXXVInterfaceTally * theTally
G4VPreCompoundModel * theDeExcitation
void SetAngularMomentum(G4ThreeVector &)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
std::map< G4String, G4AttDef > * GetInstance(const G4String &storeKey, G4bool &isNew)
G4IonTable *const theIonTable
Bool_t transparent
True if the event is transparent.
void EmitWarning(const G4String &message)
Emit a warning to G4cout.
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
G4double GetPDGCharge() const
void SetLevelDensityParameter(G4VLevelDensityParameter *aLevelDensity)
virtual void Tally(G4HadProjectile const &aTrack, G4Nucleus const &theNucleus, G4HadFinalState const &result)=0
G4INCLXXVInterfaceTally * GetTally() const
Getter for the interface tally.
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
static G4Neutron * NeutronDefinition()
void SetMomentumChange(const G4ThreeVector &aV)
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4GLOB_DLL std::ostream G4cerr
G4VLevelDensityParameter * theINCLXXLevelDensity
CLHEP::HepLorentzVector G4LorentzVector
G4FissionProbability * theINCLXXFissionProbability
G4INCLXXInterfaceStore *const theInterfaceStore