#include <G4HEAntiNeutronInelastic.hh>

Inheritance diagram for G4HEAntiNeutronInelastic:

Public Member Functions
	G4HEAntiNeutronInelastic (const G4String &name="G4HEAntiNeutronInelastic")

	~G4HEAntiNeutronInelastic ()

virtual void	ModelDescription (std::ostream &) const

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

G4int	GetNumberOfSecondaries ()

void	FirstIntInCasAntiNeutron (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)

Public Member Functions inherited from G4HEInelastic
	G4HEInelastic (const G4String &modelName="HEInelastic")

	~G4HEInelastic ()

void	SetMaxNumberOfSecondaries (const G4int maxnumber)

void	SetVerboseLevel (const G4int level)

void	ForceEnergyConservation (G4bool energyConservation)

G4bool	EnergyConservation (void)

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

G4double	Amin (G4double a, G4double b)

G4double	Amax (G4double a, G4double b)

G4int	Imin (G4int a, G4int b)

G4int	Imax (G4int a, G4int b)

void	FillParticleChange (G4HEVector pv[], G4int aVecLength)

G4double	pmltpc (G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)

G4int	Factorial (G4int n)

G4double	NuclearInelasticity (G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber)

G4double	NuclearExcitation (G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber, G4double &excitationEnergyCascade, G4double &excitationEnergyEvaporation)

void	HighEnergyCascading (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	HighEnergyClusterProduction (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	TuningOfHighEnergyCascading (G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	MediumEnergyCascading (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	MediumEnergyClusterProduction (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	QuasiElasticScattering (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)

void	ElasticScattering (G4bool &successful, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, G4double atomicWeight, G4double atomicNumber)

G4int	rtmi (G4double *x, G4double xli, G4double xri, G4double eps, G4int iend, G4double aa, G4double bb, G4double cc, G4double dd, G4double rr)

G4double	fctcos (G4double t, G4double aa, G4double bb, G4double cc, G4double dd, G4double rr)

void	StrangeParticlePairProduction (const G4double availableEnergy, const G4double centerOfMassEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)

G4double	NBodyPhaseSpace (const G4double totalEnergy, const G4bool constantCrossSection, G4HEVector pv[], G4int &vecLen)

G4double	NBodyPhaseSpace (G4int npart, G4HEVector pv[], G4double wmax, G4double wfcn, G4int maxtrial, G4int ntrial)

G4double	gpdk (G4double a, G4double b, G4double c)

void	QuickSort (G4double arr[], const G4int lidx, const G4int ridx)

G4double	Alam (G4double a, G4double b, G4double c)

G4double	CalculatePhaseSpaceWeight (G4int npart)

G4double	normal (void)

G4double	GammaRand (G4double avalue)

G4double	Erlang (G4int mvalue)

G4int	Poisson (G4double x)

void	SetParticles (void)

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

Public Attributes
G4int	vecLength

Public Attributes inherited from G4HEInelastic
G4int	verboseLevel

G4int	MAXPART

G4bool	conserveEnergy

G4HEVector	PionPlus

G4HEVector	PionZero

G4HEVector	PionMinus

G4HEVector	KaonPlus

G4HEVector	KaonZero

G4HEVector	AntiKaonZero

G4HEVector	KaonMinus

G4HEVector	KaonZeroShort

G4HEVector	KaonZeroLong

G4HEVector	Proton

G4HEVector	AntiProton

G4HEVector	Neutron

G4HEVector	AntiNeutron

G4HEVector	Lambda

G4HEVector	AntiLambda

G4HEVector	SigmaPlus

G4HEVector	SigmaZero

G4HEVector	SigmaMinus

G4HEVector	AntiSigmaPlus

G4HEVector	AntiSigmaZero

G4HEVector	AntiSigmaMinus

G4HEVector	XiZero

G4HEVector	XiMinus

G4HEVector	AntiXiZero

G4HEVector	AntiXiMinus

G4HEVector	OmegaMinus

G4HEVector	AntiOmegaMinus

G4HEVector	Deuteron

G4HEVector	Triton

G4HEVector	Alpha

G4HEVector	Gamma

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

G4bool	IsBlocked () const

void	Block ()

Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange

G4int	verboseLevel

G4double	theMinEnergy

G4double	theMaxEnergy

G4bool	isBlocked

Detailed Description

Definition at line 50 of file G4HEAntiNeutronInelastic.hh.

Constructor & Destructor Documentation

◆ G4HEAntiNeutronInelastic()

G4HEAntiNeutronInelastic::G4HEAntiNeutronInelastic ( const G4String & name = "G4HEAntiNeutronInelastic" )

Definition at line 43 of file G4HEAntiNeutronInelastic.cc.

 : G4HEInelastic(name)
{
  vecLength = 0;
  theMinEnergy = 20*GeV;
  theMaxEnergy = 10*TeV;
  MAXPART      = 2048;
  verboseLevel = 0;
  G4cout << "WARNING: model G4HEAntiNeutronInelastic is being deprecated and will\n"
         << "disappear in Geant4 version 10.0"  << G4endl;  
}

◆ ~G4HEAntiNeutronInelastic()

G4HEAntiNeutronInelastic::~G4HEAntiNeutronInelastic ( )

inline

Definition at line 54 of file G4HEAntiNeutronInelastic.hh.

54{};

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4HEAntiNeutronInelastic::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus
	)

virtual

Implements G4HadronicInteraction.

Definition at line 71 of file G4HEAntiNeutronInelastic.cc.

{
  G4HEVector* pv = new G4HEVector[MAXPART];
  const G4HadProjectile* aParticle = &aTrack;
  const G4double atomicWeight = targetNucleus.GetA_asInt();
  const G4double atomicNumber = targetNucleus.GetZ_asInt();
  G4HEVector incidentParticle(aParticle);
 
  G4int incidentCode = incidentParticle.getCode();
  G4double incidentMass = incidentParticle.getMass();
  G4double incidentTotalEnergy = incidentParticle.getEnergy();
 
  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
  // DHW 19 May 2011: variable set but not used
 
  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
 
  if (incidentKineticEnergy < 1.)
    G4cout << "GHEAntiNeutronInelastic: incident energy < 1 GeV" << G4endl;;
 
  if (verboseLevel > 1) {
    G4cout << "G4HEAntiNeutronInelastic::ApplyYourself" << G4endl;
    G4cout << "incident particle " << incidentParticle.getName()
           << "mass "              << incidentMass
           << "kinetic energy "    << incidentKineticEnergy
           << G4endl;
    G4cout << "target material with (A,Z) = (" 
           << atomicWeight << "," << atomicNumber << ")" << G4endl;
  }
 
  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
                                              atomicWeight, atomicNumber);
  if (verboseLevel > 1)
    G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
    
  incidentKineticEnergy -= inelasticity;
    
  G4double excitationEnergyGNP = 0.;
  G4double excitationEnergyDTA = 0.; 
 
  G4double excitation = NuclearExcitation(incidentKineticEnergy,
                                          atomicWeight, atomicNumber,
                                          excitationEnergyGNP,
                                          excitationEnergyDTA);
  if (verboseLevel > 1)
    G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
           << excitationEnergyDTA << G4endl;             
 
  incidentKineticEnergy -= excitation;
  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
  //                                   *(incidentTotalEnergy+incidentMass));
  //  DHW 19 May 2011: variable set but not used
 
  G4HEVector targetParticle;
  if (G4UniformRand() < atomicNumber/atomicWeight) { 
    targetParticle.setDefinition("Proton");
  } else { 
    targetParticle.setDefinition("Neutron");
  }
 
  G4double targetMass = targetParticle.getMass();
  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
                                        + targetMass*targetMass
                                        + 2.0*targetMass*incidentTotalEnergy);
  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
 
  G4bool inElastic = true;
  vecLength = 0;           
        
  if (verboseLevel > 1)
    G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
           << incidentCode << G4endl;
 
  G4bool successful = false; 
    
  FirstIntInCasAntiNeutron(inElastic, availableEnergy, pv, vecLength,
                           incidentParticle, targetParticle, atomicWeight);
 
  if (verboseLevel > 1)
    G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
 
  if ((vecLength > 0) && (availableEnergy > 1.)) 
    StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
                                  pv, vecLength,
                                  incidentParticle, targetParticle);
  HighEnergyCascading(successful, pv, vecLength,
                      excitationEnergyGNP, excitationEnergyDTA,
                      incidentParticle, targetParticle,
                      atomicWeight, atomicNumber);
  if (!successful)
    HighEnergyClusterProduction(successful, pv, vecLength,
                                excitationEnergyGNP, excitationEnergyDTA,
                                incidentParticle, targetParticle,
                                atomicWeight, atomicNumber);
  if (!successful) 
    MediumEnergyCascading(successful, pv, vecLength, 
                          excitationEnergyGNP, excitationEnergyDTA, 
                          incidentParticle, targetParticle,
                          atomicWeight, atomicNumber);
 
  if (!successful)
    MediumEnergyClusterProduction(successful, pv, vecLength,
                                  excitationEnergyGNP, excitationEnergyDTA,       
                                  incidentParticle, targetParticle,
                                  atomicWeight, atomicNumber);
  if (!successful)
    QuasiElasticScattering(successful, pv, vecLength,
                           excitationEnergyGNP, excitationEnergyDTA,
                           incidentParticle, targetParticle, 
                           atomicWeight, atomicNumber);
  if (!successful)
    ElasticScattering(successful, pv, vecLength,
                      incidentParticle,    
                      atomicWeight, atomicNumber);
 
  if (!successful)
    G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
           << G4endl;
 
  FillParticleChange(pv,  vecLength);
  delete [] pv;
  theParticleChange.SetStatusChange(stopAndKill);
  return &theParticleChange;
}

◆ FirstIntInCasAntiNeutron()

void G4HEAntiNeutronInelastic::FirstIntInCasAntiNeutron	(	G4bool &	inElastic,
		const G4double	availableEnergy,
		G4HEVector	pv[],
		G4int &	vecLen,
		const G4HEVector &	incidentParticle,
		const G4HEVector &	targetParticle,
		const G4double	atomicWeight
	)

Definition at line 200 of file G4HEAntiNeutronInelastic.cc.

{
  static const G4double expxu = 82.;     // upper bound for arg. of exp
  static const G4double expxl = -expxu;  // lower bound for arg. of exp
 
  static const G4double protb = 0.7;
  static const G4double neutb = 0.7;
  static const G4double     c = 1.25;
 
  static const G4int numMul   = 1200;
  static const G4int numMulAn = 400;
  static const G4int numSec   = 60;
 
  G4int neutronCode = Neutron.getCode();
  G4int protonCode = Proton.getCode();
 
  G4int targetCode = targetParticle.getCode();
  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
 
  static G4bool first = true;
  static G4double protmul[numMul],  protnorm[numSec];   // proton constants
  static G4double protmulAn[numMulAn],protnormAn[numSec]; 
  static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
 
  //  misc. local variables
  //  npos = number of pi+,  nneg = number of pi-,  nzero = number of pi0
 
  G4int i, counter, nt, npos, nneg, nzero;
 
   if( first ) 
     {                 // compute normalization constants, this will only be done once
       first = false;
       for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
       for( i=0; i<numSec  ; i++ ) protnorm[i] = 0.0;
       counter = -1;
       for( npos=0; npos<(numSec/3); npos++ ) 
          {
            for( nneg=std::max(0,npos-2); nneg<=npos; nneg++ ) 
               {
                 for( nzero=0; nzero<numSec/3; nzero++ ) 
                    {
                      if( ++counter < numMul ) 
                        {
                          nt = npos+nneg+nzero;
                          if( (nt>0) && (nt<=numSec) ) 
                            {
                              protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
                              protnorm[nt-1] += protmul[counter];
                            }
                        }
                    }
               }
          }
       for( i=0; i<numMul; i++ )neutmul[i]  = 0.0;
       for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
       counter = -1;
       for( npos=0; npos<numSec/3; npos++ ) 
          {
            for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 
               {
                 for( nzero=0; nzero<numSec/3; nzero++ ) 
                    {
                      if( ++counter < numMul ) 
                        {
                          nt = npos+nneg+nzero;
                          if( (nt>0) && (nt<=numSec) ) 
                            {
                               neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
                               neutnorm[nt-1] += neutmul[counter];
                            }
                        }
                    }
               }
          }
       for( i=0; i<numSec; i++ ) 
          {
            if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
            if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
          }
                                                                   // annihilation
       for( i=0; i<numMulAn  ; i++ ) protmulAn[i]  = 0.0;
       for( i=0; i<numSec    ; i++ ) protnormAn[i] = 0.0;
       counter = -1;
       for( npos=1; npos<(numSec/3); npos++ ) 
          {
            nneg = std::max(0,npos-1); 
            for( nzero=0; nzero<numSec/3; nzero++ ) 
               {
                 if( ++counter < numMulAn ) 
                   {
                     nt = npos+nneg+nzero;
                     if( (nt>1) && (nt<=numSec) ) 
                       {
                         protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
                         protnormAn[nt-1] += protmulAn[counter];
                       }
                   }
               }
          }
       for( i=0; i<numMulAn; i++ ) neutmulAn[i]  = 0.0;
       for( i=0; i<numSec;   i++ ) neutnormAn[i] = 0.0;
       counter = -1;
       for( npos=0; npos<numSec/3; npos++ ) 
          {
            nneg = npos; 
            for( nzero=0; nzero<numSec/3; nzero++ ) 
               {
                 if( ++counter < numMulAn ) 
                   {
                     nt = npos+nneg+nzero;
                     if( (nt>1) && (nt<=numSec) ) 
                       {
                          neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
                          neutnormAn[nt-1] += neutmulAn[counter];
                       }
                   }
               }
          }
       for( i=0; i<numSec; i++ ) 
          {
            if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
            if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
          }
     }                                          // end of initialization
 
         
                                              // initialize the first two places
                                              // the same as beam and target                                    
   pv[0] = incidentParticle;
   pv[1] = targetParticle;
   vecLen = 2;
 
   if( !inElastic ) 
     {                                        // nb n --> pb p
       if( targetCode == neutronCode ) 
         {
           G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
 
           G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
           if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) ) 
             {
               pv[0] = AntiProton;
               pv[1] = Proton;
             }
         }
       return;
     }
   else if (availableEnergy <= PionPlus.getMass())
       return;
 
                                                  //   inelastic scattering
 
   npos = 0, nneg = 0, nzero = 0;
   G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88, 
                      0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40, 
                      0.39, 0.36, 0.33, 0.10, 0.01};
   G4int            iplab =      G4int( incidentTotalMomentum*10.);
   if ( iplab >  9) iplab = 10 + G4int( (incidentTotalMomentum  -1.)*5. );          
   if ( iplab > 14) iplab = 15 + G4int(  incidentTotalMomentum  -2.     );
   if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.); 
                    iplab = std::min(24, iplab);
 
   if ( G4UniformRand() > anhl[iplab] )
     {
 
       G4double eab = availableEnergy;
       G4int ieab = G4int( eab*5.0 );
   
       G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
       if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) ) 
         {
                                              //   suppress high multiplicity events at low momentum
                                              //   only one additional pion will be produced
           G4double w0, wp, wm, wt, ran;
           if( targetCode == protonCode )                    // target is a proton 
             {
               w0 = - sqr(1.+protb)/(2.*c*c);
               w0 = wp = std::exp(w0);
               if( G4UniformRand() < w0/(w0+wp) )  
                 { npos = 0; nneg = 0; nzero = 1; }
               else 
                 { npos = 1; nneg = 0; nzero = 0; }       
             } 
           else 
             {                                               // target is a neutron
               w0 = -sqr(1.+neutb)/(2.*c*c);
               w0 = wp = std::exp(w0);
               wm = -sqr(-1.+neutb)/(2.*c*c);
               wm = std::exp(wm);
               wt = w0+wp+wm;
               wp = w0+wp;
               ran = G4UniformRand();
               if( ran < w0/wt)
                 { npos = 0; nneg = 0; nzero = 1; }       
               else if( ran < wp/wt)
                 { npos = 1; nneg = 0; nzero = 0; }       
               else
                 { npos = 0; nneg = 1; nzero = 0; }       
             }
         }
       else
         {
                         //  number of total particles vs. centre of mass Energy - 2*proton mass
   
           G4double aleab = std::log(availableEnergy);
           G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
                            + aleab*(0.117712+0.0136912*aleab))) - 2.0;
   
                         // normalization constant for kno-distribution.
                         // calculate first the sum of all constants, check for numerical problems.   
           G4double test, dum, anpn = 0.0;
 
           for (nt=1; nt<=numSec; nt++) {
               test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
               dum = pi*nt/(2.0*n*n);
 
               if (std::fabs(dum) < 1.0) { 
                 if( test >= 1.0e-10 )anpn += dum*test;
               } else { 
                 anpn += dum*test;
           }
       }
   
           G4double ran = G4UniformRand();
           G4double excs = 0.0;
           if (targetCode == protonCode) {
             counter = -1;
             for (npos=0; npos<numSec/3; npos++) {
               for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) {
                 for (nzero=0; nzero<numSec/3; nzero++) {
                   if (++counter < numMul) {
                     nt = npos+nneg+nzero;
                     if ( (nt>0) && (nt<=numSec) ) {
                       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
                       dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
 
                       if (std::fabs(dum) < 1.0) { 
                         if( test >= 1.0e-10 )excs += dum*test;
                       } else { 
                         excs += dum*test;
                       }
 
                       if (ran < excs) goto outOfLoop;      //----------------------->
             }  
           }    
         }     
               }                                                                                  
         }
       
                                           // 3 previous loops continued to the end
             inElastic = false;            // quasi-elastic scattering   
             return;
 
       } else {                        // target must be a neutron
             counter = -1;
             for (npos=0; npos<numSec/3; npos++) {
               for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
                 for (nzero=0; nzero<numSec/3; nzero++) {
                   if (++counter < numMul) {
                     nt = npos+nneg+nzero;
                     if ((nt>0) && (nt<=numSec) ) {
                       test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
                       dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
                       if (std::fabs(dum) < 1.0) { 
                         if( test >= 1.0e-10 )excs += dum*test;
                       } else { 
                         excs += dum*test;
                       }
 
                       if (ran < excs) goto outOfLoop;       // -------------------------->
                     }
                   }
                 }
               }
         }
                                                      // 3 previous loops continued to the end
             inElastic = false;                     // quasi-elastic scattering.
             return;
           }
         } 
       outOfLoop:           //  <------------------------------------------------------------------------   
    
       if( targetCode == protonCode)
         {
           if( npos == nneg)
             {
             }
           else if (npos == (nneg+1))
             {
               if( G4UniformRand() < 0.5)
                 {
                   pv[1] = Neutron;
                 }
               else
                 {
                   pv[0] = AntiProton;
                 }
             }
           else      
             {
               pv[0] = AntiProton;
               pv[1] = Neutron;
             } 
         }  
       else
         {
           if( npos == nneg)
             {
               if( G4UniformRand() < 0.25)
                 {
                   pv[0] = AntiProton;
                   pv[1] = Proton;
                 }
               else
                 {
                 }
             } 
           else if ( npos == (nneg-1))
             {
               pv[1] = Proton;
             }
           else
             {
               pv[0] = AntiProton;
             }
         }      
     
     }
   else                                                               // annihilation
     {  
       if ( availableEnergy > 2. * PionPlus.getMass() )
         {
 
           G4double aleab = std::log(availableEnergy);
           G4double n     = 3.62567+aleab*(0.665843+aleab*(0.336514
                            + aleab*(0.117712+0.0136912*aleab))) - 2.0;
   
                      //   normalization constant for kno-distribution.
                      //   calculate first the sum of all constants, check for numerical problems.   
           G4double test, dum, anpn = 0.0;
 
           for (nt=2; nt<=numSec; nt++) {
             test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
             dum = pi*nt/(2.0*n*n);
             if (std::fabs(dum) < 1.0) { 
               if( test >= 1.0e-10 )anpn += dum*test;
             } else { 
               anpn += dum*test;
         }
       }
   
           G4double ran = G4UniformRand();
           G4double excs = 0.0;
           if (targetCode == protonCode) {
             counter = -1;
             for (npos=1; npos<numSec/3; npos++) {
               nneg = npos-1; 
               for (nzero=0; nzero<numSec/3; nzero++) {
                 if (++counter < numMulAn) {
                   nt = npos+nneg+nzero;
                   if ( (nt>1) && (nt<=numSec) ) {
                     test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
                     dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
 
                     if (std::fabs(dum) < 1.0) { 
                       if( test >= 1.0e-10 )excs += dum*test;
                     } else { 
                       excs += dum*test;
             }
 
                     if (ran < excs) goto outOfLoopAn;      //----------------------->
           }
         }
               }
             }
                                             // 3 previous loops continued to the end
             inElastic = false;            // quasi-elastic scattering   
             return;
 
           } else {                          // target must be a neutron
             counter = -1;
             for (npos=0; npos<numSec/3; npos++) { 
               nneg = npos; 
               for (nzero=0; nzero<numSec/3; nzero++) {
                 if (++counter < numMulAn) {
                   nt = npos+nneg+nzero;
                   if ( (nt>1) && (nt<=numSec) ) {
                     test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
                     dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
 
                     if (std::fabs(dum) < 1.0) { 
                       if( test >= 1.0e-10 )excs += dum*test;
                     } else { 
                       excs += dum*test;
             }
 
                     if (ran < excs) goto outOfLoopAn;       // -------------------------->
           }
         }
           }
         }
             inElastic = false;                     // quasi-elastic scattering.
             return;
       }
       outOfLoopAn:           //  <------------------------------------------------------------------   
       vecLen = 0;
         }
     }
 
   nt = npos + nneg + nzero;
   while ( nt > 0)
       {
         G4double ran = G4UniformRand();
         if ( ran < (G4double)npos/nt)
            { 
              if( npos > 0 ) 
                { pv[vecLen++] = PionPlus;
                  npos--;
                }
            }
         else if ( ran < (G4double)(npos+nneg)/nt)
            {   
              if( nneg > 0 )
                { 
                  pv[vecLen++] = PionMinus;
                  nneg--;
                }
            }
         else
            {
              if( nzero > 0 )
                { 
                  pv[vecLen++] = PionZero;
                  nzero--;
                }
            }
         nt = npos + nneg + nzero;
       } 
   if (verboseLevel > 1)
      {
        G4cout << "Particles produced: " ;
        G4cout << pv[0].getName() << " " ;
        G4cout << pv[1].getName() << " " ;
        for (i=2; i < vecLen; i++)   
            { 
              G4cout << pv[i].getName() << " " ;
            }
         G4cout << G4endl;
      }
   return;
 }

Referenced by ApplyYourself().

◆ GetNumberOfSecondaries()

G4int G4HEAntiNeutronInelastic::GetNumberOfSecondaries ( )

inline

Definition at line 63 of file G4HEAntiNeutronInelastic.hh.

63{return vecLength;}

◆ ModelDescription()

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

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 56 of file G4HEAntiNeutronInelastic.cc.

{
  outFile << "G4HEAntiNeutronInelastic is one of the High Energy\n"
          << "Parameterized (HEP) models used to implement inelastic\n"
          << "anti-neutron scattering from nuclei.  It is a re-engineered\n"
          << "version of the GHEISHA code of H. Fesefeldt.  It divides the\n"
          << "initial collision products into backward- and forward-going\n"
          << "clusters which are then decayed into final state hadrons.\n"
          << "The model does not conserve energy on an event-by-event\n"
          << "basis.  It may be applied to anti-neutrons with initial energies\n"
          << "above 20 GeV.\n";
}

Member Data Documentation

◆ vecLength

G4int G4HEAntiNeutronInelastic::vecLength

Definition at line 58 of file G4HEAntiNeutronInelastic.hh.

Referenced by ApplyYourself(), G4HEAntiNeutronInelastic(), and GetNumberOfSecondaries().

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

Public Member Functions

Public Attributes

Additional Inherited Members