#include <G4HEKaonZeroLongInelastic.hh>

Inheritance diagram for G4HEKaonZeroLongInelastic:

Public Member Functions
	G4HEKaonZeroLongInelastic ()

	~G4HEKaonZeroLongInelastic ()

virtual void	ModelDescription (std::ostream &) const

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

G4int	GetNumberOfSecondaries ()

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

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

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 52 of file G4HEKaonZeroLongInelastic.hh.

Constructor & Destructor Documentation

◆ G4HEKaonZeroLongInelastic()

G4HEKaonZeroLongInelastic::G4HEKaonZeroLongInelastic ( )

inline

Definition at line 55 of file G4HEKaonZeroLongInelastic.hh.

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

◆ ~G4HEKaonZeroLongInelastic()

G4HEKaonZeroLongInelastic::~G4HEKaonZeroLongInelastic ( )

inline

Definition at line 65 of file G4HEKaonZeroLongInelastic.hh.

65{};

Member Function Documentation

◆ ApplyYourself()

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

virtual

Implements G4HadronicInteraction.

Definition at line 58 of file G4HEKaonZeroLongInelastic.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 << "GHEKaonZeroLongInelastic: incident energy < 1 GeV " << G4endl;
 
  if(verboseLevel > 1) {
    G4cout << "G4HEKaonZeroLongInelastic::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; 
    
  // Split K0L into K0 and K0bar
  if (G4UniformRand() < 0.5)
    FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength,
                              incidentParticle, targetParticle);
  else
    FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength,
                          incidentParticle, targetParticle, atomicWeight);
 
  // Do nuclear interaction with either K0 or K0bar 
  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;
 
  // Check for K0, K0bar and change particle types to K0L, K0S if necessary
  G4int kcode;
  for (G4int i = 0; i < vecLength; i++) {
    kcode = pv[i].getCode();
    if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) {
      if (G4UniformRand() < 0.5) 
        pv[i] = KaonZeroShort; 
      else
        pv[i] = KaonZeroLong;
    }
  } 
 
  //      ................
 
  FillParticleChange(pv,  vecLength);
  delete [] pv;
  theParticleChange.SetStatusChange(stopAndKill);
  return &theParticleChange;
}

◆ FirstIntInCasAntiKaonZero()

void G4HEKaonZeroLongInelastic::FirstIntInCasAntiKaonZero	(	G4bool &	inElastic,
		const G4double	availableEnergy,
		G4HEVector	pv[],
		G4int &	vecLen,
		const G4HEVector &	incidentParticle,
		const G4HEVector &	targetParticle
	)

Definition at line 528 of file G4HEKaonZeroLongInelastic.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   numSec = 60;
 
  G4int neutronCode = Neutron.getCode();
  G4int protonCode = Proton.getCode();
  G4int kaonMinusCode = KaonMinus.getCode();
  G4int kaonZeroCode = KaonZero.getCode();
  G4int antiKaonZeroCode = AntiKaonZero.getCode();  
 
  G4int targetCode = targetParticle.getCode();
  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
 
  static G4bool first = true;
  static G4double protmul[numMul], protnorm[numSec];  // proton constants
  static G4double neutmul[numMul], neutnorm[numSec];  // neutron constants
 
  //  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];
    }
  }                                // end of initialization
 
  // initialize the first two particles
  // the same as beam and target                                    
  pv[0] = incidentParticle;
  pv[1] = targetParticle;
  vecLen = 2;
 
  if (!inElastic || (availableEnergy <= PionPlus.getMass())) 
    return;
                                        
  // Inelastic scattering
 
  npos = 0, nneg = 0, nzero = 0;
  G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15};
  G4int iplab = G4int( incidentTotalMomentum*5.);
  if ((iplab < 10) && (G4UniformRand() < cech[iplab]) ) {
    G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.));
    G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45,
                       0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33};
    if (G4UniformRand() < cnk0[ipl]) {
      if(targetCode == protonCode) {
        return;
      } else {
        pv[0] = KaonMinus;
        pv[1] = Proton;
        return;
      }
    }
 
    G4double ran = G4UniformRand();
    if(targetCode == protonCode) {
 
      // target is a proton 
      if( ran < 0.25 ) {
        ; 
      } else if (ran < 0.50) {
        pv[0] = PionPlus;
        pv[1] = SigmaZero;
      } else if (ran < 0.75) {
        ;
      } else {
        pv[0] = PionPlus;
        pv[1] = Lambda;
      }
    } else {
 
      // target is a neutron
      if( ran < 0.25 ) { 
        pv[0] = PionMinus;
        pv[1] = SigmaPlus;
      } else if (ran < 0.50) {
        pv[0] = PionZero;
        pv[1] = SigmaZero;
      } else if (ran < 0.75) { 
        pv[0] = PionPlus;
        pv[1] = SigmaMinus;
      } else {
        pv[0] = PionZero;
        pv[1] = Lambda;
      }
    }
    return;
 
  } 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>=1) && (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 == (1+nneg))
         {
           if( G4UniformRand() < 0.5)
             {
               pv[0] = KaonMinus;
             }
           else
             {
               pv[1] = Neutron;
             }
         }
       else      
         {
           pv[0] = KaonMinus;
           pv[1] = Neutron;
         } 
     }  
  else
     {
       if( npos == nneg)
         {
           if( G4UniformRand() < 0.75)
             {
             }
           else
             {
               pv[0] = KaonMinus;
               pv[1] = Proton;
             }
         } 
       else if ( npos == (1+nneg))
         {
           pv[0] = KaonMinus;
         }
       else
         {
           pv[1] = Proton;
         }
     }      
 
 
  if( G4UniformRand() < 0.5 )   
     {
       if(    (    (pv[0].getCode() == kaonMinusCode)
                && (pv[1].getCode() == neutronCode)  )
           || (    (pv[0].getCode() == kaonZeroCode)
                && (pv[1].getCode() == protonCode)   )
           || (    (pv[0].getCode() == antiKaonZeroCode)
                && (pv[1].getCode() == protonCode)   )   )
         {
           G4double ran = G4UniformRand();
           if( pv[1].getCode() == protonCode)
             { 
               if(ran < 0.68)
                 {
                   pv[0] = PionPlus;
                   pv[1] = Lambda;
                 }
               else if (ran < 0.84)
                 {
                   pv[0] = PionZero;
                   pv[1] = SigmaPlus;
                 }
               else
                 {
                   pv[0] = PionPlus;
                   pv[1] = SigmaZero;
                 }
             }
           else
             {
               if(ran < 0.68)
                 {
                   pv[0] = PionMinus;
                   pv[1] = Lambda;
                 }
               else if (ran < 0.84)
                 {
                   pv[0] = PionMinus;
                   pv[1] = SigmaZero;
                 }
               else
                 {
                   pv[0] = PionZero;
                   pv[1] = SigmaMinus;
                 }
             }
         } 
       else
         {
           G4double ran = G4UniformRand();
           if (ran < 0.67)
              {
                pv[0] = PionZero;
                pv[1] = Lambda;
              }
           else if (ran < 0.78)
              {
                pv[0] = PionMinus;
                pv[1] = SigmaPlus;
              }
           else if (ran < 0.89)
              {
                pv[0] = PionZero;
                pv[1] = SigmaZero;
              }
           else
              {
                pv[0] = PionPlus;
                pv[1] = SigmaMinus;
              }
         }
    }
 
  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().

◆ FirstIntInCasKaonZero()

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

Definition at line 206 of file G4HEKaonZeroLongInelastic.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   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 neutmul[numMul], neutnorm[numSec];  // neutron constants
 
  // 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-1); nneg<=(npos+1); 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=npos; nneg<=(npos+2); 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];
    }
  }    // end of initialization
 
 
  // Initialize the first two particles
  // the same as beam and target
  pv[0] = incidentParticle;
  pv[1] = targetParticle;
  vecLen = 2;
 
  if( !inElastic ) {
    // quasi-elastic scattering, no pions produced
    if( targetCode == protonCode) {
      G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07};
      G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) );
      if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) {
        // charge exchange  K+ n -> K0 p
        pv[0] = KaonPlus;
        pv[1] = Neutron;
      }
    }
    return;
  } else if (availableEnergy <= PionPlus.getMass()) {
    return;
  }
 
  // Inelastic scattering
 
  npos = 0, nneg = 0, nzero = 0;
  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 == neutronCode) {
      // target is a neutron
      w0 = - sqr(1.+protb)/(2.*c*c);
      w0 = std::exp(w0);
      wm = - sqr(-1.+protb)/(2.*c*c);
      wm = std::exp(wm);
      w0 = w0/2.;
      wm = wm*1.5;
      if (G4UniformRand() < w0/(w0+wm) ) {
        npos = 0;
        nneg = 0;
        nzero = 1;
      } else {
        npos = 0;
        nneg = 1;
        nzero = 0;
      }
 
    } else {
      // target is a proton
      w0 = -sqr(1.+neutb)/(2.*c*c);
      wp = w0 = 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-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*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=npos; nneg<=(npos+2); nneg++ )
                   {
                     for (nzero=0; nzero<numSec/3; nzero++) {
                       if (++counter < numMul) {
                         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*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 == neutronCode)
     {
       if( npos == nneg)
         {
         }
       else if (npos == (nneg-1))
         {
           if( G4UniformRand() < 0.5)
             {
               pv[0] = KaonPlus;
             }
           else
             {
               pv[1] = Proton;
             }
         }
       else
         {
           pv[0] = KaonPlus;
           pv[1] = Proton;
         }
     }
   else
     {
       if( npos == nneg )
         {
           if( G4UniformRand() < 0.25)
             {
               pv[0] = KaonPlus;
               pv[1] = Neutron;
             }
           else
             {
             }
         }
       else if ( npos == (nneg+1))
         {
           pv[1] = Neutron;
         }
       else
         {
           pv[0] = KaonPlus;
         }
     }
 
  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 G4HEKaonZeroLongInelastic::GetNumberOfSecondaries ( )

inline

Definition at line 74 of file G4HEKaonZeroLongInelastic.hh.

75 { return vecLength;};

◆ ModelDescription()

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

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 43 of file G4HEKaonZeroLongInelastic.cc.

{
  outFile << "G4HEKaonZeroLongInelastic is one of the High Energy\n"
          << "Parameterized (HEP) models used to implement inelastic\n"
          << "K0L 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 K0L with initial energies\n"
          << "above 20 GeV.\n";
}

Member Data Documentation

◆ vecLength

G4int G4HEKaonZeroLongInelastic::vecLength

Definition at line 69 of file G4HEKaonZeroLongInelastic.hh.

Referenced by ApplyYourself(), and GetNumberOfSecondaries().

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

Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4HEKaonZeroLongInelastic()

◆ ~G4HEKaonZeroLongInelastic()

Member Function Documentation

◆ ApplyYourself()

◆ FirstIntInCasAntiKaonZero()

◆ FirstIntInCasKaonZero()

◆ GetNumberOfSecondaries()

◆ ModelDescription()

Member Data Documentation

◆ vecLength