INCL++ intra-nuclear cascade with G4ExcitationHandler for de-excitation. More...

#include <G4INCLXXInterface.hh>

Inheritance diagram for G4INCLXXInterface:

Public Member Functions
	G4INCLXXInterface (const G4String &name="INCL++ cascade with G4ExcitationHandler")

	~G4INCLXXInterface ()

G4int	operator== (G4INCLXXInterface &right)

G4int	operator!= (G4INCLXXInterface &right)

G4ReactionProductVector *	Propagate (G4KineticTrackVector theSecondaries, G4V3DNucleus theNucleus)

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

void	DeleteModel ()

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

virtual	~G4VIntraNuclearTransportModel ()

virtual G4ReactionProductVector *	Propagate (G4KineticTrackVector theSecondaries, G4V3DNucleus theNucleus)=0

void	SetDeExcitation (G4VPreCompoundModel *ptr)

void	Set3DNucleus (G4V3DNucleus *const value)

void	SetPrimaryProjectile (const G4HadProjectile &aPrimary)

const G4String &	GetModelName () const

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

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

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

virtual	~G4HadronicInteraction ()

virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0

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

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

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	ModelDescription (std::ostream &outFile) const

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 with G4ExcitationHandler for de-excitation.

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 93 of file G4INCLXXInterface.hh.

Constructor & Destructor Documentation

◆ G4INCLXXInterface()

G4INCLXXInterface::G4INCLXXInterface ( const G4String & name = "INCL++ cascade with G4ExcitationHandler" )

Definition at line 51 of file G4INCLXXInterface.cc.

                                                        :
  G4VIntraNuclearTransportModel(nam),
  theINCLModel(NULL),
  theInterfaceStore(G4INCLXXInterfaceStore::GetInstance()),
  complainedAboutBackupModel(false)
{
  // 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);
    theExcitationHandler = 0;
  } else {
    theExcitationHandler = new G4ExcitationHandler;
  }
 
  theBackupModel = new G4BinaryLightIonReaction;
}

◆ ~G4INCLXXInterface()

G4INCLXXInterface::~G4INCLXXInterface ( )

Definition at line 69 of file G4INCLXXInterface.cc.

{
  delete theBackupModel;
  delete theExcitationHandler;
}

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

Implements G4HadronicInteraction.

Definition at line 120 of file G4INCLXXInterface.cc.

{
  // For systems heavier than theMaxProjMassINCL, use another model (typically
  // BIC)
  const G4int theMaxProjMassINCL = theInterfaceStore->GetMaxProjMassINCL();
  if(aTrack.GetDefinition()->GetAtomicMass() > theMaxProjMassINCL
      && theNucleus.GetA_asInt() > 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.";
      G4cout << "[INCL++] Warning: " << ss.str() << G4endl;
    }
    return theBackupModel->ApplyYourself(aTrack, theNucleus);
  }
 
  const G4int maxTries = 200;
 
  // 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) {
    G4ParticleTable * const theParticleTable = G4ParticleTable::GetParticleTable();
    const G4int oldTargetA = theNucleus.GetA_asInt();
    const G4int oldTargetZ = theNucleus.GetZ_asInt();
    G4ParticleDefinition *oldTargetDef = theParticleTable->GetIon(oldTargetZ, oldTargetA, 0.0);
    const G4ParticleDefinition *oldProjectileDef = aTrack.GetDefinition();
 
    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);
        const G4double theProjectileMass = theParticleTable->GetIonTable()->GetIonMass(oldTargetZ, oldTargetA);
        toInverseKinematics = new G4LorentzRotation(aTrack.Get4Momentum().boostVector());
        G4LorentzVector theProjectile4Momentum(0.0, 0.0, 0.0, theProjectileMass);
        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;
 
  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();
 
    if(theInterfaceStore->GetDumpInput()) {
      G4cout << theINCLModel->configToString() << G4endl;
    }
    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());
      }
 
      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);
      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 " << aTrack.GetDefinition()->GetParticleName() << " + "
            << G4ParticleTable::GetParticleTable()->GetIon(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0.0)->GetParticleName()
            << ", 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());
        }
 
    G4Fragment remnant(A, Z, fourMomentum);
        remnant.SetAngularMomentum(spin);
    remnants.push_back(remnant);
      }
    }
    nTries++;
  } while(!eventIsOK && nTries < maxTries);
 
  // 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 " << aTrack.GetDefinition()->GetParticleName()
    << " + " << G4ParticleTable::GetParticleTable()->GetIon(theNucleus.GetZ_asInt(), theNucleus.GetA_asInt(), 0.0)->GetParticleName()
    << " 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(theExcitationHandler != 0) {
    for(std::list<G4Fragment>::const_iterator i = remnants.begin();
    i != remnants.end(); i++) {
      G4ReactionProductVector *deExcitationResult = theExcitationHandler->BreakItUp((*i));
    
      for(G4ReactionProductVector::iterator fragment = deExcitationResult->begin();
      fragment != deExcitationResult->end(); ++fragment) {
    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();
 
  return &theResult;
}

◆ DeleteModel()

void G4INCLXXInterface::DeleteModel ( )

inline

Definition at line 116 of file G4INCLXXInterface.hh.

                     {
    delete theINCLModel;
    theINCLModel = NULL;
  }

◆ operator!=()

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

inline

Definition at line 102 of file G4INCLXXInterface.hh.

                                             {
    return (this != &right);
  }

◆ operator==()

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

inline

Definition at line 98 of file G4INCLXXInterface.hh.

                                             {
    return (this == &right);
  }

◆ Propagate()

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

virtual

Implements G4VIntraNuclearTransportModel.

Definition at line 371 of file G4INCLXXInterface.cc.

                                                                                            {
  return 0;
}

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

◆ operator!=()

◆ operator==()

◆ Propagate()