#include <G4INCLCascade.hh>

Public Member Functions
	INCL (Config const *const config)
	~INCL ()
	INCL (const INCL &rhs)
	Dummy copy constructor to silence Coverity warning.
INCL &	operator= (const INCL &rhs)
	Dummy assignment operator to silence Coverity warning.
G4bool	prepareReaction (const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z, const G4int S)
G4bool	initializeTarget (const G4int A, const G4int Z, const G4int S, AnnihilationType theAType)
G4double	initUniverseRadiusForAntiprotonAtRest (const G4int A, const G4int Z)
const EventInfo &	processEvent ()
const EventInfo &	processEvent (ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4int targetA, const G4int targetZ, const G4int targetS)
void	finalizeGlobalInfo (Random::SeedVector const &initialSeeds)
const GlobalInfo &	getGlobalInfo () const

Detailed Description

Definition at line 56 of file G4INCLCascade.hh.

Constructor & Destructor Documentation

◆ INCL() [1/2]

G4INCL::INCL::INCL ( Config const *const config )

Definition at line 84 of file G4INCLCascade.cc.

    :propagationModel(0), theA(208), theZ(82), theS(0),
    targetInitSuccess(false),
    maxImpactParameter(0.),
    maxUniverseRadius(0.),
    maxInteractionDistance(0.),
    fixedImpactParameter(0.),
    theConfig(config),
    nucleus(NULL),
    forceTransparent(false),
    minRemnantSize(4)
  {
    // Set the logger object.
#ifdef INCLXX_IN_GEANT4_MODE
    Logger::initVerbosityLevelFromEnvvar();
#else // INCLXX_IN_GEANT4_MODE
    Logger::initialize(theConfig);
#endif // INCLXX_IN_GEANT4_MODE
 
    // Set the random number generator algorithm. The system can support
    // multiple different generator algorithms in a completely
    // transparent way.
    Random::initialize(theConfig);
 
    // Select the Pauli and CDPP blocking algorithms
    Pauli::initialize(theConfig);
 
    // Set the cross-section set
    CrossSections::initialize(theConfig);
 
    // Set the phase-space generator
    PhaseSpaceGenerator::initialize(theConfig);
 
    // Select the Coulomb-distortion algorithm:
    CoulombDistortion::initialize(theConfig);
 
    // Select the clustering algorithm:
    Clustering::initialize(theConfig);
 
    // Initialize the INCL particle table:
    ParticleTable::initialize(theConfig);
 
    // Initialize the value of cutNN in BinaryCollisionAvatar
    BinaryCollisionAvatar::setCutNN(theConfig->getCutNN());
 
    // Initialize the value of strange cross section bias
    BinaryCollisionAvatar::setBias(theConfig->getBias());
 
    // Propagation model is responsible for finding avatars and
    // transporting the particles. In principle this step is "hidden"
    // behind an abstract interface and the rest of the system does not
    // care how the transportation and avatar finding is done. This
    // should allow us to "easily" experiment with different avatar
    // finding schemes and even to support things like curved
    // trajectories in the future.
    propagationModel = new StandardPropagationModel(theConfig->getLocalEnergyBBType(),theConfig->getLocalEnergyPiType(),theConfig->getHadronizationTime());
    if(theConfig->getCascadeActionType() == AvatarDumpActionType)
      cascadeAction = new AvatarDumpAction();
    else
      cascadeAction = new CascadeAction();
    cascadeAction->beforeRunAction(theConfig);
 
    theGlobalInfo.cascadeModel = theConfig->getVersionString();
    theGlobalInfo.deexcitationModel = theConfig->getDeExcitationString();
#ifdef INCL_ROOT_USE
    theGlobalInfo.rootSelection = theConfig->getROOTSelectionString();
#endif
 
#ifndef INCLXX_IN_GEANT4_MODE
    // Fill in the global information
    theGlobalInfo.At = theConfig->getTargetA();
    theGlobalInfo.Zt = theConfig->getTargetZ();
    theGlobalInfo.St = theConfig->getTargetS();
    const ParticleSpecies theSpecies = theConfig->getProjectileSpecies();
    theGlobalInfo.Ap = theSpecies.theA;
    theGlobalInfo.Zp = theSpecies.theZ;
    theGlobalInfo.Sp = theSpecies.theS;
    theGlobalInfo.Ep = theConfig->getProjectileKineticEnergy();
    theGlobalInfo.biasFactor = theConfig->getBias();
#endif
 
    fixedImpactParameter = theConfig->getImpactParameter();
  }

Referenced by INCL(), and operator=().

◆ ~INCL()

G4INCL::INCL::~INCL ( )

Definition at line 168 of file G4INCLCascade.cc.

              {
    InteractionAvatar::deleteBackupParticles();
#ifndef INCLXX_IN_GEANT4_MODE
    NuclearMassTable::deleteTable();
#endif
    PhaseSpaceGenerator::deletePhaseSpaceGenerator();
    CrossSections::deleteCrossSections();
    Pauli::deleteBlockers();
    CoulombDistortion::deleteCoulomb();
    Random::deleteGenerator();
    Clustering::deleteClusteringModel();
#ifndef INCLXX_IN_GEANT4_MODE
    Logger::deleteLoggerSlave();
#endif
    NuclearDensityFactory::clearCache();
    NuclearPotential::clearCache();
    cascadeAction->afterRunAction();
    delete cascadeAction;
    delete propagationModel;
    delete theConfig;
  }

◆ INCL() [2/2]

G4INCL::INCL::INCL ( const INCL & rhs )

Dummy copy constructor to silence Coverity warning.

Member Function Documentation

◆ finalizeGlobalInfo()

void G4INCL::INCL::finalizeGlobalInfo ( Random::SeedVector const & initialSeeds )

Definition at line 1390 of file G4INCLCascade.cc.

                                                                    {
    const G4double normalisationFactor = theGlobalInfo.geometricCrossSection /
      ((G4double) theGlobalInfo.nShots);
    theGlobalInfo.nucleonAbsorptionCrossSection = normalisationFactor *
      ((G4double) theGlobalInfo.nNucleonAbsorptions);
    theGlobalInfo.pionAbsorptionCrossSection = normalisationFactor *
      ((G4double) theGlobalInfo.nPionAbsorptions);
    theGlobalInfo.reactionCrossSection = normalisationFactor *
      ((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
    theGlobalInfo.errorReactionCrossSection = normalisationFactor *
      std::sqrt((G4double) (theGlobalInfo.nShots - theGlobalInfo.nTransparents));
    theGlobalInfo.forcedCNCrossSection = normalisationFactor *
      ((G4double) theGlobalInfo.nForcedCompoundNucleus);
    theGlobalInfo.errorForcedCNCrossSection = normalisationFactor *
      std::sqrt((G4double) (theGlobalInfo.nForcedCompoundNucleus));
    theGlobalInfo.completeFusionCrossSection = normalisationFactor *
      ((G4double) theGlobalInfo.nCompleteFusion);
    theGlobalInfo.errorCompleteFusionCrossSection = normalisationFactor *
      std::sqrt((G4double) (theGlobalInfo.nCompleteFusion));
    theGlobalInfo.energyViolationInteractionCrossSection = normalisationFactor *
      ((G4double) theGlobalInfo.nEnergyViolationInteraction);
 
    theGlobalInfo.initialRandomSeeds.assign(initialSeeds.begin(), initialSeeds.end());
 
    Random::SeedVector theSeeds = Random::getSeeds();
    theGlobalInfo.finalRandomSeeds.assign(theSeeds.begin(), theSeeds.end());
  }

◆ getGlobalInfo()

const GlobalInfo & G4INCL::INCL::getGlobalInfo ( ) const

inline

Definition at line 90 of file G4INCLCascade.hh.

90{ return theGlobalInfo; }

◆ initializeTarget()

G4bool G4INCL::INCL::initializeTarget	(	const G4int	A,
		const G4int	Z,
		const G4int	S,
		AnnihilationType	theAType )

Definition at line 395 of file G4INCLCascade.cc.

                                                                                                      { 
    delete nucleus;
 
    if (theAType==PType || theAType==NType || theAType==DNbarNPbarPType || theAType==DNbarNPbarNType ||theAType==DNbarPPbarPType || theAType==DNbarPPbarNType) {
      G4double newmaxUniverseRadius=0.;
      if (theAType==PType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+1, Z+1);
      else if (theAType==NType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+1, Z);
      else if (theAType==DNbarPPbarPType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+2, Z+2);
      else if (theAType==DNbarNPbarNType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+2, Z);
      else if (theAType==DNbarNPbarPType || theAType==DNbarPPbarNType) newmaxUniverseRadius=initUniverseRadiusForAntiprotonAtRest(A+2, Z+1);
      nucleus = new Nucleus(A, Z, S, theConfig, newmaxUniverseRadius, theAType);
    }
    else{
      nucleus = new Nucleus(A, Z, S, theConfig, maxUniverseRadius, theAType);
    }
    nucleus->getStore()->getBook().reset();
    nucleus->initializeParticles();
    propagationModel->setNucleus(nucleus);
    return true;
  }

Referenced by prepareReaction().

◆ initUniverseRadiusForAntiprotonAtRest()

G4double G4INCL::INCL::initUniverseRadiusForAntiprotonAtRest	(	const G4int	A,
		const G4int	Z )

Definition at line 1534 of file G4INCLCascade.cc.

                                                                                   {
    G4double rMax = 0.0;
    if(A==0) {
      IsotopicDistribution const &anIsotopicDistribution =
        ParticleTable::getNaturalIsotopicDistribution(Z);
      IsotopeVector theIsotopes = anIsotopicDistribution.getIsotopes();
      for(IsotopeIter i=theIsotopes.begin(), e=theIsotopes.end(); i!=e; ++i) {
        const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, i->theA, Z);
        const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, i->theA, Z);
        const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
        rMax = std::max(maximumRadius, rMax);
      }
    } else {
      const G4double pMaximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
      const G4double nMaximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
      const G4double maximumRadius = std::max(pMaximumRadius, nMaximumRadius);
      rMax = std::max(maximumRadius, rMax);
    }
    return rMax;                
    }

Referenced by initializeTarget().

◆ operator=()

INCL & G4INCL::INCL::operator= ( const INCL & rhs )

Dummy assignment operator to silence Coverity warning.

◆ prepareReaction()

G4bool G4INCL::INCL::prepareReaction	(	const ParticleSpecies &	projectileSpecies,
		const G4double	kineticEnergy,
		const G4int	A,
		const G4int	Z,
		const G4int	S )

Definition at line 190 of file G4INCLCascade.cc.

                                                                                                                                                  {
    if(A < 0 || A > 300 || Z < 1 || Z > 200) {
      INCL_ERROR("Unsupported target: A = " << A << " Z = " << Z << " S = " << S << '\n'
                 << "Target configuration rejected." << '\n');
      return false;
    }
    if((projectileSpecies.theType==Composite || projectileSpecies.theType == antiComposite)&&
       (projectileSpecies.theZ==projectileSpecies.theA || projectileSpecies.theZ==0)) {
      INCL_ERROR("Unsupported projectile: A = " << projectileSpecies.theA << " Z = " << projectileSpecies.theZ << " S = " << projectileSpecies.theS << '\n'
                 << "Projectile configuration rejected." << '\n');
      return false;
    }
 
    // Reset the forced-transparent flag
    forceTransparent = false;
 
    // Initialise the maximum universe radius
    initUniverseRadius(projectileSpecies, kineticEnergy, A, Z);
    // Initialise the nucleus
 
//D
    //reset
    G4bool ProtonIsTheVictim = false; 
    G4bool NeutronIsTheVictim = false;
    G4bool DNbProtonIsTheVictim = false; 
    G4bool DPbProtonIsTheVictim = false;
    theEventInfo.annihilationP = false;
    theEventInfo.annihilationN = false;
    G4bool isModelA = true; //Antideuteron 
 
    //G4double AnnihilationBarrier = kineticEnergy;
    if((projectileSpecies.theType == antiProton && kineticEnergy <= theConfig->getAtrestThreshold()) || (projectileSpecies.theType == antiNeutron && kineticEnergy <= theConfig->getnbAtrestThreshold())){
      double SpOverSn;
      if(projectileSpecies.theType == antiProton)
        SpOverSn = 1.331;//from experiments with deuteron (E.Klempt)
      else if(projectileSpecies.theType == antiNeutron)
        SpOverSn = 1./1.331; //Opposite for antineutron
      else{
        SpOverSn = 1;
        INCL_ERROR("Neither antiProton nor antiNeutron annihilated");
      }
      //INCL_WARN("theA number set to A-1 from " << A <<'\n');
 
      G4double neutronprob;
      if(theConfig->isNaturalTarget()){ // A = 0 in this case
        theA = ParticleTable::drawRandomNaturalIsotope(Z) - 1; //43 and 61 are ok (Technetium and Promethium)
        neutronprob = (theA + 1 - Z)/(theA + 1 - Z + SpOverSn*Z);  
      }
      else{
        theA = A - 1;
        neutronprob = (A - Z)/(A - Z + SpOverSn*Z);  //from experiments with deuteron (E.Klempt)
      }
 
      theS = S;
      
      G4double rndm = Random::shoot();
      if(rndm >= neutronprob){     //proton is annihilated
        theEventInfo.annihilationP = true;
        theZ = Z - 1;
        ProtonIsTheVictim = true;
        //INCL_WARN("theZ number set to Z-1 from " << Z << '\n');
      }  
      else{        //neutron is annihilated
        theEventInfo.annihilationN = true;
        theZ = Z;
        NeutronIsTheVictim = true;
      }  
    } else if(projectileSpecies.theType == antiComposite && kineticEnergy <= theConfig->getdbAtrestThreshold()){ 
      if(Z > 30)
        isModelA=false;
      else if(Z > 9 && Z <=30){ //Maybe change and add another dependance than Z
        double rndmA = Random::shoot();
        if(rndmA > 0.5)//Random threshold : should be improved to take into account the orbit in which the separation takes place.
          isModelA=false;
      }
      if(isModelA){
        //Antideuteron Model A case : 2 annihilation at the same time
        double pbarSpOverSn = 1.331;
        double nbarSpOverSn = 1./1.331;
        double pbarneutronprob;
        double nbarneutronprob;
        if(theConfig->isNaturalTarget()){
          theA = ParticleTable::drawRandomNaturalIsotope(Z) - 2;
          nbarneutronprob = (theA + 2 - Z)/(theA + 2 - Z + nbarSpOverSn*Z);  
          pbarneutronprob = (theA + 2 - Z)/(theA + 2 - Z + pbarSpOverSn*Z);  
        }
        else{
          theA = A - 2;
          nbarneutronprob = (A - Z)/(A - Z + nbarSpOverSn*Z); 
          pbarneutronprob = (A - Z)/(A - Z + pbarSpOverSn*Z);
    }
        theS = S;
        G4double rndm = Random::shoot(); //for nbar
        G4double rndm2 = Random::shoot(); //for pbar
 
        if (rndm >= nbarneutronprob){ // nbarp
          DNbProtonIsTheVictim = true;
          if(rndm2 >= pbarneutronprob){ // pbarp
            theZ = Z - 2;
            DPbProtonIsTheVictim = true;
          } else if(rndm2 < pbarneutronprob){ //pbarn
            theZ = Z - 1;
          }
        } else if(rndm < nbarneutronprob){//nbarn
          if(rndm2 >= pbarneutronprob){ // pbarp
            theZ = Z - 1;
            DPbProtonIsTheVictim = true;
          } else if(rndm2 < pbarneutronprob){ //pbarn
            theZ = Z;
          }
        }
      } else if (!isModelA){//Model B : Antiproton is detached from antideuteron
        double SpOverSn = 1.331;
        double neutronprob;
        if(theConfig->isNaturalTarget()){
          theA = ParticleTable::drawRandomNaturalIsotope(Z) - 1; 
          neutronprob = (theA + 1 - Z)/(theA + 1 - Z + SpOverSn*Z);  
        }
        else{
          theA = A - 1;
          neutronprob = (A - Z)/(A - Z + SpOverSn*Z);
        }
 
        theS = S;
      
        double rndm = Random::shoot();
        if(rndm >= neutronprob){     //proton is annihilated
          theEventInfo.annihilationP = true;
          theZ = Z - 1;
          ProtonIsTheVictim = true;
        }  
        else{        //neutron is annihilated
          theEventInfo.annihilationN = true;
          theZ = Z;
          NeutronIsTheVictim = true;
        }  
      }
    }
    else{ // not annihilation of pbar, nbar, dbar
      theZ = Z;
      theS = S;
      if(theConfig->isNaturalTarget())
        theA = ParticleTable::drawRandomNaturalIsotope(Z); //change order
      else
        theA = A;
    }
 
    AnnihilationType theAType = Def;
    if(ProtonIsTheVictim == true && NeutronIsTheVictim == false)
    theAType = PType;
    if(NeutronIsTheVictim == true && ProtonIsTheVictim == false)
    theAType = NType;
    if(projectileSpecies.theType == antiComposite && kineticEnergy <= theConfig->getdbAtrestThreshold() && isModelA){
      if(DNbProtonIsTheVictim == true && DPbProtonIsTheVictim ==true)
        theAType = DNbarPPbarPType;
      else if(DNbProtonIsTheVictim == false && DPbProtonIsTheVictim ==true)
        theAType = DNbarNPbarPType;
      else if(DNbProtonIsTheVictim == false && DPbProtonIsTheVictim == false)
        theAType = DNbarNPbarNType;
      else if(DNbProtonIsTheVictim == true && DPbProtonIsTheVictim ==false)
        theAType = DNbarPPbarNType;
    } else if (projectileSpecies.theType == antiComposite && kineticEnergy <= theConfig->getdbAtrestThreshold() && !isModelA){
      if(ProtonIsTheVictim == true && NeutronIsTheVictim == false)
        theAType = PType;
      if(NeutronIsTheVictim == true && ProtonIsTheVictim == false)
        theAType = NType;
    }
 
//D
 
    initializeTarget(theA, theZ, theS, theAType);
    
    // Set the maximum impact parameter
    maxImpactParameter = CoulombDistortion::maxImpactParameter(projectileSpecies, kineticEnergy, nucleus);
    INCL_DEBUG("Maximum impact parameter initialised: " << maxImpactParameter << '\n');
 
    // For forced CN events
    initMaxInteractionDistance(projectileSpecies, kineticEnergy);
// Set the geometric cross sectiony section
    if((projectileSpecies.theType == antiProton && kineticEnergy <= theConfig->getAtrestThreshold()) || (projectileSpecies.theType == antiNeutron && kineticEnergy <= theConfig->getnbAtrestThreshold())
      || (projectileSpecies.theType == antiComposite && kineticEnergy <= theConfig->getdbAtrestThreshold()) ){
      G4int currentA = A;
      if(theConfig->isNaturalTarget()){
        currentA = ParticleTable::drawRandomNaturalIsotope(Z);
      }
      G4double kineticEnergy2=kineticEnergy;
      if (kineticEnergy2 <= 0.) kineticEnergy2=0.001;
      theGlobalInfo.geometricCrossSection = 9.7* //normalization factor from Corradini
        Math::pi*std::pow((1.840 + 1.120*std::pow(currentA,(1./3.))),2)*
        (1. + (Z*G4INCL::PhysicalConstants::eSquared*(currentA+1))/(currentA*kineticEnergy2*(1.840 + 1.120*std::pow(currentA,(1./3.))))); 
         //xsection formula was borrowed from Corradini et al. https://doi.org/10.1016/j.physletb.2011.09.069    
    }
    else{
      theGlobalInfo.geometricCrossSection =
        Math::tenPi*std::pow(maxImpactParameter,2);
    }
 
    // Set the minimum remnant size
    if(projectileSpecies.theA > 0)
      minRemnantSize = std::min(theA, 4);
    else
      minRemnantSize = std::min(theA-1, 4);
    return true;
  }

Referenced by processEvent().

◆ processEvent() [1/2]

const EventInfo & G4INCL::INCL::processEvent ( )