INCL++ intra-nuclear cascade. More...

#include <G4INCLXXInterface.hh>

Inheritance diagram for G4INCLXXInterface:

Collaboration diagram for G4INCLXXInterface:

Public Member Functions
	G4INCLXXInterface (G4VPreCompoundModel *const aPreCompound=0)

	~G4INCLXXInterface ()

G4int	operator== (G4INCLXXInterface &right)

G4int	operator!= (G4INCLXXInterface &right)

G4ReactionProductVector *	Propagate (G4KineticTrackVector theSecondaries, G4V3DNucleus theNucleus)

G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)

void	DeleteModel ()

virtual void	ModelDescription (std::ostream &outFile) const

G4String const &	GetDeExcitationModelName () const

Public Member Functions inherited from G4VIntraNuclearTransportModel
	G4VIntraNuclearTransportModel (const G4String &modelName="CascadeModel", G4VPreCompoundModel *ptr=0)

virtual	~G4VIntraNuclearTransportModel ()

virtual G4ReactionProductVector *	PropagateNuclNucl (G4KineticTrackVector theSecondaries, G4V3DNucleus theNucleus, G4V3DNucleus *theProjectileNucleus)

void	SetDeExcitation (G4VPreCompoundModel *ptr)

void	Set3DNucleus (G4V3DNucleus *const value)

void	SetPrimaryProjectile (const G4HadProjectile &aPrimary)

const G4String &	GetModelName () const

virtual void	PropagateModelDescription (std::ostream &outFile) const

Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")

virtual	~G4HadronicInteraction ()

virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)

virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)

G4double	GetMinEnergy () const

G4double	GetMinEnergy (const G4Material aMaterial, const G4Element anElement) const

void	SetMinEnergy (G4double anEnergy)

void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)

void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)

G4double	GetMaxEnergy () const

G4double	GetMaxEnergy (const G4Material aMaterial, const G4Element anElement) const

void	SetMaxEnergy (const G4double anEnergy)

void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)

void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)

const G4HadronicInteraction *	GetMyPointer () const

virtual G4int	GetVerboseLevel () const

virtual void	SetVerboseLevel (G4int value)

const G4String &	GetModelName () const

void	DeActivateFor (const G4Material *aMaterial)

void	ActivateFor (const G4Material *aMaterial)

void	DeActivateFor (const G4Element *anElement)

void	ActivateFor (const G4Element *anElement)

G4bool	IsBlocked (const G4Material *aMaterial) const

G4bool	IsBlocked (const G4Element *anElement) const

void	SetRecoilEnergyThreshold (G4double val)

G4double	GetRecoilEnergyThreshold () const

G4bool	operator== (const G4HadronicInteraction &right) const

G4bool	operator!= (const G4HadronicInteraction &right) const

virtual const std::pair< G4double, G4double >	GetFatalEnergyCheckLevels () const

virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const

void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)

virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

Private Member Functions
G4bool	AccurateProjectile (const G4HadProjectile &aTrack, const G4Nucleus &theTargetNucleus) const

	G4INCLXXInterface (const G4INCLXXInterface &rhs)
	Dummy copy constructor to shut up Coverity warnings. More...

G4INCLXXInterface &	operator= (G4INCLXXInterface const &rhs)
	Dummy assignment operator to shut up Coverity warnings. More...

G4INCL::ParticleType	toINCLParticleType (G4ParticleDefinition const *const) const
	Convert G4ParticleDefinition to corresponding INCL particle type. More...

G4INCL::ParticleSpecies	toINCLParticleSpecies (G4HadProjectile const &) const
	Convert G4HadProjectile to corresponding INCL particle species. More...

G4double	toINCLKineticEnergy (G4HadProjectile const &) const
	Convert G4HadProjectile to corresponding INCL particle kinetic energy. More...

G4DynamicParticle *	toG4Particle (G4int A, G4int Z, G4double kinE, G4double px, G4double py, G4double pz) const
	Convert an INCL particle to a G4DynamicParticle. More...

G4ParticleDefinition *	toG4ParticleDefinition (G4int A, G4int Z) const
	Convert A and Z to a G4ParticleDefinition. More...

G4double	remnant4MomentumScaling (G4double mass, G4double kineticE, G4double px, G4double py, G4double pz) const
	Rescale remnant momentum if necessary. More...

Private Attributes
G4INCL::INCL *	theINCLModel

G4VPreCompoundModel *	thePreCompoundModel

G4HadFinalState	theResult

G4HadronicInteraction *	theBackupModel

G4HadronicInteraction *	theBackupModelNucleon

G4INCLXXInterfaceStore *const	theInterfaceStore

G4INCLXXVInterfaceTally *	theTally

G4bool	complainedAboutBackupModel

G4bool	complainedAboutPreCompound

G4IonTable *const	theIonTable

G4bool	dumpRemnantInfo

G4VLevelDensityParameter *	theINCLXXLevelDensity

G4FissionProbability *	theINCLXXFissionProbability

Additional Inherited Members
Protected Member Functions inherited from G4VIntraNuclearTransportModel
G4V3DNucleus *	Get3DNucleus () const

G4VPreCompoundModel *	GetDeExcitation () const

const G4HadProjectile *	GetPrimaryProjectile () const

Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)

G4bool	IsBlocked () const

void	Block ()

Protected Attributes inherited from G4VIntraNuclearTransportModel
G4String	theTransportModelName

G4V3DNucleus *	the3DNucleus

G4VPreCompoundModel *	theDeExcitation

const G4HadProjectile *	thePrimaryProjectile

Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange

G4int	verboseLevel

G4double	theMinEnergy

G4double	theMaxEnergy

G4bool	isBlocked

Detailed Description

INCL++ intra-nuclear cascade.

Interface for INCL++. This interface handles basic hadron bullet particles (protons, neutrons, pions), as well as light ions.

Example usage in case of protons:

G4INCLXXInterface* inclModel = new G4INCLXXInterface;
inclModel -> SetMinEnergy(0.0 * MeV); // Set the energy limits
inclModel -> SetMaxEnergy(3.0 * GeV);
G4ProtonInelasticProcess* protonInelasticProcess = new G4ProtonInelasticProcess();
G4ProtonInelasticCrossSection* protonInelasticCrossSection =  new G4ProtonInelasticCrossSection();
protonInelasticProcess -> RegisterMe(inclModel);
protonInelasticProcess -> AddDataSet(protonInelasticCrossSection);
particle = G4Proton::Proton();
processManager = particle -> GetProcessManager();
processManager -> AddDiscreteProcess(protonInelasticProcess);

The same setup procedure is needed for neutron, pion and generic-ion inelastic processes as well.

Definition at line 103 of file G4INCLXXInterface.hh.

Constructor & Destructor Documentation

◆ G4INCLXXInterface() [1/2]

G4INCLXXInterface::G4INCLXXInterface ( G4VPreCompoundModel *const aPreCompound = 0 )

Definition at line 59 of file G4INCLXXInterface.cc.

                                                                              :
   G4VIntraNuclearTransportModel(G4INCLXXInterfaceStore::GetInstance()->getINCLXXVersionName()),
   theINCLModel(NULL),
   thePreCompoundModel(aPreCompound),
   theInterfaceStore(G4INCLXXInterfaceStore::GetInstance()),
   theTally(NULL),
   complainedAboutBackupModel(false),
   complainedAboutPreCompound(false),
   theIonTable(G4IonTable::GetIonTable()),
   theINCLXXLevelDensity(NULL),
   theINCLXXFissionProbability(NULL)
 {
   if(!thePreCompoundModel) {
     G4HadronicInteraction* p =
       G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
     thePreCompoundModel = static_cast<G4VPreCompoundModel*>(p);
     if(!thePreCompoundModel) { thePreCompoundModel = new G4PreCompoundModel; }
   }
 
   // Use the environment variable G4INCLXX_NO_DE_EXCITATION to disable de-excitation
   if(getenv("G4INCLXX_NO_DE_EXCITATION")) {
     G4String message = "de-excitation is completely disabled!";
     theInterfaceStore->EmitWarning(message);
     theDeExcitation = 0;
   } else {
     G4HadronicInteraction* p =
       G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
     theDeExcitation = static_cast<G4VPreCompoundModel*>(p);
     if(!theDeExcitation) { theDeExcitation = new G4PreCompoundModel; }
 
     // set the fission parameters for G4ExcitationHandler
     G4VEvaporationChannel * const theFissionChannel =
       theDeExcitation->GetExcitationHandler()->GetEvaporation()->GetFissionChannel();
     G4CompetitiveFission * const theFissionChannelCast = dynamic_cast<G4CompetitiveFission *>(theFissionChannel);
     if(theFissionChannelCast) {
       theINCLXXLevelDensity = new G4FissionLevelDensityParameterINCLXX;
       theFissionChannelCast->SetLevelDensityParameter(theINCLXXLevelDensity);
       theINCLXXFissionProbability = new G4FissionProbability;
       theINCLXXFissionProbability->SetFissionLevelDensityParameter(theINCLXXLevelDensity);
       theFissionChannelCast->SetEmissionStrategy(theINCLXXFissionProbability);
       theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler uses its own level-density parameter for fission");
     } else {
       theInterfaceStore->EmitBigWarning("INCL++/G4ExcitationHandler could not use its own level-density parameter for fission");
     }
   }
 
   // use the envvar G4INCLXX_DUMP_REMNANT to dump information about the
   // remnants on stdout
   if(getenv("G4INCLXX_DUMP_REMNANT"))
     dumpRemnantInfo = true;
   else
     dumpRemnantInfo = false;
 
   theBackupModel = new G4BinaryLightIonReaction;
   theBackupModelNucleon = new G4BinaryCascade;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ~G4INCLXXInterface()

G4INCLXXInterface::~G4INCLXXInterface ( )

Definition at line 116 of file G4INCLXXInterface.cc.

 {
   delete theINCLXXLevelDensity;
   delete theINCLXXFissionProbability;
 }

◆ G4INCLXXInterface() [2/2]

G4INCLXXInterface::G4INCLXXInterface ( const G4INCLXXInterface & rhs )

private

Dummy copy constructor to shut up Coverity warnings.

Member Function Documentation

◆ AccurateProjectile()

G4bool G4INCLXXInterface::AccurateProjectile	(	const G4HadProjectile &	aTrack,
		const G4Nucleus &	theTargetNucleus
	)		const

private

Definition at line 122 of file G4INCLXXInterface.cc.

                                                                                                              {
   // Use direct kinematics if the projectile is a nucleon or a pion
   const G4ParticleDefinition *projectileDef = aTrack.GetDefinition();
   if(projectileDef == G4Proton::Proton()
      || projectileDef == G4Neutron::Neutron()
      || projectileDef == G4PionPlus::PionPlus()
      || projectileDef == G4PionZero::PionZero()
      || projectileDef == G4PionMinus::PionMinus())
     return false;
 
   // Here all projectiles should be nuclei
   const G4int pA = projectileDef->GetAtomicMass();
   if(pA<=0) {
     std::stringstream ss;
     ss << "the model does not know how to handle a collision between a "
       << projectileDef->GetParticleName() << " projectile and a Z="
       << theNucleus.GetZ_asInt() << ", A=" << theNucleus.GetA_asInt();
     theInterfaceStore->EmitBigWarning(ss.str());
     return true;
   }
 
   // If either nucleus is a LCP (A<=4), run the collision as light on heavy
   const G4int tA = theNucleus.GetA_asInt();
   if(tA<=4 || pA<=4) {
     if(pA<tA)
       return false;
     else
       return true;
   }
 
   // If one of the nuclei is heavier than theMaxProjMassINCL, run the collision
   // as light on heavy.
   // Note that here we are sure that either the projectile or the target is
   // smaller than theMaxProjMass; otherwise theBackupModel would have been
   // called.
   const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
   if(pA > theMaxProjMassINCL)
     return true;
   else if(tA > theMaxProjMassINCL)
     return false;
   else
     // In all other cases, use the global setting
     return theInterfaceStore->GetAccurateProjectile();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ApplyYourself()

G4HadFinalState * G4INCLXXInterface::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	theNucleus
	)

virtual

Main method to apply the INCL physics model.

Parameters

aTrack	the projectile particle
theNucleus	target nucleus

Returns: the output of the INCL physics model

Implements G4HadronicInteraction.

Definition at line 167 of file G4INCLXXInterface.cc.

 {
   G4ParticleDefinition const * const trackDefinition = aTrack.GetDefinition();
   const G4bool isIonTrack = trackDefinition->GetParticleType()==G4GenericIon::GenericIon()->GetParticleType();
   const G4int trackA = trackDefinition->GetAtomicMass();
   const G4int trackZ = (G4int) trackDefinition->GetPDGCharge();
   const G4int nucleusA = theNucleus.GetA_asInt();
   const G4int nucleusZ = theNucleus.GetZ_asInt();
 
   // For reactions induced by weird projectiles (e.g. He2), bail out
   if((isIonTrack && (trackZ<=0 || trackA<=trackZ)) || (nucleusA>1 && (nucleusZ<=0 || nucleusA<=nucleusZ))) {
     theResult.Clear();
     theResult.SetStatusChange(isAlive);
     theResult.SetEnergyChange(aTrack.GetKineticEnergy());
     theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
     return &theResult;
   }
 
   // For reactions on nucleons, use the backup model (without complaining)
   if(trackA<=1 && nucleusA<=1) {
     return theBackupModelNucleon->ApplyYourself(aTrack, theNucleus);
   }
 
   // For systems heavier than theMaxProjMassINCL, use another model (typically
   // BIC)
   const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
   if(trackA > theMaxProjMassINCL && nucleusA > theMaxProjMassINCL) {
     if(!complainedAboutBackupModel) {
       complainedAboutBackupModel = true;
       std::stringstream ss;
       ss << "INCL++ refuses to handle reactions between nuclei with A>"
         << theMaxProjMassINCL
         << ". A backup model ("
         << theBackupModel->GetModelName()
         << ") will be used instead.";
       theInterfaceStore->EmitBigWarning(ss.str());
     }
     return theBackupModel->ApplyYourself(aTrack, theNucleus);
   }
 
   // For energies lower than cascadeMinEnergyPerNucleon, use PreCompound
   const G4double cascadeMinEnergyPerNucleon = theInterfaceStore->GetCascadeMinEnergyPerNucleon();
   const G4double trackKinE = aTrack.GetKineticEnergy();
   if((trackDefinition==G4Neutron::NeutronDefinition() || trackDefinition==G4Proton::ProtonDefinition())
       && trackKinE < cascadeMinEnergyPerNucleon) {
     if(!complainedAboutPreCompound) {
       complainedAboutPreCompound = true;
       std::stringstream ss;
       ss << "INCL++ refuses to handle nucleon-induced reactions below "
         << cascadeMinEnergyPerNucleon / MeV
         << " MeV. A PreCoumpound model ("
         << thePreCompoundModel->GetModelName()
         << ") will be used instead.";
       theInterfaceStore->EmitBigWarning(ss.str());
     }
     return thePreCompoundModel->ApplyYourself(aTrack, theNucleus);
   }
 
   // Calculate the total four-momentum in the entrance channel
   const G4double theNucleusMass = theIonTable->GetIonMass(nucleusZ, nucleusA);
   const G4double theTrackMass = trackDefinition->GetPDGMass();
   const G4double theTrackEnergy = trackKinE + theTrackMass;
   const G4double theTrackMomentumAbs2 = theTrackEnergy*theTrackEnergy - theTrackMass*theTrackMass;
   const G4double theTrackMomentumAbs = ((theTrackMomentumAbs2>0.0) ? std::sqrt(theTrackMomentumAbs2) : 0.0);
   const G4ThreeVector theTrackMomentum = aTrack.Get4Momentum().getV().unit() * theTrackMomentumAbs;
 
   G4LorentzVector goodTrack4Momentum(theTrackMomentum, theTrackEnergy);
   G4LorentzVector fourMomentumIn;
   fourMomentumIn.setE(theTrackEnergy + theNucleusMass);
   fourMomentumIn.setVect(theTrackMomentum);
 
   // Check if inverse kinematics should be used
   const G4bool inverseKinematics = AccurateProjectile(aTrack, theNucleus);
 
   // If we are running in inverse kinematics, redefine aTrack and theNucleus
   G4LorentzRotation *toInverseKinematics = NULL;
   G4LorentzRotation *toDirectKinematics = NULL;
   G4HadProjectile const *aProjectileTrack = &aTrack;
   G4Nucleus *theTargetNucleus = &theNucleus;
   if(inverseKinematics) {
     G4ParticleDefinition *oldTargetDef = theIonTable->GetIon(nucleusZ, nucleusA, 0);
     const G4ParticleDefinition *oldProjectileDef = trackDefinition;
 
     if(oldProjectileDef != 0 && oldTargetDef != 0) {
       const G4int newTargetA = oldProjectileDef->GetAtomicMass();
       const G4int newTargetZ = oldProjectileDef->GetAtomicNumber();
 
       if(newTargetA > 0 && newTargetZ > 0) {
         // This should give us the same energy per nucleon
         theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ);
         toInverseKinematics = new G4LorentzRotation(goodTrack4Momentum.boostVector());
         G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theNucleusMass);
         G4DynamicParticle swappedProjectileParticle(oldTargetDef, (*toInverseKinematics) * theProjectile4Momentum);
         aProjectileTrack = new G4HadProjectile(swappedProjectileParticle);
       } else {
         G4String message = "badly defined target after swapping. Falling back to normal (non-swapped) mode.";
         theInterfaceStore->EmitWarning(message);
         toInverseKinematics = new G4LorentzRotation;
       }
     } else {
       G4String message = "oldProjectileDef or oldTargetDef was null";
       theInterfaceStore->EmitWarning(message);
       toInverseKinematics = new G4LorentzRotation;
     }
   }
 
   // INCL assumes the projectile particle is going in the direction of
   // the Z-axis. Here we construct proper rotation to convert the
   // momentum vectors of the outcoming particles to the original
   // coordinate system.
   G4LorentzVector projectileMomentum = aProjectileTrack->Get4Momentum();
 
   // INCL++ assumes that the projectile is going in the direction of
   // the z-axis. In principle, if the coordinate system used by G4
   // hadronic framework is defined differently we need a rotation to
   // transform the INCL++ reaction products to the appropriate
   // frame. Please note that it isn't necessary to apply this
   // transform to the projectile because when creating the INCL++
   // projectile particle; \see{toINCLKineticEnergy} needs to use only the
   // projectile energy (direction is simply assumed to be along z-axis).
   G4RotationMatrix toZ;
   toZ.rotateZ(-projectileMomentum.phi());
   toZ.rotateY(-projectileMomentum.theta());
   G4RotationMatrix toLabFrame3 = toZ.inverse();
   G4LorentzRotation toLabFrame(toLabFrame3);
   // However, it turns out that the projectile given to us by G4
   // hadronic framework is already going in the direction of the
   // z-axis so this rotation is actually unnecessary. Both toZ and
   // toLabFrame turn out to be unit matrices as can be seen by
   // uncommenting the folowing two lines:
   // G4cout <<"toZ = " << toZ << G4endl;
   // G4cout <<"toLabFrame = " << toLabFrame << G4endl;
 
   theResult.Clear(); // Make sure the output data structure is clean.
   theResult.SetStatusChange(stopAndKill);
 
   std::list<G4Fragment> remnants;
 
   const G4int maxTries = 200;
   G4int nTries = 0;
   // INCL can generate transparent events. However, this is meaningful
   // only in the standalone code. In Geant4 we must "force" INCL to
   // produce a valid cascade.
   G4bool eventIsOK = false;
   do {
     const G4INCL::ParticleSpecies theSpecies = toINCLParticleSpecies(*aProjectileTrack);
     const G4double kineticEnergy = toINCLKineticEnergy(*aProjectileTrack);
 
     // The INCL model will be created at the first use
     theINCLModel = G4INCLXXInterfaceStore::GetInstance()->GetINCLModel();
 
     const G4INCL::EventInfo eventInfo = theINCLModel->processEvent(theSpecies, kineticEnergy, theTargetNucleus->GetA_asInt(), theTargetNucleus->GetZ_asInt());
     //    eventIsOK = !eventInfo.transparent && nTries < maxTries;
     eventIsOK = !eventInfo.transparent;
     if(eventIsOK) {
 
       // If the collision was run in inverse kinematics, prepare to boost back
       // all the reaction products
       if(inverseKinematics) {
         toDirectKinematics = new G4LorentzRotation(toInverseKinematics->inverse());
       }
 
       G4LorentzVector fourMomentumOut;
 
       for(G4int i = 0; i < eventInfo.nParticles; ++i) {
     G4int A = eventInfo.A[i];
     G4int Z = eventInfo.Z[i];
     //  G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
     G4double kinE = eventInfo.EKin[i];
     G4double px = eventInfo.px[i];
     G4double py = eventInfo.py[i];
     G4double pz = eventInfo.pz[i];
     G4DynamicParticle *p = toG4Particle(A, Z , kinE, px, py, pz);
     if(p != 0) {
       G4LorentzVector momentum = p->Get4Momentum();
 
           // Apply the toLabFrame rotation
           momentum *= toLabFrame;
           // Apply the inverse-kinematics boost
           if(inverseKinematics) {
             momentum *= *toDirectKinematics;
             momentum.setVect(-momentum.vect());
           }
 
       // Set the four-momentum of the reaction products
       p->Set4Momentum(momentum);
           fourMomentumOut += momentum;
       theResult.AddSecondary(p);
 
     } else {
       G4String message = "the model produced a particle that couldn't be converted to Geant4 particle.";
           theInterfaceStore->EmitWarning(message);
     }
       }
 
       for(G4int i = 0; i < eventInfo.nRemnants; ++i) {
     const G4int A = eventInfo.ARem[i];
     const G4int Z = eventInfo.ZRem[i];
     //  G4cout <<"INCL particle A = " << A << " Z = " << Z << G4endl;
     const G4double kinE = eventInfo.EKinRem[i];
     const G4double px = eventInfo.pxRem[i];
     const G4double py = eventInfo.pyRem[i];
     const G4double pz = eventInfo.pzRem[i];
         G4ThreeVector spin(
             eventInfo.jxRem[i]*hbar_Planck,
             eventInfo.jyRem[i]*hbar_Planck,
             eventInfo.jzRem[i]*hbar_Planck
             );
     const G4double excitationE = eventInfo.EStarRem[i];
     const G4double nuclearMass = G4NucleiProperties::GetNuclearMass(A, Z) + excitationE;
         const G4double scaling = remnant4MomentumScaling(nuclearMass,
             kinE,
             px, py, pz);
     G4LorentzVector fourMomentum(scaling * px, scaling * py, scaling * pz,
                      nuclearMass + kinE);
     if(std::abs(scaling - 1.0) > 0.01) {
           std::stringstream ss;
           ss << "momentum scaling = " << scaling
             << "\n                Lorentz vector = " << fourMomentum
             << ")\n                A = " << A << ", Z = " << Z
             << "\n                E* = " << excitationE << ", nuclearMass = " << nuclearMass
             << "\n                remnant i=" << i << ", nRemnants=" << eventInfo.nRemnants
             << "\n                Reaction was: " << aTrack.GetKineticEnergy()/MeV
             << "-MeV " << trackDefinition->GetParticleName() << " + "
             << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
             << ", in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
           theInterfaceStore->EmitWarning(ss.str());
     }
 
         // Apply the toLabFrame rotation
         fourMomentum *= toLabFrame;
         spin *= toLabFrame3;
         // Apply the inverse-kinematics boost
         if(inverseKinematics) {
           fourMomentum *= *toDirectKinematics;
           fourMomentum.setVect(-fourMomentum.vect());
         }
 
         fourMomentumOut += fourMomentum;
     G4Fragment remnant(A, Z, fourMomentum);
         remnant.SetAngularMomentum(spin);
         if(dumpRemnantInfo) {
           G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << "  spin: " << spin << G4endl;
         }
     remnants.push_back(remnant);
       }
 
       // Check four-momentum conservation
       const G4LorentzVector violation4Momentum = fourMomentumOut - fourMomentumIn;
       const G4double energyViolation = std::abs(violation4Momentum.e());
       const G4double momentumViolation = violation4Momentum.rho();
       if(energyViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
         std::stringstream ss;
         ss << "energy conservation violated by " << energyViolation/MeV << " MeV in "
           << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
           << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
           << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
         theInterfaceStore->EmitWarning(ss.str());
         eventIsOK = false;
         const G4int nSecondaries = theResult.GetNumberOfSecondaries();
         for(G4int j=0; j<nSecondaries; ++j)
           delete theResult.GetSecondary(j)->GetParticle();
         theResult.Clear();
         theResult.SetStatusChange(stopAndKill);
         remnants.clear();
       } else if(momentumViolation > G4INCLXXInterfaceStore::GetInstance()->GetConservationTolerance()) {
         std::stringstream ss;
         ss << "momentum conservation violated by " << momentumViolation/MeV << " MeV in "
           << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
           << " + " << theIonTable->GetIonName(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0)
           << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics. Will resample.";
         theInterfaceStore->EmitWarning(ss.str());
         eventIsOK = false;
         const G4int nSecondaries = theResult.GetNumberOfSecondaries();
         for(G4int j=0; j<nSecondaries; ++j)
           delete theResult.GetSecondary(j)->GetParticle();
         theResult.Clear();
         theResult.SetStatusChange(stopAndKill);
         remnants.clear();
       }
     }
     nTries++;
   } while(!eventIsOK && nTries < maxTries); /* Loop checking, 10.07.2015, D.Mancusi */
 
   // Clean up the objects that we created for the inverse kinematics
   if(inverseKinematics) {
     delete aProjectileTrack;
     delete theTargetNucleus;
     delete toInverseKinematics;
     delete toDirectKinematics;
   }
 
   if(!eventIsOK) {
     std::stringstream ss;
     ss << "maximum number of tries exceeded for the proposed "
     << aTrack.GetKineticEnergy()/MeV << "-MeV " << trackDefinition->GetParticleName()
     << " + " << theIonTable->GetIonName(nucleusZ, nucleusA, 0)
     << " inelastic reaction, in " << (inverseKinematics ? "inverse" : "direct") << " kinematics.";
     theInterfaceStore->EmitWarning(ss.str());
     theResult.SetStatusChange(isAlive);
     theResult.SetEnergyChange(aTrack.GetKineticEnergy());
     theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
     return &theResult;
   }
 
   // De-excitation:
 
   if(theDeExcitation != 0) {
     for(std::list<G4Fragment>::iterator i = remnants.begin();
     i != remnants.end(); ++i) {
       G4ReactionProductVector *deExcitationResult = theDeExcitation->DeExcite((*i));
 
       for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
       fragment != deExcitationResult->end(); ++fragment) {
     const G4ParticleDefinition *def = (*fragment)->GetDefinition();
     if(def != 0) {
       G4DynamicParticle *theFragment = new G4DynamicParticle(def, (*fragment)->GetMomentum());
       theResult.AddSecondary(theFragment);
     }
       }
 
       for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
       fragment != deExcitationResult->end(); ++fragment) {
     delete (*fragment);
       }
       deExcitationResult->clear();
       delete deExcitationResult;
     }
   }
 
   remnants.clear();
 
   if((theTally = theInterfaceStore->GetTally()))
     theTally->Tally(aTrack, theNucleus, theResult);
 
   return &theResult;
 }

Here is the caller graph for this function:

◆ DeleteModel()

void G4INCLXXInterface::DeleteModel ( )

inline

Definition at line 128 of file G4INCLXXInterface.hh.

                      {
     delete theINCLModel;
     theINCLModel = NULL;
   }

Here is the call graph for this function:

◆ GetDeExcitationModelName()

G4String const & G4INCLXXInterface::GetDeExcitationModelName ( ) const

Definition at line 601 of file G4INCLXXInterface.cc.

                                                                   {
   return theDeExcitation->GetModelName();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ModelDescription()

void G4INCLXXInterface::ModelDescription ( std::ostream & outFile ) const

virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 586 of file G4INCLXXInterface.cc.

                                                                   {
    outFile
      << "The Liège Intranuclear Cascade (INCL++) is a model for reactions induced\n"
      << "by nucleons, pions and light ion on any nucleus. The reaction is\n"
      << "described as an avalanche of binary nucleon-nucleon collisions, which can\n"
      << "lead to the emission of energetic particles and to the formation of an\n"
      << "excited thermalised nucleus (remnant). The de-excitation of the remnant is\n"
      << "outside the scope of INCL++ and is typically described by another model.\n\n"
      << "INCL++ has been reasonably well tested for nucleon (~50 MeV to ~15 GeV),\n"
      << "pion (idem) and light-ion projectiles (up to A=18, ~10A MeV to 1A GeV).\n"
      << "Most tests involved target nuclei close to the stability valley, with\n"
      << "numbers between 4 and 250.\n\n"
      << "Reference: D. Mancusi et al., Phys. Rev. C90 (2014) 054602\n\n";
 }

Here is the caller graph for this function:

◆ operator!=()

G4int G4INCLXXInterface::operator!= ( G4INCLXXInterface & right )

inline

Definition at line 112 of file G4INCLXXInterface.hh.

                                              {
     return (this != &right);
   }

Here is the call graph for this function:

◆ operator=()

G4INCLXXInterface& G4INCLXXInterface::operator= ( G4INCLXXInterface const & rhs )

private

Dummy assignment operator to shut up Coverity warnings.

Here is the caller graph for this function:

◆ operator==()

G4int G4INCLXXInterface::operator== ( G4INCLXXInterface & right )

inline

Definition at line 108 of file G4INCLXXInterface.hh.

                                              {
     return (this == &right);
   }

◆ Propagate()

G4ReactionProductVector * G4INCLXXInterface::Propagate	(	G4KineticTrackVector *	theSecondaries,
		G4V3DNucleus *	theNucleus
	)

virtual

Implements G4VIntraNuclearTransportModel.

Definition at line 505 of file G4INCLXXInterface.cc.

                                                                                             {
   return 0;
 }

Here is the caller graph for this function:

◆ remnant4MomentumScaling()

G4double G4INCLXXInterface::remnant4MomentumScaling	(	G4double	mass,
		G4double	kineticE,
		G4double	px,
		G4double	py,
		G4double	pz
	)		const

private

Rescale remnant momentum if necessary.

Definition at line 573 of file G4INCLXXInterface.cc.

                                              {
   const G4double p2 = px*px + py*py + pz*pz;
   if(p2 > 0.0) {
     const G4double pnew2 = kineticE*kineticE + 2.0*kineticE*mass;
     return std::sqrt(pnew2)/std::sqrt(p2);
   } else {
     return 1.0;
   }
 }

Here is the caller graph for this function:

◆ toG4Particle()

G4DynamicParticle * G4INCLXXInterface::toG4Particle	(	G4int	A,
		G4int	Z,
		G4double	kinE,
		G4double	px,
		G4double	py,
		G4double	pz
	)		const

private

Convert an INCL particle to a G4DynamicParticle.

Definition at line 558 of file G4INCLXXInterface.cc.

                                                                                  {
   const G4ParticleDefinition *def = toG4ParticleDefinition(A, Z);
   if(def == 0) { // Check if we have a valid particle definition
     return 0;
   }
   const G4double energy = kinE * MeV;
   const G4ThreeVector momentum(px, py, pz);
   const G4ThreeVector momentumDirection = momentum.unit();
   G4DynamicParticle *p = new G4DynamicParticle(def, momentumDirection, energy);
   return p;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ toG4ParticleDefinition()

G4ParticleDefinition * G4INCLXXInterface::toG4ParticleDefinition	(	G4int	A,
		G4int	Z
	)		const

private

Convert A and Z to a G4ParticleDefinition.

Definition at line 541 of file G4INCLXXInterface.cc.

                                                                                       {
   if     (A == 1 && Z == 1)  return G4Proton::Proton();
   else if(A == 1 && Z == 0)  return G4Neutron::Neutron();
   else if(A == 0 && Z == 1)  return G4PionPlus::PionPlus();
   else if(A == 0 && Z == -1) return G4PionMinus::PionMinus();
   else if(A == 0 && Z == 0)  return G4PionZero::PionZero();
   else if(A == 2 && Z == 1)  return G4Deuteron::Deuteron();
   else if(A == 3 && Z == 1)  return G4Triton::Triton();
   else if(A == 3 && Z == 2)  return G4He3::He3();
   else if(A == 4 && Z == 2)  return G4Alpha::Alpha();
   else if(A > 0 && Z > 0 && A > Z) { // Returns ground state ion definition
     return theIonTable->GetIon(Z, A, 0);
   } else { // Error, unrecognized particle
     return 0;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ toINCLKineticEnergy()

G4double G4INCLXXInterface::toINCLKineticEnergy ( G4HadProjectile const & aTrack ) const

private

Convert G4HadProjectile to corresponding INCL particle kinetic energy.

Definition at line 537 of file G4INCLXXInterface.cc.

                                                                                    {
   return aTrack.GetKineticEnergy();
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ toINCLParticleSpecies()

G4INCL::ParticleSpecies G4INCLXXInterface::toINCLParticleSpecies ( G4HadProjectile const & aTrack ) const

private

Convert G4HadProjectile to corresponding INCL particle species.

Definition at line 523 of file G4INCLXXInterface.cc.

                                                                                                   {
   const G4ParticleDefinition *pdef = aTrack.GetDefinition();
   const G4INCL::ParticleType theType = toINCLParticleType(pdef);
   if(theType!=G4INCL::Composite)
     return G4INCL::ParticleSpecies(theType);
   else {
     G4INCL::ParticleSpecies theSpecies;
     theSpecies.theType=theType;
     theSpecies.theA=pdef->GetAtomicMass();
     theSpecies.theZ=pdef->GetAtomicNumber();
     return theSpecies;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ toINCLParticleType()

G4INCL::ParticleType G4INCLXXInterface::toINCLParticleType ( G4ParticleDefinition const * const pdef ) const

private

Convert G4ParticleDefinition to corresponding INCL particle type.

Definition at line 509 of file G4INCLXXInterface.cc.

                                                                                                       {
   if(     pdef == G4Proton::Proton())           return G4INCL::Proton;
   else if(pdef == G4Neutron::Neutron())         return G4INCL::Neutron;
   else if(pdef == G4PionPlus::PionPlus())       return G4INCL::PiPlus;
   else if(pdef == G4PionMinus::PionMinus())     return G4INCL::PiMinus;
   else if(pdef == G4PionZero::PionZero())       return G4INCL::PiZero;
   else if(pdef == G4Deuteron::Deuteron())       return G4INCL::Composite;
   else if(pdef == G4Triton::Triton())           return G4INCL::Composite;
   else if(pdef == G4He3::He3())                 return G4INCL::Composite;
   else if(pdef == G4Alpha::Alpha())             return G4INCL::Composite;
   else if(pdef->GetParticleType() == G4GenericIon::GenericIon()->GetParticleType()) return G4INCL::Composite;
   else                                            return G4INCL::UnknownParticle;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

◆ complainedAboutBackupModel

G4bool G4INCLXXInterface::complainedAboutBackupModel

private

Definition at line 178 of file G4INCLXXInterface.hh.

◆ complainedAboutPreCompound

G4bool G4INCLXXInterface::complainedAboutPreCompound

private

Definition at line 179 of file G4INCLXXInterface.hh.

◆ dumpRemnantInfo

G4bool G4INCLXXInterface::dumpRemnantInfo

private

Definition at line 183 of file G4INCLXXInterface.hh.

◆ theBackupModel

G4HadronicInteraction* G4INCLXXInterface::theBackupModel

private

Definition at line 172 of file G4INCLXXInterface.hh.

◆ theBackupModelNucleon

G4HadronicInteraction* G4INCLXXInterface::theBackupModelNucleon

private

Definition at line 173 of file G4INCLXXInterface.hh.

◆ theINCLModel

G4INCL::INCL* G4INCLXXInterface::theINCLModel

private

Definition at line 166 of file G4INCLXXInterface.hh.

◆ theINCLXXFissionProbability

G4FissionProbability* G4INCLXXInterface::theINCLXXFissionProbability

private

Definition at line 186 of file G4INCLXXInterface.hh.

◆ theINCLXXLevelDensity

G4VLevelDensityParameter* G4INCLXXInterface::theINCLXXLevelDensity

private

Definition at line 185 of file G4INCLXXInterface.hh.

◆ theInterfaceStore

G4INCLXXInterfaceStore* const G4INCLXXInterface::theInterfaceStore

private

Definition at line 175 of file G4INCLXXInterface.hh.

◆ theIonTable

G4IonTable* const G4INCLXXInterface::theIonTable

private

Definition at line 181 of file G4INCLXXInterface.hh.

◆ thePreCompoundModel

G4VPreCompoundModel* G4INCLXXInterface::thePreCompoundModel

private

Definition at line 168 of file G4INCLXXInterface.hh.

◆ theResult

G4HadFinalState G4INCLXXInterface::theResult

private

Definition at line 170 of file G4INCLXXInterface.hh.

◆ theTally

G4INCLXXVInterfaceTally* G4INCLXXInterface::theTally

private

Definition at line 176 of file G4INCLXXInterface.hh.

The documentation for this class was generated from the following files:

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4INCLXXInterface() [1/2]

◆ ~G4INCLXXInterface()

◆ G4INCLXXInterface() [2/2]

Member Function Documentation

◆ AccurateProjectile()

◆ ApplyYourself()

◆ DeleteModel()

◆ GetDeExcitationModelName()

◆ ModelDescription()

◆ operator!=()

◆ operator=()

◆ operator==()

◆ Propagate()

◆ remnant4MomentumScaling()

◆ toG4Particle()

◆ toG4ParticleDefinition()

◆ toINCLKineticEnergy()

◆ toINCLParticleSpecies()

◆ toINCLParticleType()

Member Data Documentation

◆ complainedAboutBackupModel

◆ complainedAboutPreCompound

◆ dumpRemnantInfo

◆ theBackupModel

◆ theBackupModelNucleon

◆ theINCLModel

◆ theINCLXXFissionProbability

◆ theINCLXXLevelDensity

◆ theInterfaceStore

◆ theIonTable

◆ thePreCompoundModel

◆ theResult

◆ theTally