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

#include <G4INCLXXInterface.hh>

Inheritance diagram for G4INCLXXInterface:

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

	~G4INCLXXInterface ()

G4bool	operator== (G4INCLXXInterface &right)

G4bool	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

void	SetDeExcitation (G4VPreCompoundModel *ptr)

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

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

	G4VIntraNuclearTransportModel (const G4VIntraNuclearTransportModel &right)=delete

const G4VIntraNuclearTransportModel &	operator= (const G4VIntraNuclearTransportModel &right)=delete

G4bool	operator== (const G4VIntraNuclearTransportModel &right) const =delete

G4bool	operator!= (const G4VIntraNuclearTransportModel &right) const =delete

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 &aTrack, G4Nucleus &targetNucleus)

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)

G4int	GetVerboseLevel () const

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

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 &)

virtual void	InitialiseModel ()

	G4HadronicInteraction (const G4HadronicInteraction &right)=delete

const G4HadronicInteraction &	operator= (const G4HadronicInteraction &right)=delete

G4bool	operator== (const G4HadronicInteraction &right) const =delete

G4bool	operator!= (const G4HadronicInteraction &right) const =delete

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);
 
G4HadronInelasticProcess* protonInelasticProcess = new G4HadronInelasticProcess( "protonInelastic", G4Proton::Definition() );
G4VCrossSectionDataSet* protonInelasticCrossSection = new G4BGGNucleonInelasticXS( G4Proton::Proton() );
 
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()

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

Definition at line 70 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),
  secID(-1)
{
  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(std::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(std::getenv("G4INCLXX_DUMP_REMNANT"))
    dumpRemnantInfo = true;
  else
    dumpRemnantInfo = false;
 
  theBackupModel = new G4BinaryLightIonReaction;
  theBackupModelNucleon = new G4BinaryCascade;
  secID = G4PhysicsModelCatalog::GetModelID( "model_INCLXXCascade" );
}

◆ ~G4INCLXXInterface()

G4INCLXXInterface::~G4INCLXXInterface ( )

Definition at line 129 of file G4INCLXXInterface.cc.

{
  delete theINCLXXLevelDensity;
  delete theINCLXXFissionProbability;
}

Member Function Documentation

◆ 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

Reimplemented from G4HadronicInteraction.

Definition at line 176 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 trackL = trackDefinition->GetNumberOfLambdasInHypernucleus();
  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 && trackL==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),
  // except for anti_proton projectile (in this case, INCLXX is used).
  if(trackA<=1 && nucleusA<=1 && (trackZ>=0 || trackA==0)) {
    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();
      const G4int newTargetL = oldProjectileDef->GetNumberOfLambdasInHypernucleus();
      if(newTargetA > 0 && newTargetZ > 0) {
        // This should give us the same energy per nucleon
        theTargetNucleus = new G4Nucleus(newTargetA, newTargetZ, newTargetL);
        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(),
                   -theTargetNucleus->GetL());  // Strangeness has opposite sign
    //    eventIsOK = !eventInfo.transparent && nTries < maxTries;                              // of the number of Lambdas 
    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];
        G4int S = eventInfo.S[i];  // Strangeness
        G4int PDGCode = eventInfo.PDGCode[i];
        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, S, PDGCode, 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;
 
    // Propagate the particle's parent resonance information
          G4HadSecondary secondary(p, 1.0, secID);
          G4ParticleDefinition* parentResonanceDef = nullptr;
          if ( eventInfo.parentResonancePDGCode[i] != 0 ) {
            parentResonanceDef = G4ParticleTable::GetParticleTable()->FindParticle(eventInfo.parentResonancePDGCode[i]);
          }
          secondary.SetParentResonanceDef(parentResonanceDef);
          secondary.SetParentResonanceID(eventInfo.parentResonanceID[i]);
    
          theResult.AddSecondary(secondary);
 
        } 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];
        const G4int S = eventInfo.SRem[i];
        // Check that the remnant is a physical bound state: if not, resample the collision.
        if(( Z == 0  &&  S == 0  &&  A > 1 ) ||                 // No bound states for nn, nnn, nnnn, ...
           ( Z == 0  &&  S != 0  &&  A < 4 ) ||                 // No bound states for nl, ll, nnl, nll, lll
           ( Z != 0  &&  S != 0  &&  A == Z + std::abs(S) )) {  // No bound states for pl, ppl, pll, ...
          std::stringstream ss;
          ss << "unphysical residual fragment : Z=" << Z << "  S=" << S << "  A=" << A 
             << "  skipping it and resampling the collision";
          theInterfaceStore->EmitWarning(ss.str());
          eventIsOK = false;
          continue;
        }
        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];
        G4double nuclearMass = excitationE;
        if ( S == 0 ) {
          nuclearMass += G4NucleiProperties::GetNuclearMass(A, Z);
        } else {
    // Assumed that the opposite of the strangeness of the remnant gives the number of Lambdas inside it
          nuclearMass += G4HyperNucleiProperties::GetNuclearMass(A, Z, std::abs(S));
        }
        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 << ", S = " << S
             << "\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, std::abs(S), fourMomentum);  // Assumed that -strangeness gives the number of Lambdas
        remnant.SetAngularMomentum(spin);
        remnant.SetCreatorModelID(secID);
        if(dumpRemnantInfo) {
          G4cerr << "G4INCLXX_DUMP_REMNANT: " << remnant << "  spin: " << spin << G4endl;
        }
        remnants.push_back(remnant);
      }
 
      // Give up is the event is not ok (e.g. unphysical residual)
      if(!eventIsOK) {
        const G4int nSecondaries = (G4int)theResult.GetNumberOfSecondaries();
        for(G4int j=0; j<nSecondaries; ++j) delete theResult.GetSecondary(j)->GetParticle();
        theResult.Clear();
        theResult.SetStatusChange(stopAndKill);
        remnants.clear();
      } else {
        // 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 = (G4int)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 = (G4int)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, (*fragment)->GetCreatorModelID());
  }
      }
 
      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;
}

◆ DeleteModel()

void G4INCLXXInterface::DeleteModel ( )

inline

Definition at line 128 of file G4INCLXXInterface.hh.

                     {
    delete theINCLModel;
    theINCLModel = NULL;
  }

◆ GetDeExcitationModelName()

G4String const & G4INCLXXInterface::GetDeExcitationModelName ( ) const

Definition at line 678 of file G4INCLXXInterface.cc.

                                                                  {
  return theDeExcitation->GetModelName();
}

◆ ModelDescription()

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

virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 663 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";
}

◆ operator!=()

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

inline

Definition at line 112 of file G4INCLXXInterface.hh.

                                              {
    return (this != &right);
  }

◆ operator==()

G4bool 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 554 of file G4INCLXXInterface.cc.

                                                                                            {
  return 0;
}

◆ SetDeExcitation()

void G4VIntraNuclearTransportModel::SetDeExcitation ( G4VPreCompoundModel * ptr )

inline

Definition at line 81 of file G4VIntraNuclearTransportModel.hh.

{
  // previous pre-compound model will be deleted at the end of job 
  theDeExcitation = value;
}

Referenced by G4INCLXXInterfaceStore::UseAblaDeExcitation().

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

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4INCLXXInterface()

◆ ~G4INCLXXInterface()

Member Function Documentation

◆ ApplyYourself()

◆ DeleteModel()

◆ GetDeExcitationModelName()

◆ ModelDescription()

◆ operator!=()

◆ operator==()

◆ Propagate()

◆ SetDeExcitation()