#include <G4NucleiModel.hh>

Public Types
typedef std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > >	modelLists

Public Member Functions
G4double	absorptionCrossSection (G4double e, G4int type) const

G4bool	empty () const

G4bool	forceFirst (const G4CascadParticle &cparticle) const

	G4NucleiModel ()

	G4NucleiModel (G4int a, G4int z)

	G4NucleiModel (G4InuclNuclei *nuclei)

void	generateModel (G4int a, G4int z)

void	generateModel (G4InuclNuclei *nuclei)

G4InuclElementaryParticle	generateNucleon (G4int type, G4int zone) const

G4LorentzVector	generateNucleonMomentum (G4int type, G4int zone) const

void	generateParticleFate (G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)

G4double	getDensity (G4int ip, G4int izone) const

G4double	getFermiKinetic (G4int ip, G4int izone) const

G4double	getFermiMomentum (G4int ip, G4int izone) const

G4int	getNumberOfNeutrons () const

G4int	getNumberOfProtons () const

G4int	getNumberOfZones () const

G4double	getPotential (G4int ip, G4int izone) const

G4double	getRadius () const

G4double	getRadius (G4int izone) const

G4double	getRadiusUnits () const

std::pair< G4int, G4int >	getTypesOfNucleonsInvolved () const

G4double	getVolume (G4int izone) const

G4int	getZone (G4double r) const

G4CascadParticle	initializeCascad (G4InuclElementaryParticle *particle)

void	initializeCascad (G4InuclNuclei bullet, G4InuclNuclei target, modelLists &output)

G4bool	isProjectile (const G4CascadParticle &cparticle) const

void	printModel () const

void	reset (G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)

void	setVerboseLevel (G4int verbose)

G4bool	stillInside (const G4CascadParticle &cparticle)

G4double	totalCrossSection (G4double ke, G4int rtype) const

G4bool	worthToPropagate (const G4CascadParticle &cparticle) const

virtual	~G4NucleiModel ()

Static Public Member Functions
static G4bool	useQuasiDeuteron (G4int ptype, G4int qdtype=0)

Protected Types
typedef std::pair< G4InuclElementaryParticle, G4double >	partner

Protected Member Functions
void	boundaryTransition (G4CascadParticle &cparticle)

void	choosePointAlongTraj (G4CascadParticle &cparticle)

void	fillBindingEnergies ()

void	fillPotentials (G4int type, G4double tot_vol)

void	fillZoneRadii (G4double nuclearRadius)

G4double	fillZoneVolumes (G4double nuclearRadius)

G4double	generateInteractionLength (const G4CascadParticle &cparticle, G4double path, G4double invmfp) const

void	generateInteractionPartners (G4CascadParticle &cparticle)

G4InuclElementaryParticle	generateQuasiDeuteron (G4int type1, G4int type2, G4int zone) const

G4double	getCurrentDensity (G4int ip, G4int izone) const

G4double	getRatio (G4int ip) const

G4double	inverseMeanFreePath (const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)

G4bool	passFermi (const std::vector< G4InuclElementaryParticle > &particles, G4int zone)

G4bool	passTrailing (const G4ThreeVector &hit_position)

G4double	zoneIntegralGaussian (G4double ur1, G4double ur2, G4double nuclearRadius) const

G4double	zoneIntegralWoodsSaxon (G4double ur1, G4double ur2, G4double nuclearRadius) const

Static Protected Member Functions
static G4bool	sortPartners (const partner &p1, const partner &p2)

Protected Attributes
std::vector< partner >	thePartners

Private Types
enum	PotentialType { WoodsSaxon =0 , Gaussian =1 }

Private Attributes
G4int	A

std::vector< G4double >	acsecs

std::vector< G4double >	binding_energies

std::vector< G4ThreeVector >	collisionPts

std::vector< G4ThreeVector >	coordinates

const G4double	crossSectionUnits

G4int	current_nucl1

G4int	current_nucl2

G4LorentzConvertor	dummy_convertor

G4CollisionOutput	EPCoutput

std::vector< std::vector< G4double > >	fermi_momenta

const G4double	fermiMomentum

G4CascadeInterpolator< 30 >	gammaQDinterp

const G4double	gammaQDscale

std::vector< G4LorentzVector >	momentums

const G4InuclElementaryParticle	neutronEP

G4int	neutronNumber

G4int	neutronNumberCurrent

G4double	nuclei_radius

G4double	nuclei_volume

std::vector< std::vector< G4double > >	nucleon_densities

G4int	number_of_zones

std::vector< G4double >	pf

const G4double	potentialThickness

const G4InuclElementaryParticle	protonEP

G4int	protonNumber

G4int	protonNumberCurrent

std::vector< G4InuclElementaryParticle >	qdeutrons

const G4double	R_nucleon

const G4double	radiusForSmall

const G4double	radiusScale

const G4double	radiusScale2

const G4double	radiusUnits

const G4double	radScaleAlpha

std::vector< G4InuclElementaryParticle >	raw_particles

std::vector< G4double >	rod

const G4double	skinDepth

G4InuclNuclei *	theNucleus

G4double	ur [7]

G4double	v [6]

G4double	v1 [6]

G4int	verboseLevel

std::vector< G4double >	vz

G4int	Z

std::vector< std::vector< G4double > >	zone_potentials

std::vector< G4double >	zone_radii

std::vector< G4double >	zone_volumes

Static Private Attributes
static const G4double	alfa3 [3] = { 0.7, 0.3, 0.01 }

static const G4double	alfa6 [6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 }

static const G4double	hyperon_vp = 0.030

static const G4double	kaon_vp = 0.015

static const G4double	large = 1000.

static const G4double	pion_vp = 0.007

static const G4double	pion_vp_small = 0.007

static const G4double	piTimes4thirds = pi*4./3.

static const G4double	small = 1.0e-9

Detailed Description

Definition at line 91 of file G4NucleiModel.hh.

Member Typedef Documentation

◆ modelLists

typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > G4NucleiModel::modelLists

Definition at line 161 of file G4NucleiModel.hh.

◆ partner

typedef std::pair<G4InuclElementaryParticle, G4double> G4NucleiModel::partner

protected

Definition at line 203 of file G4NucleiModel.hh.

Member Enumeration Documentation

◆ PotentialType

enum G4NucleiModel::PotentialType

private

Enumerator
WoodsSaxon
Gaussian

Definition at line 295 of file G4NucleiModel.hh.

295{ WoodsSaxon=0, Gaussian=1 };

G4NucleiModel::WoodsSaxon

@ WoodsSaxon

Definition: G4NucleiModel.hh:295

G4NucleiModel::Gaussian

@ Gaussian

Definition: G4NucleiModel.hh:295

Constructor & Destructor Documentation

◆ G4NucleiModel() [1/3]

G4NucleiModel::G4NucleiModel ( )

Definition at line 232 of file G4NucleiModel.cc.

  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
    A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
    current_nucl2(0), gammaQDinterp(kebins),
    crossSectionUnits(G4CascadeParameters::xsecScale()),
    radiusUnits(G4CascadeParameters::radiusScale()),
    skinDepth(0.611207*radiusUnits),
    radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
    radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
    radiusForSmall(G4CascadeParameters::radiusSmall()),
    radScaleAlpha(G4CascadeParameters::radiusAlpha()),
    fermiMomentum(G4CascadeParameters::fermiScale()),
    R_nucleon(G4CascadeParameters::radiusTrailing()),
    gammaQDscale(G4CascadeParameters::gammaQDScale()),
    potentialThickness(1.0),
    neutronEP(neutron), protonEP(proton) {}

◆ G4NucleiModel() [2/3]

G4NucleiModel::G4NucleiModel	(	G4int	a,
		G4int	z
	)

Definition at line 250 of file G4NucleiModel.cc.

  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
    A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
    current_nucl2(0), gammaQDinterp(kebins),
    crossSectionUnits(G4CascadeParameters::xsecScale()),
    radiusUnits(G4CascadeParameters::radiusScale()),
    skinDepth(0.611207*radiusUnits),
    radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
    radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
    radiusForSmall(G4CascadeParameters::radiusSmall()),
    radScaleAlpha(G4CascadeParameters::radiusAlpha()),
    fermiMomentum(G4CascadeParameters::fermiScale()),
    R_nucleon(G4CascadeParameters::radiusTrailing()),
    gammaQDscale(G4CascadeParameters::gammaQDScale()),
    potentialThickness(1.0),
    neutronEP(neutron), protonEP(proton) {
  generateModel(a,z);
}

References generateModel().

◆ G4NucleiModel() [3/3]

G4NucleiModel::G4NucleiModel ( G4InuclNuclei * nuclei )

explicit

Definition at line 270 of file G4NucleiModel.cc.

  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
    A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
    current_nucl2(0), gammaQDinterp(kebins),
    crossSectionUnits(G4CascadeParameters::xsecScale()),
    radiusUnits(G4CascadeParameters::radiusScale()),
    skinDepth(0.611207*radiusUnits),
    radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
    radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
    radiusForSmall(G4CascadeParameters::radiusSmall()),
    radScaleAlpha(G4CascadeParameters::radiusAlpha()),
    fermiMomentum(G4CascadeParameters::fermiScale()),
    R_nucleon(G4CascadeParameters::radiusTrailing()),
    gammaQDscale(G4CascadeParameters::gammaQDScale()),
    potentialThickness(1.0),
    neutronEP(neutron), protonEP(proton) {
  generateModel(nuclei);
}

References generateModel(), and G4InuclParticleNames::nuclei.

◆ ~G4NucleiModel()

G4NucleiModel::~G4NucleiModel ( )

virtual

Definition at line 290 of file G4NucleiModel.cc.

                              {
  delete theNucleus;
  theNucleus = 0;
}

References theNucleus.

Member Function Documentation

◆ absorptionCrossSection()

G4double G4NucleiModel::absorptionCrossSection	(	G4double	e,
		G4int	type
	)		const

Definition at line 1922 of file G4NucleiModel.cc.

                                                                            {
  if (!useQuasiDeuteron(type)) {
    G4cerr << "absorptionCrossSection() only valid for incident pions or gammas" 
           << G4endl;
    return 0.;
  }
 
  G4double csec = 0.;
 
  // Pion absorption is parametrized for low vs. medium energy
  // ... use for muon capture as well
  if (type == pionPlus || type == pionMinus || type == pionZero ||
      type == muonMinus) {
    if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
              + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
    else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);     
  }
 
  if (type == photon) {
    csec = gammaQDinterp.interpolate(ke, gammaQDxsec) * gammaQDscale;
  }
 
  if (csec < 0.0) csec = 0.0;
 
  if (verboseLevel > 2) {
    G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;   
  }
 
  return crossSectionUnits * csec;
}

References crossSectionUnits, G4cerr, G4cout, G4endl, gammaQDinterp, gammaQDscale, anonymous_namespace{G4NucleiModel.cc}::gammaQDxsec, G4CascadeInterpolator< NBINS >::interpolate(), G4InuclParticleNames::muonMinus, G4InuclParticleNames::photon, G4InuclParticleNames::pionMinus, G4InuclParticleNames::pionPlus, G4InuclParticleNames::pionZero, useQuasiDeuteron(), and verboseLevel.

Referenced by inverseMeanFreePath().

◆ boundaryTransition()

void G4NucleiModel::boundaryTransition ( G4CascadParticle & cparticle )

protected

Definition at line 1115 of file G4NucleiModel.cc.

                                                                  {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
  }
 
  G4int zone = cparticle.getCurrentZone();
 
  if (cparticle.movingInsideNuclei() && zone == 0) {
    if (verboseLevel) G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
    return;
  }
 
  G4LorentzVector mom = cparticle.getMomentum();
  G4ThreeVector pos = cparticle.getPosition();
  
  G4int type = cparticle.getParticle().type();
  
  G4double r = pos.mag();
  G4double p = mom.vect().mag();           // NAT
  G4double pr = pos.dot(mom.vect()) / r;
  G4double pperp2 = p*p - pr*pr;           // NAT
 
  G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
 
  // NOTE:  dv is the height of the wall seen by the particle  
  G4double dv = getPotential(type,next_zone) - getPotential(type,zone);
  if (verboseLevel > 3) {
    G4cout << "Potentials for type " << type << " = " 
       << getPotential(type,zone) << " , "
       << getPotential(type,next_zone) << G4endl;
  }
 
  G4double qv = dv*dv + 2.0*dv*mom.e() + pr*pr;
  //  G4double qv = dv*dv + 2.0*dv*mom.m() + pr*pr;    // more correct? NAT
  G4double p1r = 0.; 
 
  // Perpendicular contribution to pr^2 after penetrating     // NAT
  // potential, to leading order in thickness                 // NAT
  G4double qperp = 2.0*pperp2*potentialThickness/r;           // NAT 
 
  if (verboseLevel > 3) {
    G4cout << " type " << type << " zone " << zone << " next " << next_zone
       << " qv " << qv << " dv " << dv << G4endl;
  }
 
  G4bool adjustpperp = false;                                   // NAT
  G4double smallish = 0.001;                                    // NAT
 
//  if (qv <= 0.0) {    // reflection
  if (qv <= 0.0 && qv+qperp <=0 ) {     // reflection         // NAT
    if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
    p1r = -pr;
    cparticle.incrementReflectionCounter();
//  } else {            // transition
 
  } else if (qv > 0.0) {        // standard transition  // NAT
    if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
    p1r = std::sqrt(qv);
    if (pr < 0.0) p1r = -p1r;
    cparticle.updateZone(next_zone);
    cparticle.resetReflection();
 
  } else {   // transition via transverse kinetic energy (allowed for thick walls)  // NAT
    if (verboseLevel > 3) G4cout << " passes thru boundary due to angular momentum" << G4endl;
    p1r = smallish * pr; // don't want exactly tangent momentum
    adjustpperp = true;
 
    cparticle.updateZone(next_zone);
    cparticle.resetReflection();
  }
  
  G4double prr = (p1r - pr)/r;  // change to radial momentum, divided by r
  
  if (verboseLevel > 3) {
    G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
       << prr*pos.y()  << " pz " << prr*pos.z() << " mag "
       << std::fabs(prr*r) << G4endl;
  }
 
  if (adjustpperp) {                       // NAT
    G4ThreeVector old_pperp = mom.vect() - pos*(pr/r);
    G4double new_pperp_mag = std::sqrt(std::max(0.0, pperp2 + qv - p1r*p1r) );
    // new total momentum found by rescaling p_perp
    mom.setVect(old_pperp * new_pperp_mag/std::sqrt(pperp2));
    // add a small radial component to make sure that we propagate into new zone
    mom.setVect(mom.vect() + pos*p1r/r);
  } else {
    mom.setVect(mom.vect() + pos*prr);
  }
 
  cparticle.updateParticleMomentum(mom);
}

References CLHEP::HepLorentzVector::e(), G4cerr, G4cout, G4endl, G4CascadParticle::getCurrentZone(), G4CascadParticle::getMomentum(), G4CascadParticle::getParticle(), G4CascadParticle::getPosition(), getPotential(), G4CascadParticle::incrementReflectionCounter(), CLHEP::Hep3Vector::mag(), G4INCL::Math::max(), G4CascadParticle::movingInsideNuclei(), pos, potentialThickness, G4CascadParticle::resetReflection(), CLHEP::HepLorentzVector::setVect(), G4InuclElementaryParticle::type(), G4CascadParticle::updateParticleMomentum(), G4CascadParticle::updateZone(), CLHEP::HepLorentzVector::vect(), and verboseLevel.

Referenced by generateParticleFate().

◆ choosePointAlongTraj()

void G4NucleiModel::choosePointAlongTraj ( G4CascadParticle & cparticle )

protected

Definition at line 1212 of file G4NucleiModel.cc.

                                                                    {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::choosePointAlongTraj" << G4endl;
 
  // Get trajectory through nucleus by computing exit point of line,
  // assuming that current position is on surface
 
  // FIXME:  These need to be reusable (static) buffers
  G4ThreeVector pos  = cparticle.getPosition();
  G4ThreeVector rhat = pos.unit();
 
  G4ThreeVector phat = cparticle.getMomentum().vect().unit();
  if (cparticle.getMomentum().vect().mag() < small) phat.set(0.,0.,1.);
 
  if (verboseLevel > 3)
    G4cout << " pos " << pos << " phat " << phat << " rhat " << rhat << G4endl;
 
  G4ThreeVector posout = pos;
  G4double prang = rhat.angle(-phat);
 
  if (prang < 1e-6) posout = -pos;      // Check for radial incidence
  else {
    G4double posrot = 2.*prang - pi;
    posout.rotate(posrot, phat.cross(rhat));
    if (verboseLevel > 3) G4cout << " posrot " << posrot/deg << " deg";
  }
 
  if (verboseLevel > 3) G4cout << " posout " << posout << G4endl;
 
  // Get list of zone crossings along trajectory
  G4ThreeVector posmid = (pos+posout)/2.;   // Midpoint of trajectory
  G4double r2mid = posmid.mag2();
  G4double lenmid = (posout-pos).mag()/2.;  // Half-length of trajectory
 
  G4int zoneout = number_of_zones-1;
  G4int zonemid = getZone(posmid.mag());    // Middle zone traversed
 
  // Every zone is entered then exited, so preallocate vector
  G4int ncross = (number_of_zones-zonemid)*2;
 
  if (verboseLevel > 3) {
    G4cout << " posmid " << posmid << " lenmid " << lenmid
       << " zoneout " << zoneout << " zonemid " << zonemid
       << " ncross " << ncross << G4endl;
  }
 
  // FIXME:  These need to be reusable (static) buffers
  std::vector<G4double> wtlen(ncross,0.);   // CDF from entry point
  std::vector<G4double> len(ncross,0.);     // Distance from entry point
 
  // Work from outside in, to accumulate CDF steps properly
  G4int i;              // Loop variable, used multiple times
  for (i=0; i<ncross/2; i++) {
    G4int iz = zoneout-i;
    G4double ds = std::sqrt(zone_radii[iz]*zone_radii[iz]-r2mid);
 
    len[i] = lenmid - ds;       // Distance from entry to crossing
    len[ncross-1-i] = lenmid + ds;  // Distance to outbound crossing
 
    if (verboseLevel > 3) {
      G4cout << " i " << i << " iz " << iz << " ds " << ds
         << " len " << len[i] << G4endl;
    }
  }
 
  // Compute weights for each zone along trajectory and integrate
  for (i=1; i<ncross; i++) {
    G4int iz = (i<ncross/2) ? zoneout-i+1 : zoneout-ncross+i+1;
 
    G4double dlen = len[i]-len[i-1];    // Distance from previous crossing
 
    // Uniform probability across entire zone
    //*** G4double wt = dlen*(getDensity(1,iz)+getDensity(2,iz));
 
    // Probability based on interaction length through zone
    G4double invmfp = (inverseMeanFreePath(cparticle, neutronEP, iz)
               + inverseMeanFreePath(cparticle, protonEP, iz));
 
    // Integral of exp(-len/mfp) from start of zone to end
    G4double wt = (G4Exp(-len[i-1]*invmfp)-G4Exp(-len[i]*invmfp)) / invmfp;
 
    wtlen[i] = wtlen[i-1] + wt;
 
    if (verboseLevel > 3) {
      G4cout << " i " << i << " iz " << iz << " avg.mfp " << 1./invmfp
         << " dlen " << dlen  << " wt " << wt << " wtlen " << wtlen[i]
         << G4endl;
    }
  }
 
  // Normalize CDF to unit integral
  std::transform(wtlen.begin(), wtlen.end(), wtlen.begin(),
         std::bind(std::divides<G4double>(), std::placeholders::_1, wtlen.back()));
  
  if (verboseLevel > 3) {
    G4cout << " weights";
    for (i=0; i<ncross; i++) G4cout << " " << wtlen[i];
    G4cout << G4endl;
  }
  
  // Choose random point along trajectory, weighted by density
  G4double rand = G4UniformRand();
  G4int ir = std::upper_bound(wtlen.begin(),wtlen.end(),rand) - wtlen.begin();
 
  G4double frac = (rand-wtlen[ir-1]) / (wtlen[ir]-wtlen[ir-1]);
  G4double drand = (1.-frac)*len[ir-1] + frac*len[ir];
 
  if (verboseLevel > 3) {
    G4cout << " rand " << rand << " ir " << ir << " frac " << frac
       << " drand " << drand << G4endl;
  }
 
  pos += drand * phat;      // Distance from entry point along trajectory
 
  // Update input particle with new location
  cparticle.updatePosition(pos);
  cparticle.updateZone(getZone(pos.mag()));
 
  if (verboseLevel > 2) {
    G4cout << " moved particle to zone " << cparticle.getCurrentZone() 
       << " @ " << pos << G4endl;
  }
}

References CLHEP::Hep3Vector::angle(), CLHEP::Hep3Vector::cross(), deg, anonymous_namespace{G4QuasiElRatios.cc}::ds, G4cout, G4endl, G4Exp(), G4UniformRand, G4CascadParticle::getCurrentZone(), G4CascadParticle::getMomentum(), G4CascadParticle::getPosition(), getZone(), inverseMeanFreePath(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), neutronEP, number_of_zones, pi, pos, protonEP, CLHEP::Hep3Vector::rotate(), CLHEP::Hep3Vector::set(), small, G4coutFormatters::anonymous_namespace{G4coutFormatters.cc}::transform(), CLHEP::Hep3Vector::unit(), G4CascadParticle::updatePosition(), G4CascadParticle::updateZone(), CLHEP::HepLorentzVector::vect(), verboseLevel, and zone_radii.

Referenced by initializeCascad().

◆ empty()

G4bool G4NucleiModel::empty ( ) const

inline

Definition at line 150 of file G4NucleiModel.hh.

                       { 
    return neutronNumberCurrent < 1 && protonNumberCurrent < 1; 
  }

References neutronNumberCurrent, and protonNumberCurrent.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ fillBindingEnergies()

void G4NucleiModel::fillBindingEnergies ( )

protected

Definition at line 390 of file G4NucleiModel.cc.

                                        {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
 
  G4double dm = bindingEnergy(A,Z);
 
  // Binding energy differences for proton and neutron loss, respectively
  // FIXME:  Why is fabs() used here instead of the signed difference?
  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
}

References A, binding_energies, G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4endl, GeV, verboseLevel, and Z.

Referenced by generateModel().

◆ fillPotentials()

void G4NucleiModel::fillPotentials	(	G4int	type,
		G4double	tot_vol
	)

protected

Definition at line 479 of file G4NucleiModel.cc.

                                                               {
  if (verboseLevel > 1) 
    G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
 
  if (type != proton && type != neutron) return;
 
  const G4double mass = G4InuclElementaryParticle::getParticleMass(type);
 
  // FIXME:  This is the fabs() binding energy difference, not signed
  const G4double dm = binding_energies[type-1];
 
  rod.clear(); rod.reserve(number_of_zones);
  pf.clear();  pf.reserve(number_of_zones);
  vz.clear();  vz.reserve(number_of_zones);
 
  G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
  G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
 
  for (G4int i = 0; i < number_of_zones; i++) {
    G4double rd = dd0 * v[i] / v1[i];
    rod.push_back(rd);
    G4double pff = fermiMomentum * G4cbrt(rd);
    pf.push_back(pff);
    vz.push_back(0.5 * pff * pff / mass + dm);
  }
  
  nucleon_densities.push_back(rod);
  fermi_momenta.push_back(pf);
  zone_potentials.push_back(vz);
}

References binding_energies, fermi_momenta, fermiMomentum, G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, G4InuclElementaryParticle::getParticleMass(), G4InuclParticleNames::neutron, neutronNumber, nucleon_densities, number_of_zones, pf, piTimes4thirds, G4InuclParticleNames::proton, protonNumber, rod, v, v1, verboseLevel, vz, and zone_potentials.

Referenced by generateModel().

◆ fillZoneRadii()

void G4NucleiModel::fillZoneRadii ( G4double nuclearRadius )

protected

Definition at line 404 of file G4NucleiModel.cc.

                                                        {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
 
  G4double skinRatio = nuclearRadius/skinDepth;
  G4double skinDecay = G4Exp(-skinRatio);    
 
  if (A < 5) {          // Light ions treated as simple balls
    zone_radii.push_back(nuclearRadius);
    ur[0] = 0.;
    ur[1] = 1.;
  } else if (A < 12) {      // Small nuclei have Gaussian potential
    G4double rSq = nuclearRadius * nuclearRadius;
    G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
    
    ur[0] = 0.0;
    for (G4int i = 0; i < number_of_zones; i++) {
      G4double y = std::sqrt(-G4Log(alfa3[i]));
      zone_radii.push_back(gaussRadius * y);
      ur[i+1] = y;
    }
  } else if (A < 100) {     // Intermediate nuclei have three-zone W-S
    ur[0] = -skinRatio;
    for (G4int i = 0; i < number_of_zones; i++) {
      G4double y = G4Log((1.0 + skinDecay)/alfa3[i] - 1.0);
      zone_radii.push_back(nuclearRadius + skinDepth * y);
      ur[i+1] = y;
    }
  } else {          // Heavy nuclei have six-zone W-S
    ur[0] = -skinRatio;
    for (G4int i = 0; i < number_of_zones; i++) {
      G4double y = G4Log((1.0 + skinDecay)/alfa6[i] - 1.0);
      zone_radii.push_back(nuclearRadius + skinDepth * y);
      ur[i+1] = y;
    }
  }
}

References A, alfa3, alfa6, G4cout, G4endl, G4Exp(), G4Log(), number_of_zones, skinDepth, ur, verboseLevel, and zone_radii.

Referenced by generateModel().

◆ fillZoneVolumes()

G4double G4NucleiModel::fillZoneVolumes ( G4double nuclearRadius )

protected

Definition at line 444 of file G4NucleiModel.cc.

                                                              {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
 
  G4double tot_vol = 0.;    // Return value omitting 4pi/3 for sphere!
 
  if (A < 5) {          // Light ions are treated as simple balls
    v[0] = v1[0] = 1.;
    tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
    zone_volumes.push_back(tot_vol*piTimes4thirds); // True volume of sphere
 
    return tot_vol;
  }
 
  PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
 
  for (G4int i = 0; i < number_of_zones; i++) {
    if (usePotential == WoodsSaxon) {
      v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
    } else {
      v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
    }
    
    tot_vol += v[i];
    
    v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
    if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
 
    zone_volumes.push_back(v1[i]*piTimes4thirds);   // True volume of shell
  }
 
  return tot_vol;   // Sum of zone integrals, not geometric volume
}

References A, G4cout, G4endl, Gaussian, number_of_zones, piTimes4thirds, ur, v, v1, verboseLevel, WoodsSaxon, zone_radii, zone_volumes, zoneIntegralGaussian(), and zoneIntegralWoodsSaxon().

Referenced by generateModel().

◆ forceFirst()

G4bool G4NucleiModel::forceFirst ( const G4CascadParticle & cparticle ) const

Definition at line 1338 of file G4NucleiModel.cc.

                                                                        {
  return (isProjectile(cparticle) && 
      (cparticle.getParticle().isPhoton() ||
       cparticle.getParticle().isMuon())
      );
}

References G4CascadParticle::getParticle(), G4InuclElementaryParticle::isMuon(), G4InuclElementaryParticle::isPhoton(), and isProjectile().

Referenced by generateInteractionLength(), and initializeCascad().

◆ generateInteractionLength()

G4double G4NucleiModel::generateInteractionLength	(	const G4CascadParticle &	cparticle,
		G4double	path,
		G4double	invmfp
	)		const

protected

Definition at line 1890 of file G4NucleiModel.cc.

                                                           {
  // Delay interactions of newly formed secondaries (minimum int. length)
  const G4double young_cut = std::sqrt(10.0) * 0.25;
  const G4double huge_num = 50.0;   // Argument to exponential
 
  G4double spath = large;       // Buffer for return value
 
  if (invmfp < small) return spath; // No interaction, avoid unnecessary work
 
  G4double pw = -path * invmfp;     // Ratio of path in zone to MFP
  if (pw < -huge_num) pw = -huge_num;
  pw = 1.0 - G4Exp(pw);
    
  if (verboseLevel > 2) 
    G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
    
  // Primary particle(s) should always interact at least once
  if (forceFirst(cparticle) || (inuclRndm() < pw)) {
    spath = -G4Log(1.0 - pw * inuclRndm()) / invmfp;
    if (cparticle.young(young_cut, spath)) spath = large;
    
    if (verboseLevel > 2) 
      G4cout << " spath " << spath << " path " << path << G4endl;
  }
 
  return spath;
}

References forceFirst(), G4cout, G4endl, G4Exp(), G4Log(), G4InuclSpecialFunctions::inuclRndm(), large, small, verboseLevel, and G4CascadParticle::young().

Referenced by generateInteractionPartners().

◆ generateInteractionPartners()

void G4NucleiModel::generateInteractionPartners ( G4CascadParticle & cparticle )

protected

Definition at line 694 of file G4NucleiModel.cc.

                                                                      {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
  }
 
  thePartners.clear();      // Reset buffer for next cycle
 
  G4int ptype = cparticle.getParticle().type();
  G4int zone = cparticle.getCurrentZone();
 
  G4double r_in;        // Destination radius if inbound
  G4double r_out;       // Destination radius if outbound
 
  if (zone == number_of_zones) {
    r_in = nuclei_radius;
    r_out = 0.0;
  } else if (zone == 0) {           // particle is outside core
    r_in = 0.0;
    r_out = zone_radii[0];
  } else {
    r_in = zone_radii[zone - 1];
    r_out = zone_radii[zone];
  }  
 
  G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
 
  if (verboseLevel > 2) {
    if (isProjectile(cparticle)) G4cout << " incident particle: ";
    G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
       << G4endl;
  }
 
  if (path < -small) {          // something wrong
    if (verboseLevel)
      G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
    return;
  }
 
  if (std::fabs(path) < small) {    // Not moving, or just at boundary
    if (cparticle.getMomentum().vect().mag() > small) {
      if (verboseLevel > 3)
    G4cout << " generateInteractionPartners-> zero path" << G4endl;
      
      thePartners.push_back(partner()); // Dummy list terminator with zero path
      return;
    }
 
    if (zone >= number_of_zones)    // Place captured-at-rest in nucleus
      zone = number_of_zones-1;
  }
 
  G4double invmfp = 0.;         // Buffers for interaction probability
  G4double spath  = 0.;
  for (G4int ip = 1; ip < 3; ip++) {
    // Only process nucleons which remain active in target
    if (ip==proton  && protonNumberCurrent < 1) continue;
    if (ip==neutron && neutronNumberCurrent < 1) continue;
    if (ip==neutron && ptype==muonMinus) continue;  // mu-/n forbidden
 
    // All nucleons are assumed to be at rest when colliding
    G4InuclElementaryParticle particle = generateNucleon(ip, zone);
    invmfp = inverseMeanFreePath(cparticle, particle);
    spath  = generateInteractionLength(cparticle, path, invmfp);
 
    if (path<small || spath < path) {
      if (verboseLevel > 3) {
    G4cout << " adding partner[" << thePartners.size() << "]: "
           << particle << G4endl;
      }
      thePartners.push_back(partner(particle, spath));
    }
  } // for (G4int ip...
  
  if (verboseLevel > 2) {
    G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
  }
 
  // Absorption possible for pions or photons interacting with dibaryons
  if (useQuasiDeuteron(cparticle.getParticle().type())) {
    if (verboseLevel > 2) {
      G4cout << " trying quasi-deuterons with bullet: "
         << cparticle.getParticle() << G4endl;
    }
  
    // Initialize buffers for quasi-deuteron results
    qdeutrons.clear();
    acsecs.clear();
  
    G4double tot_invmfp = 0.0;      // Total inv. mean-free-path for all QDs
 
    // Proton-proton state interacts with pi-, mu- or neutrals
    if (protonNumberCurrent >= 2 && ptype != pip) {
      G4InuclElementaryParticle ppd = generateQuasiDeuteron(pro, pro, zone);
      if (verboseLevel > 2)
    G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
      
      invmfp = inverseMeanFreePath(cparticle, ppd);      
      tot_invmfp += invmfp;
      acsecs.push_back(invmfp);
      qdeutrons.push_back(ppd);
    }
    
    // Proton-neutron state interacts with any pion type or photon
    if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
      G4InuclElementaryParticle npd = generateQuasiDeuteron(pro, neu, zone);
      if (verboseLevel > 2) 
    G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
      
      invmfp = inverseMeanFreePath(cparticle, npd);
      tot_invmfp += invmfp;
      acsecs.push_back(invmfp);
      qdeutrons.push_back(npd);
    }
    
    // Neutron-neutron state interacts with pi+ or neutrals
    if (neutronNumberCurrent >= 2 && ptype != pim && ptype != mum) {
      G4InuclElementaryParticle nnd = generateQuasiDeuteron(neu, neu, zone);
      if (verboseLevel > 2)
    G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
      
      invmfp = inverseMeanFreePath(cparticle, nnd);
      tot_invmfp += invmfp;
      acsecs.push_back(invmfp);
      qdeutrons.push_back(nnd);
    } 
    
    // Select quasideuteron interaction from non-zero cross-section choices
    if (verboseLevel > 2) {
      for (size_t i=0; i<qdeutrons.size(); i++) {
    G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
           << "] " << acsecs[i];
      }
      G4cout << G4endl;
    }
  
    if (tot_invmfp > small) {       // Must have absorption cross-section
      G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
      
      if (path<small || apath < path) {     // choose the qdeutron
    G4double sl = inuclRndm() * tot_invmfp;
    G4double as = 0.0;
    
    for (size_t i = 0; i < qdeutrons.size(); i++) {
      as += acsecs[i];
      if (sl < as) { 
        if (verboseLevel > 2)
          G4cout << " deut type " << qdeutrons[i] << G4endl; 
 
        thePartners.push_back(partner(qdeutrons[i], apath));
        break;
      }
    }   // for (qdeutrons...
      }     // if (apath < path)
    }       // if (tot_invmfp > small)
  }     // if (useQuasiDeuteron) [pion, muon or photon]
  
  if (verboseLevel > 2) {
    G4cout << " after deuterons " << thePartners.size() << " partners"
       << G4endl;
  }
  
  if (thePartners.size() > 1) {     // Sort list by path length
    std::sort(thePartners.begin(), thePartners.end(), sortPartners);
  }
  
  G4InuclElementaryParticle particle;       // Total path at end of list
  thePartners.push_back(partner(particle, path));
}

Referenced by generateParticleFate().

◆ generateModel() [1/2]

void G4NucleiModel::generateModel	(	G4int	a,
		G4int	z
	)

Definition at line 315 of file G4NucleiModel.cc.

                                                  {
  if (verboseLevel) {
    G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
       << G4endl;
  }
 
  // If model already built, just return; otherwise initialize everything
  if (a == A && z == Z) {
    if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
    reset();
    return;
  }
 
  A = a;
  Z = z;
  delete theNucleus;
  theNucleus = new G4InuclNuclei(A,Z);      // For conservation checking
 
  neutronNumber = A - Z;
  protonNumber = Z;
  reset();
 
  if (verboseLevel > 3) {
    G4cout << "  crossSectionUnits = " << crossSectionUnits << G4endl
       << "  radiusUnits = " << radiusUnits << G4endl
       << "  skinDepth = " << skinDepth << G4endl
       << "  radiusScale = " << radiusScale << G4endl
       << "  radiusScale2 = " << radiusScale2 << G4endl
       << "  radiusForSmall = " << radiusForSmall << G4endl
       << "  radScaleAlpha  = " << radScaleAlpha << G4endl
       << "  fermiMomentum = " << fermiMomentum << G4endl
       << "  piTimes4thirds = " << piTimes4thirds << G4endl;
  }
 
  G4double nuclearRadius;       // Nuclear radius computed from A
  if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
  else     nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
 
  // This will be used to pre-allocate lots of arrays below
  number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
 
  // Clear all parameters arrays for reloading
  binding_energies.clear();
  nucleon_densities.clear();
  zone_potentials.clear();
  fermi_momenta.clear();
  zone_radii.clear();
  zone_volumes.clear();
 
  fillBindingEnergies();
  fillZoneRadii(nuclearRadius);
 
  G4double tot_vol = fillZoneVolumes(nuclearRadius);    // Woods-Saxon integral
 
  fillPotentials(proton, tot_vol);      // Protons
  fillPotentials(neutron, tot_vol);     // Neutrons
 
  // Additional flat zone potentials for other hadrons
  const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
  const std::vector<G4double> kp(number_of_zones, kaon_vp);
  const std::vector<G4double> hp(number_of_zones, hyperon_vp);
 
  zone_potentials.push_back(vp);
  zone_potentials.push_back(kp);
  zone_potentials.push_back(hp);
 
  nuclei_radius = zone_radii.back();
  nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
 
  if (verboseLevel > 3) printModel();
}

References A, binding_energies, crossSectionUnits, fermi_momenta, fermiMomentum, fillBindingEnergies(), fillPotentials(), fillZoneRadii(), fillZoneVolumes(), G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, hyperon_vp, kaon_vp, G4InuclParticleNames::neutron, neutronNumber, nuclei_radius, nuclei_volume, nucleon_densities, number_of_zones, pion_vp, pion_vp_small, piTimes4thirds, printModel(), G4InuclParticleNames::proton, protonNumber, radiusForSmall, radiusScale, radiusScale2, radiusUnits, radScaleAlpha, reset(), skinDepth, theNucleus, verboseLevel, Z, zone_potentials, zone_radii, and zone_volumes.

◆ generateModel() [2/2]

void G4NucleiModel::generateModel ( G4InuclNuclei * nuclei )

Definition at line 311 of file G4NucleiModel.cc.

                                                       {
  generateModel(nuclei->getA(), nuclei->getZ());
}

References generateModel(), and G4InuclParticleNames::nuclei.

Referenced by G4NucleiModel(), generateModel(), and G4IntraNucleiCascader::initialize().

◆ generateNucleon()

G4InuclElementaryParticle G4NucleiModel::generateNucleon	(	G4int	type,
		G4int	zone
	)		const

Definition at line 657 of file G4NucleiModel.cc.

                                                           {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
  }
 
  G4LorentzVector mom = generateNucleonMomentum(type, zone);
  return G4InuclElementaryParticle(mom, type);
}

References G4cout, G4endl, generateNucleonMomentum(), and verboseLevel.

Referenced by generateInteractionPartners().

◆ generateNucleonMomentum()

G4LorentzVector G4NucleiModel::generateNucleonMomentum	(	G4int	type,
		G4int	zone
	)		const

Definition at line 648 of file G4NucleiModel.cc.

                                                                   {
  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
  G4double mass = G4InuclElementaryParticle::getParticleMass(type);
 
  return generateWithRandomAngles(pmod, mass);
}

References G4InuclSpecialFunctions::G4cbrt(), G4InuclSpecialFunctions::generateWithRandomAngles(), getFermiMomentum(), G4InuclElementaryParticle::getParticleMass(), and G4InuclSpecialFunctions::inuclRndm().

Referenced by generateNucleon(), and generateQuasiDeuteron().

◆ generateParticleFate()

void G4NucleiModel::generateParticleFate	(	G4CascadParticle &	cparticle,
		G4ElementaryParticleCollider *	theEPCollider,
		std::vector< G4CascadParticle > &	cascade
	)

Definition at line 864 of file G4NucleiModel.cc.

                                                               {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
 
  if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
 
  // Create four-vector checking
#ifdef G4CASCADE_CHECK_ECONS
  G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel");  // Second arg is in GeV
  balance.setVerboseLevel(verboseLevel);
#endif
 
  outgoing_cparticles.clear();      // Clear return buffer for this event
  generateInteractionPartners(cparticle);   // Fills "thePartners" data
 
  if (thePartners.empty()) { // smth. is wrong -> needs special treatment
    if (verboseLevel)
      G4cerr << " generateParticleFate-> got empty interaction-partners list "
         << G4endl;
    return;
  }
 
  G4int npart = thePartners.size(); // Last item is a total-path placeholder
 
  if (npart == 1) {         // cparticle is on the next zone entry
    if (verboseLevel > 1) 
      G4cout << " no interactions; moving to next zone" << G4endl;
    
    cparticle.propagateAlongThePath(thePartners[0].second);
    cparticle.incrementCurrentPath(thePartners[0].second);
    boundaryTransition(cparticle);
    outgoing_cparticles.push_back(cparticle);
 
    if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
 
    //A.R. 19-Jun-2013: Fixed rare cases of non-reproducibility.
    current_nucl1 = 0;
    current_nucl2 = 0;
 
    return;
  } // if (npart == 1)
 
  if (verboseLevel > 1)
    G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
  
  G4ThreeVector old_position = cparticle.getPosition();
  G4InuclElementaryParticle& bullet = cparticle.getParticle();
  G4bool no_interaction = true;
  G4int zone = cparticle.getCurrentZone();
  
  for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder
    if (i > 0) cparticle.updatePosition(old_position); 
    
    G4InuclElementaryParticle& target = thePartners[i].first; 
    
    if (verboseLevel > 3) {
      if (target.quasi_deutron()) G4cout << " try absorption: ";
      G4cout << " target " << target.type() << " bullet " << bullet.type()
         << G4endl;
    }
    
    EPCoutput.reset();
    // Pass current (A,Z) configuration for possible recoils
    G4int massNumberCurrent = protonNumberCurrent+neutronNumberCurrent;
    theEPCollider->setNucleusState(massNumberCurrent, protonNumberCurrent);
    theEPCollider->collide(&bullet, &target, EPCoutput);
    
    // If collision failed, exit loop over partners
    if (EPCoutput.numberOfOutgoingParticles() == 0) break;
    
    if (verboseLevel > 2) {
      EPCoutput.printCollisionOutput();
#ifdef G4CASCADE_CHECK_ECONS
      balance.collide(&bullet, &target, EPCoutput);
      balance.okay();       // Do checks, but ignore result
#endif
    }
    
    // Get list of outgoing particles for evaluation
    std::vector<G4InuclElementaryParticle>& outgoing_particles = 
      EPCoutput.getOutgoingParticles();
    
    if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
    
    // Trailing effect: reject interaction at previously hit nucleon
    cparticle.propagateAlongThePath(thePartners[i].second);
    const G4ThreeVector& new_position = cparticle.getPosition();
    
    if (!passTrailing(new_position)) continue;
    collisionPts.push_back(new_position);
    
    // Sort particles according to beta (fastest first)
    std::sort(outgoing_particles.begin(), outgoing_particles.end(),
          G4ParticleLargerBeta() );
    
    if (verboseLevel > 2)
      G4cout << " adding " << outgoing_particles.size()
         << " output particles" << G4endl;
    
    // NOTE:  Embedded temporary is optimized away (no copying gets done)
    G4int nextGen = cparticle.getGeneration()+1;
    for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 
      outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
                             new_position, zone,
                             0.0, nextGen));
    }
    
    no_interaction = false;
    current_nucl1 = 0;
    current_nucl2 = 0;
    
#ifdef G4CASCADE_DEBUG_CHARGE
    {
      G4double out_charge = 0.0;
      
      for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) 
    out_charge += outgoing_particles[ip].getCharge();
      
      G4cout << " multiplicity " << outgoing_particles.size()
         << " bul type " << bullet.type()
         << " targ type " << target.type()
         << "\n initial charge "
         << bullet.getCharge() + target.getCharge() 
         << " out charge " << out_charge << G4endl;  
    }
#endif
    
    if (verboseLevel > 2)
      G4cout << " partner type " << target.type() << G4endl;
    
    if (target.nucleon()) {
      current_nucl1 = target.type();
    } else {
      if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
      current_nucl1 = (target.type() - 100) / 10;
      current_nucl2 = target.type() - 100 - 10 * current_nucl1;
    }   
    
    if (current_nucl1 == 1) {
      if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
      protonNumberCurrent--;
    } else {
      if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
      neutronNumberCurrent--;
    } 
    
    if (current_nucl2 == 1) {
      if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
      protonNumberCurrent--;
    } else if (current_nucl2 == 2) {
      if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
      neutronNumberCurrent--;
    }
    
    break;
  }     // loop over partners
  
  if (no_interaction) {         // still no interactions
    if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
    
    // For conservation checking (below), get particle before updating
    static G4ThreadLocal G4InuclElementaryParticle *prescatCP_G4MT_TLS_ = 0;
    if (!prescatCP_G4MT_TLS_) {
      prescatCP_G4MT_TLS_ = new G4InuclElementaryParticle;
      G4AutoDelete::Register(prescatCP_G4MT_TLS_);
    }
    G4InuclElementaryParticle &prescatCP = *prescatCP_G4MT_TLS_;    // Avoid memory churn
    prescatCP = cparticle.getParticle();
    
    // Last "partner" is just a total-path placeholder
    cparticle.updatePosition(old_position); 
    cparticle.propagateAlongThePath(thePartners[npart-1].second);
    cparticle.incrementCurrentPath(thePartners[npart-1].second);
    boundaryTransition(cparticle);
    outgoing_cparticles.push_back(cparticle);
    
    // Check conservation for simple scattering (ignore target nucleus!)
#ifdef G4CASCADE_CHECK_ECONS
    if (verboseLevel > 2) {
      balance.collide(&prescatCP, 0, outgoing_cparticles);
      balance.okay();       // Report violations, but don't act on them
    }
#endif
  } // if (no_interaction)
 
  return;
}

References boundaryTransition(), G4CascadeCheckBalance::collide(), G4ElementaryParticleCollider::collide(), collisionPts, current_nucl1, current_nucl2, EPCoutput, G4cerr, G4cout, G4endl, G4ThreadLocal, generateInteractionPartners(), G4InuclParticle::getCharge(), G4CascadParticle::getCurrentZone(), G4CascadParticle::getGeneration(), G4CollisionOutput::getOutgoingParticles(), G4CascadParticle::getParticle(), G4CascadParticle::getPosition(), G4CascadParticle::incrementCurrentPath(), neutronNumberCurrent, G4InuclElementaryParticle::nucleon(), G4CollisionOutput::numberOfOutgoingParticles(), G4CascadeCheckBalance::okay(), passFermi(), passTrailing(), G4CollisionOutput::printCollisionOutput(), G4CascadParticle::propagateAlongThePath(), protonNumberCurrent, G4InuclElementaryParticle::quasi_deutron(), G4AutoDelete::Register(), G4CollisionOutput::reset(), second, G4ElementaryParticleCollider::setNucleusState(), G4VCascadeCollider::setVerboseLevel(), thePartners, G4InuclElementaryParticle::type(), G4CascadParticle::updatePosition(), and verboseLevel.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ generateQuasiDeuteron()

G4InuclElementaryParticle G4NucleiModel::generateQuasiDeuteron	(	G4int	type1,
		G4int	type2,
		G4int	zone
	)		const

protected

Definition at line 668 of file G4NucleiModel.cc.

                                      {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::generateQuasiDeuteron" << G4endl;
  }
 
  // Quasideuteron consists of an unbound but associated nucleon pair
  
  // FIXME:  Why generate two separate nucleon momenta (randomly!) and
  //         add them, instead of just throwing a net momentum for the
  //         dinulceon state?  And why do I have to capture the two
  //         return values into local variables?
  G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
  G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
  G4LorentzVector dmom = mom1+mom2;
 
  G4int dtype = 0;
       if (type1*type2 == pro*pro) dtype = 111;
  else if (type1*type2 == pro*neu) dtype = 112;
  else if (type1*type2 == neu*neu) dtype = 122;
 
  return G4InuclElementaryParticle(dmom, dtype);
}

References G4cout, G4endl, generateNucleonMomentum(), G4InuclParticleNames::neu, G4InuclParticleNames::pro, and verboseLevel.

Referenced by generateInteractionPartners().

◆ getCurrentDensity()

G4double G4NucleiModel::getCurrentDensity	(	G4int	ip,
		G4int	izone
	)		const

protected

Definition at line 1397 of file G4NucleiModel.cc.

                                                                     {
  const G4double pn_spec = 1.0;     // Scale factor for pn vs. pp/nn
  //const G4double pn_spec = 0.5;
 
  G4double dens = 0.;
 
  if (ip < 100) dens = getDensity(ip,izone);
  else {    // For dibaryons, remove extra 1/volume term in density product
    switch (ip) {
    case diproton:  
      dens = getDensity(proton,izone) * getDensity(proton,izone);
      break;
    case unboundPN: 
      dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
      break;
    case dineutron:
      dens = getDensity(neutron,izone) * getDensity(neutron,izone);
      break;
    default: dens = 0.;
    }
    dens *= getVolume(izone);
  }
 
  return getRatio(ip) * dens;
}

References G4InuclParticleNames::dineutron, G4InuclParticleNames::diproton, getDensity(), getRatio(), getVolume(), G4InuclParticleNames::neutron, G4InuclParticleNames::proton, and G4InuclParticleNames::unboundPN.

Referenced by inverseMeanFreePath().

◆ getDensity()

G4double G4NucleiModel::getDensity	(	G4int	ip,
		G4int	izone
	)		const

inline

Definition at line 110 of file G4NucleiModel.hh.

                                                   {
    return nucleon_densities[ip - 1][izone];
  }

References nucleon_densities.

Referenced by getCurrentDensity(), and printModel().

◆ getFermiKinetic()

G4double G4NucleiModel::getFermiKinetic	(	G4int	ip,
		G4int	izone
	)		const

Definition at line 634 of file G4NucleiModel.cc.

                                                                   {
  G4double ekin = 0.0;
  
  if (ip < 3 && izone < number_of_zones) {  // ip for proton/neutron only
    G4double pfermi = fermi_momenta[ip - 1][izone]; 
    G4double mass = G4InuclElementaryParticle::getParticleMass(ip);
    
    ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
  }  
  return ekin;
}

References fermi_momenta, G4InuclElementaryParticle::getParticleMass(), and number_of_zones.

Referenced by worthToPropagate().

◆ getFermiMomentum()

G4double G4NucleiModel::getFermiMomentum	(	G4int	ip,
		G4int	izone
	)		const

inline

Definition at line 114 of file G4NucleiModel.hh.

                                                         {
    return fermi_momenta[ip - 1][izone];
  }

References fermi_momenta.

Referenced by generateNucleonMomentum(), and printModel().

◆ getNumberOfNeutrons()

G4int G4NucleiModel::getNumberOfNeutrons ( ) const

inline

Definition at line 147 of file G4NucleiModel.hh.

147{ return neutronNumberCurrent; }

References neutronNumberCurrent.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getNumberOfProtons()

G4int G4NucleiModel::getNumberOfProtons ( ) const

inline

Definition at line 148 of file G4NucleiModel.hh.

148{ return protonNumberCurrent; }

References protonNumberCurrent.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getNumberOfZones()

G4int G4NucleiModel::getNumberOfZones ( ) const

inline

Definition at line 141 of file G4NucleiModel.hh.

141{ return number_of_zones; }

References number_of_zones.

◆ getPotential()

G4double G4NucleiModel::getPotential	(	G4int	ip,
		G4int	izone
	)		const

inline

Definition at line 120 of file G4NucleiModel.hh.

                                                     {
    if (ip == 9 || ip < 0) return 0.0;      // Photons and leptons
    G4int ip0 = ip < 3 ? ip - 1 : 2;
    if (ip > 10 && ip < 18) ip0 = 3;
    if (ip > 20) ip0 = 4;
    return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
  }

References number_of_zones, and zone_potentials.

Referenced by boundaryTransition(), and printModel().

◆ getRadius() [1/2]

G4double G4NucleiModel::getRadius ( ) const

inline

Definition at line 131 of file G4NucleiModel.hh.

131{ return nuclei_radius; }

References nuclei_radius.

◆ getRadius() [2/2]

G4double G4NucleiModel::getRadius ( G4int izone ) const

inline

Definition at line 132 of file G4NucleiModel.hh.

                                        {
    return ( (izone<0) ? 0.
         : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
  }

References nuclei_radius, number_of_zones, and zone_radii.

◆ getRadiusUnits()

G4double G4NucleiModel::getRadiusUnits ( ) const

inline

Definition at line 129 of file G4NucleiModel.hh.

129{ return radiusUnits*CLHEP::fermi; }

CLHEP::fermi

static constexpr double fermi

Definition: SystemOfUnits.h:84

References CLHEP::fermi, and radiusUnits.

Referenced by G4IntraNucleiCascader::processSecondary().

◆ getRatio()

G4double G4NucleiModel::getRatio ( G4int ip ) const

protected

Definition at line 1380 of file G4NucleiModel.cc.

                                               {
  if (verboseLevel > 4) {
    G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
  }
 
  switch (ip) {
  case proton:    return G4double(protonNumberCurrent)/G4double(protonNumber);
  case neutron:   return G4double(neutronNumberCurrent)/G4double(neutronNumber);
  case diproton:  return getRatio(proton)*getRatio(proton);
  case unboundPN: return getRatio(proton)*getRatio(neutron);
  case dineutron: return getRatio(neutron)*getRatio(neutron);
  default:        return 0.;
  }
 
  return 0.;
}

References G4InuclParticleNames::dineutron, G4InuclParticleNames::diproton, G4cout, G4endl, getRatio(), G4InuclParticleNames::neutron, neutronNumber, neutronNumberCurrent, G4InuclParticleNames::proton, protonNumber, protonNumberCurrent, G4InuclParticleNames::unboundPN, and verboseLevel.

Referenced by getCurrentDensity(), and getRatio().

◆ getTypesOfNucleonsInvolved()

std::pair< G4int, G4int > G4NucleiModel::getTypesOfNucleonsInvolved ( ) const

inline

Definition at line 167 of file G4NucleiModel.hh.

                                                           {
    return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
  }

References current_nucl1, and current_nucl2.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ getVolume()

G4double G4NucleiModel::getVolume ( G4int izone ) const

inline

Definition at line 136 of file G4NucleiModel.hh.

                                        {
    return ( (izone<0) ? 0.
         : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
  }

References nuclei_volume, number_of_zones, and zone_volumes.

Referenced by getCurrentDensity().

◆ getZone()

G4int G4NucleiModel::getZone ( G4double r ) const

inline

Definition at line 142 of file G4NucleiModel.hh.

                                  {
    for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
    return number_of_zones;
  }

References number_of_zones, and zone_radii.

Referenced by choosePointAlongTraj(), and G4IntraNucleiCascader::processSecondary().

◆ initializeCascad() [1/2]

G4CascadParticle G4NucleiModel::initializeCascad ( G4InuclElementaryParticle * particle )

Definition at line 1425 of file G4NucleiModel.cc.

                                                                   {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
  }
 
  // FIXME:  Previous version generated random sin(theta), then used -cos(theta)
  //         Using generateWithRandomAngles changes result!
  // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
  G4double costh = std::sqrt(1.0 - inuclRndm());
  G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
 
  // Start particle outside nucleus, unless capture-at-rest
  G4int zone = number_of_zones;
  if (particle->getKineticEnergy() < small) zone--;
 
  G4CascadParticle cpart(*particle, pos, zone, large, 0);
 
  // SPECIAL CASE:  Inbound photons are emplanted along through-path
  if (forceFirst(cpart)) choosePointAlongTraj(cpart);
 
  if (verboseLevel > 2) G4cout << cpart << G4endl;
 
  return cpart;
}

References choosePointAlongTraj(), forceFirst(), G4cout, G4endl, G4InuclSpecialFunctions::generateWithFixedTheta(), G4InuclParticle::getKineticEnergy(), G4InuclSpecialFunctions::inuclRndm(), large, nuclei_radius, number_of_zones, pos, small, and verboseLevel.

Referenced by G4IntraNucleiCascader::setupCascade().

◆ initializeCascad() [2/2]

void G4NucleiModel::initializeCascad	(	G4InuclNuclei *	bullet,
		G4InuclNuclei *	target,
		modelLists &	output
	)

Definition at line 1450 of file G4NucleiModel.cc.

                                         {
  if (verboseLevel) {
    G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
       << G4endl;
  }
 
  const G4double max_a_for_cascad = 5.0;
  const G4double ekin_cut = 2.0;
  const G4double small_ekin = 1.0e-6;
  const G4double r_large2for3 = 62.0;
  const G4double r0forAeq3 = 3.92;
  const G4double s3max = 6.5;
  const G4double r_large2for4 = 69.14;
  const G4double r0forAeq4 = 4.16;
  const G4double s4max = 7.0;
  const G4int itry_max = 100;
 
  // Convenient references to input buffer contents
  std::vector<G4CascadParticle>& casparticles = output.first;
  std::vector<G4InuclElementaryParticle>& particles = output.second;
 
  casparticles.clear();     // Reset input buffer to fill this event
  particles.clear();
 
  // first decide whether it will be cascad or compound final nuclei
  G4int ab = bullet->getA();
  G4int zb = bullet->getZ();
  G4int at = target->getA();
  G4int zt = target->getZ();
 
  G4double massb = bullet->getMass();   // For creating LorentzVectors below
 
  if (ab < max_a_for_cascad) {
 
    G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
    G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
    G4double ben = benb < bent ? bent : benb;
 
    if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
      G4int itryg = 0;
 
    /* Loop checking 08.06.2015 MHK */
      while (casparticles.size() == 0 && itryg < itry_max) {      
    itryg++;
    particles.clear();
      
    //    nucleons coordinates and momenta in nuclei rest frame
    coordinates.clear();
    momentums.clear();
     
    if (ab < 3) { // deuteron, simplest case
      G4double r = 2.214 - 3.4208 * G4Log(1.0 - 0.981 * inuclRndm());
      G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
      coordinates.push_back(coord1);
      coordinates.push_back(-coord1);
 
      G4double p = 0.0;
      G4bool bad = true;
      G4int itry = 0;
 
      while (bad && itry < itry_max) {  /* Loop checking 08.06.2015 MHK */
        itry++;
        p = 456.0 * inuclRndm();
 
        if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
           p * r > 312.0) bad = false;
      }
 
      if (itry == itry_max)
        if (verboseLevel > 2){ 
          G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;   
        }
 
      p = 0.0005 * p;
 
      if (verboseLevel > 2){ 
        G4cout << " p nuc " << p << G4endl;
      }
 
      G4LorentzVector mom = generateWithRandomAngles(p, massb);
 
      momentums.push_back(mom);
      mom.setVect(-mom.vect());
      momentums.push_back(-mom);
    } else {
      G4ThreeVector coord1;
 
      G4bool badco = true;
 
      G4int itry = 0;
        
      if (ab == 3) {
        while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
          if (itry > 0) coordinates.clear();
          itry++;   
          G4int i(0);    
 
          for (i = 0; i < 2; i++) {
        G4int itry1 = 0;
        G4double ss, u, rho; 
        G4double fmax = G4Exp(-0.5) / std::sqrt(0.5);
 
        while (itry1 < itry_max) {  /* Loop checking 08.06.2015 MHK */
          itry1++;
          ss = -G4Log(inuclRndm());
          u = fmax * inuclRndm();
          rho = std::sqrt(ss) * G4Exp(-ss);
 
          if (rho > u && ss < s3max) {
            ss = r0forAeq3 * std::sqrt(ss);
            coord1 = generateWithRandomAngles(ss).vect();
            coordinates.push_back(coord1);
 
            if (verboseLevel > 2){
              G4cout << " i " << i << " r " << coord1.mag() << G4endl;
            }
            break;
          }
        }
 
        if (itry1 == itry_max) { // bad case
          coord1.set(10000.,10000.,10000.);
          coordinates.push_back(coord1);
          break;
        }
          }
 
          coord1 = -coordinates[0] - coordinates[1]; 
          if (verboseLevel > 2) {
        G4cout << " 3  r " << coord1.mag() << G4endl;
          }
 
          coordinates.push_back(coord1);        
        
          G4bool large_dist = false;
 
          for (i = 0; i < 2; i++) {
        for (G4int j = i+1; j < 3; j++) {
          G4double r2 = (coordinates[i]-coordinates[j]).mag2();
 
          if (verboseLevel > 2) {
            G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
          }
 
          if (r2 > r_large2for3) {
            large_dist = true;
 
            break; 
          }      
        }
 
        if (large_dist) break;
          } 
 
          if(!large_dist) badco = false;
 
        }
 
      } else { // a >= 4
        G4double b = 3./(ab - 2.0);
        G4double b1 = 1.0 - b / 2.0;
        G4double u = b1 + std::sqrt(b1 * b1 + b);
        G4double fmax = (1.0 + u/b) * u * G4Exp(-u);
      
        while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
 
          if (itry > 0) coordinates.clear();
          itry++;
          G4int i(0);
        
          for (i = 0; i < ab-1; i++) {
        G4int itry1 = 0;
        G4double ss; 
 
        while (itry1 < itry_max) {  /* Loop checking 08.06.2015 MHK */
          itry1++;
          ss = -G4Log(inuclRndm());
          u = fmax * inuclRndm();
 
          if (std::sqrt(ss) * G4Exp(-ss) * (1.0 + ss/b) > u
              && ss < s4max) {
            ss = r0forAeq4 * std::sqrt(ss);
            coord1 = generateWithRandomAngles(ss).vect();
            coordinates.push_back(coord1);
 
            if (verboseLevel > 2) {
              G4cout << " i " << i << " r " << coord1.mag() << G4endl;
            }
 
            break;
          }
        }
 
        if (itry1 == itry_max) { // bad case
          coord1.set(10000.,10000.,10000.);
          coordinates.push_back(coord1);
          break;
        }
          }
 
          coord1 *= 0.0;    // Cheap way to reset
          for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
 
          coordinates.push_back(coord1);   
 
          if (verboseLevel > 2){
        G4cout << " last r " << coord1.mag() << G4endl;
          }
        
          G4bool large_dist = false;
 
          for (i = 0; i < ab-1; i++) {
        for (G4int j = i+1; j < ab; j++) {
         
          G4double r2 = (coordinates[i]-coordinates[j]).mag2();
 
          if (verboseLevel > 2){
            G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
          }
 
          if (r2 > r_large2for4) {
            large_dist = true;
 
            break; 
          }      
        }
 
        if (large_dist) break;
          } 
 
          if (!large_dist) badco = false;
        }
      } 
 
      if(badco) {
        G4cout << " can not generate the nucleons coordinates for a "
           << ab << G4endl; 
 
        casparticles.clear();   // Return empty buffer on error
        particles.clear();
        return;
 
      } else { // momentums
        G4double p, u, x;
        G4LorentzVector mom;
        //G4bool badp = True;
 
        for (G4int i = 0; i < ab - 1; i++) {
          G4int itry2 = 0;
 
          while(itry2 < itry_max) { /* Loop checking 08.06.2015 MHK */
        itry2++;
        u = -G4Log(0.879853 - 0.8798502 * inuclRndm());
        x = u * G4Exp(-u);
 
        if(x > inuclRndm()) {
          p = std::sqrt(0.01953 * u);
          mom = generateWithRandomAngles(p, massb);
          momentums.push_back(mom);
 
          break;
        }
          } // while (itry2 < itry_max)
 
          if(itry2 == itry_max) {
        G4cout << " can not generate proper momentum for a "
               << ab << G4endl;
 
        casparticles.clear();   // Return empty buffer on error
        particles.clear();
        return;
          } 
        }   // for (i=0 ...
        // last momentum
 
        mom *= 0.;      // Cheap way to reset
        mom.setE(bullet->getEnergy()+target->getEnergy());
 
        for(G4int j=0; j< ab-1; j++) mom -= momentums[j]; 
 
        momentums.push_back(mom);
      } 
    }
 
    // Coordinates and momenta at rest are generated, now back to the lab
    G4double rb = 0.0;
    G4int i(0);
 
    for(i = 0; i < G4int(coordinates.size()); i++) {      
      G4double rp = coordinates[i].mag2();
 
      if(rp > rb) rb = rp;
    }
 
    // nuclei i.p. as a whole
    G4double s1 = std::sqrt(inuclRndm()); 
    G4double phi = randomPHI();
    G4double rz = (nuclei_radius + rb) * s1;
    G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
                 -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
 
    for (i = 0; i < G4int(coordinates.size()); i++) {
      coordinates[i] += global_pos;
    }  
 
    // all nucleons at rest
    raw_particles.clear();
 
    for (G4int ipa = 0; ipa < ab; ipa++) {
      G4int knd = ipa < zb ? 1 : 2;
      raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
    } 
      
    G4InuclElementaryParticle dummy(small_ekin, 1);
    G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
    toTheBulletRestFrame.toTheTargetRestFrame();
 
    particleIterator ipart;
 
    for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
      ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 
    }
 
    // fill cascad particles and outgoing particles
 
    for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
      G4LorentzVector mom = raw_particles[ip].getMomentum();
      G4double pmod = mom.rho();
      G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
      G4double det = t0 * t0 + nuclei_radius * nuclei_radius
                   - coordinates[ip].mag2();
      G4double tr = -1.0;
 
      if(det > 0.0) {
        G4double t1 = t0 + std::sqrt(det);
        G4double t2 = t0 - std::sqrt(det);
 
        if(std::fabs(t1) <= std::fabs(t2)) {     
          if(t1 > 0.0) {
        if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
          }
 
          if(tr < 0.0 && t2 > 0.0) {
 
        if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
          }
 
        } else {
          if(t2 > 0.0) {
 
        if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
          }
 
          if(tr < 0.0 && t1 > 0.0) {
        if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
          }
        } 
 
      }
 
      if(tr >= 0.0) { // cascad particle
        coordinates[ip] += mom.vect()*tr / pmod;
        casparticles.push_back(G4CascadParticle(raw_particles[ip], 
                            coordinates[ip], 
                            number_of_zones, large, 0));
 
      } else {
        particles.push_back(raw_particles[ip]); 
      } 
    }
      }    
 
      if(casparticles.size() == 0) {
    particles.clear();
 
    G4cout << " can not generate proper distribution for " << itry_max
           << " steps " << G4endl;
      }    
    }
  }
 
  if(verboseLevel > 2){
    G4int ip(0);
 
    G4cout << " cascad particles: " << casparticles.size() << G4endl;
    for(ip = 0; ip < G4int(casparticles.size()); ip++)
      G4cout << casparticles[ip] << G4endl;
 
    G4cout << " outgoing particles: " << particles.size() << G4endl;
    for(ip = 0; ip < G4int(particles.size()); ip++)
      G4cout << particles[ip] << G4endl;
  }
 
  return;   // Buffer has been filled
}

References ab, G4LorentzConvertor::backToTheLab(), G4InuclSpecialFunctions::bindingEnergy(), coordinates, CLHEP::Hep3Vector::dot(), G4cout, G4endl, G4Exp(), G4Log(), G4InuclSpecialFunctions::generateWithRandomAngles(), G4InuclNuclei::getA(), G4InuclParticle::getEnergy(), G4InuclParticle::getKineticEnergy(), G4InuclParticle::getMass(), G4InuclNuclei::getZ(), GeV, G4InuclSpecialFunctions::inuclRndm(), large, CLHEP::Hep3Vector::mag(), momentums, nuclei_radius, number_of_zones, G4InuclSpecialFunctions::randomPHI(), raw_particles, CLHEP::HepLorentzVector::rho(), CLHEP::Hep3Vector::set(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setVect(), G4LorentzConvertor::toTheTargetRestFrame(), CLHEP::HepLorentzVector::vect(), verboseLevel, and CLHEP::HepLorentzVector::z().

◆ inverseMeanFreePath()

G4double G4NucleiModel::inverseMeanFreePath	(	const G4CascadParticle &	cparticle,
		const G4InuclElementaryParticle &	target,
		G4int	zone = `-1`
	)

protected

Definition at line 1851 of file G4NucleiModel.cc.

                               {
  G4int ptype = cparticle.getParticle().type();
  G4int ip = target.type();
 
  // Ensure that zone specified is within nucleus, for array lookups
  if (zone<0) zone = cparticle.getCurrentZone();
  if (zone>=number_of_zones) zone = number_of_zones-1;
 
  // Special cases: neutrinos, and muon-on-neutron, have infinite path
  if (cparticle.getParticle().isNeutrino()) return 0.;
  if (ptype == muonMinus && ip == neutron) return 0.;
 
  // Configure frame transformation to get kinetic energy for lookups
  dummy_convertor.setBullet(cparticle.getParticle());
  dummy_convertor.setTarget(&target);
  dummy_convertor.toTheCenterOfMass();      // Fill internal kinematics
  G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
 
  // Get cross section for interacting with target (dibaryons are absorptive)
  G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
                             : absorptionCrossSection(ekin, ptype);
 
  if (verboseLevel > 2) {
    G4cout << " ip " << ip << " zone " << zone << " ekin " << ekin
       << " dens " << getCurrentDensity(ip, zone)
       << " csec " << csec << G4endl;
  }
 
  if (csec <= 0.) return 0.;    // No interaction, avoid unnecessary work
 
  // Get nuclear density and compute mean free path
  return csec * getCurrentDensity(ip, zone);
}

References absorptionCrossSection(), dummy_convertor, G4cout, G4endl, getCurrentDensity(), G4CascadParticle::getCurrentZone(), G4LorentzConvertor::getKinEnergyInTheTRS(), G4CascadParticle::getParticle(), G4InuclElementaryParticle::isNeutrino(), G4InuclParticleNames::muonMinus, G4InuclParticleNames::neutron, number_of_zones, G4LorentzConvertor::setBullet(), G4LorentzConvertor::setTarget(), totalCrossSection(), G4LorentzConvertor::toTheCenterOfMass(), G4InuclElementaryParticle::type(), and verboseLevel.

Referenced by choosePointAlongTraj(), and generateInteractionPartners().

◆ isProjectile()

G4bool G4NucleiModel::isProjectile ( const G4CascadParticle & cparticle ) const

Definition at line 1345 of file G4NucleiModel.cc.

                                                                          {
  return (cparticle.getGeneration() == 0);  // Only initial state particles
}

References G4CascadParticle::getGeneration().

Referenced by forceFirst(), and generateInteractionPartners().

◆ passFermi()

G4bool G4NucleiModel::passFermi	(	const std::vector< G4InuclElementaryParticle > &	particles,
		G4int	zone
	)

protected

Definition at line 1069 of file G4NucleiModel.cc.

                            {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
  }
 
  // Only check Fermi momenta for nucleons
  for (G4int i = 0; i < G4int(particles.size()); i++) {
    if (!particles[i].nucleon()) continue;
 
    G4int type   = particles[i].type();
    G4double mom = particles[i].getMomModule();
    G4double pfermi  = fermi_momenta[type-1][zone];
 
    if (verboseLevel > 2)
      G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
    
    if (mom < pfermi) {
    if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
    return false;
    }
  }
  return true; 
}

References fermi_momenta, G4cout, G4endl, G4InuclParticleNames::nucleon(), and verboseLevel.

Referenced by generateParticleFate().

◆ passTrailing()

G4bool G4NucleiModel::passTrailing ( const G4ThreeVector & hit_position )

protected

Definition at line 1098 of file G4NucleiModel.cc.

                                                                    {
  if (verboseLevel > 1)
    G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
 
  G4double dist;
  for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
    dist = (collisionPts[i] - hit_position).mag();
    if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
    if (dist < R_nucleon) {
      if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
      return false;
    }
  }
  return true;      // New point far enough away to be used
}

References collisionPts, G4cout, G4endl, R_nucleon, and verboseLevel.

Referenced by generateParticleFate().

◆ printModel()

void G4NucleiModel::printModel ( ) const

Definition at line 612 of file G4NucleiModel.cc.

                                     {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::printModel" << G4endl;
  }
 
  G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
     << " proton binding energy " << binding_energies[0]
     << " neutron binding energy " << binding_energies[1] << G4endl
     << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
     << " number of zones " << number_of_zones << G4endl;
 
  for (G4int i = 0; i < number_of_zones; i++)
    G4cout << " zone " << i+1 << " radius " << zone_radii[i] 
       << " volume " << zone_volumes[i] << G4endl
       << " protons: density " << getDensity(1,i) << " PF " << 
      getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
       << " neutrons: density " << getDensity(2,i) << " PF " << 
      getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
       << " pions: VP " << getPotential(3,i) << G4endl;
}

References A, binding_energies, G4cout, G4endl, getDensity(), getFermiMomentum(), getPotential(), nuclei_radius, nuclei_volume, number_of_zones, verboseLevel, Z, zone_radii, and zone_volumes.

Referenced by generateModel().

◆ reset()

void G4NucleiModel::reset	(	G4int	nHitNeutrons = `0`,
		G4int	nHitProtons = `0`,
		const std::vector< G4ThreeVector > *	hitPoints = `0`
	)

Definition at line 298 of file G4NucleiModel.cc.

                                                         {
  neutronNumberCurrent = neutronNumber - nHitNeutrons;
  protonNumberCurrent  = protonNumber - nHitProtons;
  
  // zero or copy collision point array for trailing effect
  if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
  else collisionPts = *hitPoints;
}

References collisionPts, neutronNumber, neutronNumberCurrent, protonNumber, and protonNumberCurrent.

Referenced by G4IntraNucleiCascader::copyWoundedNucleus(), generateModel(), and G4IntraNucleiCascader::newCascade().

◆ setVerboseLevel()

void G4NucleiModel::setVerboseLevel ( G4int verbose )

inline

Definition at line 99 of file G4NucleiModel.hh.

99{ verboseLevel = verbose; }

References verboseLevel.

Referenced by G4IntraNucleiCascader::setVerboseLevel().

◆ sortPartners()

static G4bool G4NucleiModel::sortPartners	(	const partner &	p1,
		const partner &	p2
	)

inlinestaticprotected

Definition at line 209 of file G4NucleiModel.hh.

                                                                   {
    return (p2.second > p1.second);
  }

Referenced by generateInteractionPartners().

◆ stillInside()

G4bool G4NucleiModel::stillInside ( const G4CascadParticle & cparticle )

inline

Definition at line 154 of file G4NucleiModel.hh.

                                                        {
    return cparticle.getCurrentZone() < number_of_zones;
  }

References G4CascadParticle::getCurrentZone(), and number_of_zones.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ totalCrossSection()

G4double G4NucleiModel::totalCrossSection	(	G4double	ke,
		G4int	rtype
	)		const

Definition at line 1954 of file G4NucleiModel.cc.

{
  // All scattering cross-sections are available from tables
  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
  if (!xsecTable) {
    G4cerr << " unknown collison type = " << rtype << G4endl;
    return 0.;
  }
 
  return (crossSectionUnits * xsecTable->getCrossSection(ke));
}

References crossSectionUnits, G4cerr, G4endl, G4CascadeChannel::getCrossSection(), and G4CascadeChannelTables::GetTable().

Referenced by inverseMeanFreePath().

◆ useQuasiDeuteron()

G4bool G4NucleiModel::useQuasiDeuteron	(	G4int	ptype,
		G4int	qdtype = `0`
	)

static

Definition at line 1057 of file G4NucleiModel.cc.

                                                                {
  if (qdtype==pn || qdtype==0)      // All absorptive particles
    return (ptype==pi0 || ptype==pip || ptype==pim || ptype==gam || ptype==mum);
  else if (qdtype==pp)          // Negative or neutral only
    return (ptype==pi0 || ptype==pim || ptype==gam || ptype==mum);
  else if (qdtype==nn)          // Positive or neutral only
    return (ptype==pi0 || ptype==pip || ptype==gam);
 
  return false;     // Not a quasideuteron target
}

References G4InuclParticleNames::gam, G4InuclParticleNames::mum, G4InuclParticleNames::nn, G4InuclParticleNames::pi0, G4InuclParticleNames::pim, G4InuclParticleNames::pip, G4InuclParticleNames::pn, and G4InuclParticleNames::pp.

Referenced by absorptionCrossSection(), G4ElementaryParticleCollider::collide(), and generateInteractionPartners().

◆ worthToPropagate()

G4bool G4NucleiModel::worthToPropagate ( const G4CascadParticle & cparticle ) const

Definition at line 1349 of file G4NucleiModel.cc.

                                                                              {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
  }
 
  const G4double ekin_scale = 2.0;
 
  G4bool worth = true;
 
  if (cparticle.reflectedNow()) {   // Just reflected -- keep going?
    G4int zone = cparticle.getCurrentZone();
    G4int ip = cparticle.getParticle().type();
 
    // NOTE:  Temporarily backing out use of potential for non-nucleons
    G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
      getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
 
    worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
 
    if (verboseLevel > 3) {
      G4cout << " type=" << ip
         << " ekin=" << cparticle.getParticle().getKineticEnergy()
         << " potential=" << ekin_cut
         << " : worth? " << worth << G4endl;
    }
  }
 
  return worth;
}

References G4cout, G4endl, G4CascadParticle::getCurrentZone(), getFermiKinetic(), G4InuclParticle::getKineticEnergy(), G4CascadParticle::getParticle(), G4InuclElementaryParticle::nucleon(), G4CascadParticle::reflectedNow(), G4InuclElementaryParticle::type(), and verboseLevel.

Referenced by G4IntraNucleiCascader::generateCascade().

◆ zoneIntegralGaussian()

G4double G4NucleiModel::zoneIntegralGaussian	(	G4double	ur1,
		G4double	ur2,
		G4double	nuclearRadius
	)		const

protected

Definition at line 564 of file G4NucleiModel.cc.

                                                {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
  }
 
  G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
 
  const G4double epsilon = 1.0e-3;
  const G4int itry_max = 1000;
 
  G4double dr = r2 - r1;
  G4double fr1 = r1 * r1 * G4Exp(-r1 * r1);
  G4double fr2 = r2 * r2 * G4Exp(-r2 * r2);
  G4double fi = (fr1 + fr2) / 2.;
  G4double fun1 = fi * dr;
  G4double fun;
  G4int jc = 1;
  G4double dr1 = dr;
  G4int itry = 0;
 
  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
    dr /= 2.;
    itry++;
    G4double r = r1 - dr;
    fi = 0.0;
 
    for (G4int i = 0; i < jc; i++) { 
      r += dr1; 
      fi += r * r * G4Exp(-r * r);
    }
 
    fun = 0.5 * fun1 + fi * dr;  
 
    if (std::fabs((fun - fun1) / fun) <= epsilon) break;
 
    jc *= 2;
    dr1 = dr;
    fun1 = fun;
  } // while (itry < itry_max)
 
  if (verboseLevel > 2 && itry == itry_max)
    G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
 
  return gaussRadius*gaussRadius*gaussRadius * fun;
}

References A, epsilon(), G4cerr, G4cout, G4endl, G4Exp(), and verboseLevel.

Referenced by fillZoneVolumes().

◆ zoneIntegralWoodsSaxon()

G4double G4NucleiModel::zoneIntegralWoodsSaxon	(	G4double	ur1,
		G4double	ur2,
		G4double	nuclearRadius
	)		const

protected

Definition at line 511 of file G4NucleiModel.cc.

                                                         {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
  }
 
  const G4double epsilon = 1.0e-3;
  const G4int itry_max = 1000;
 
  G4double skinRatio = nuclearRadius / skinDepth;
 
  G4double d2 = 2.0 * skinRatio;
  G4double dr = r2 - r1;
  G4double fr1 = r1 * (r1 + d2) / (1.0 + G4Exp(r1));
  G4double fr2 = r2 * (r2 + d2) / (1.0 + G4Exp(r2));
  G4double fi = (fr1 + fr2) / 2.;
  G4double fun1 = fi * dr;
  G4double fun;
  G4int jc = 1;
  G4double dr1 = dr;
  G4int itry = 0;
 
  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
    dr /= 2.;
    itry++;
 
    G4double r = r1 - dr;
    fi = 0.0;
 
    for (G4int i = 0; i < jc; i++) { 
      r += dr1; 
      fi += r * (r + d2) / (1.0 + G4Exp(r));
    }
 
    fun = 0.5 * fun1 + fi * dr;
 
    if (std::fabs((fun - fun1) / fun) <= epsilon) break;
 
    jc *= 2;
    dr1 = dr;
    fun1 = fun;
  } // while (itry < itry_max)
 
  if (verboseLevel > 2 && itry == itry_max)
    G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
 
  G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
 
  return skinDepth3 * (fun + skinRatio*skinRatio*G4Log((1.0 + G4Exp(-r1)) / (1.0 + G4Exp(-r2))));
}