#include <G4NucleiModel.hh>

Collaboration diagram for G4NucleiModel:

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

Public Member Functions
	G4NucleiModel ()

	G4NucleiModel (G4int a, G4int z)

	G4NucleiModel (G4InuclNuclei *nuclei)

virtual	~G4NucleiModel ()

void	setVerboseLevel (G4int verbose)

void	generateModel (G4InuclNuclei *nuclei)

void	generateModel (G4int a, G4int z)

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

void	printModel () const

G4double	getDensity (G4int ip, G4int izone) const

G4double	getFermiMomentum (G4int ip, G4int izone) const

G4double	getFermiKinetic (G4int ip, G4int izone) const

G4double	getPotential (G4int ip, G4int izone) const

G4double	getRadiusUnits () const

G4double	getRadius () const

G4double	getRadius (G4int izone) const

G4double	getVolume (G4int izone) const

G4int	getNumberOfZones () const

G4int	getZone (G4double r) const

G4int	getNumberOfNeutrons () const

G4int	getNumberOfProtons () const

G4bool	empty () const

G4bool	stillInside (const G4CascadParticle &cparticle)

G4CascadParticle	initializeCascad (G4InuclElementaryParticle *particle)

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

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

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

G4bool	forceFirst (const G4CascadParticle &cparticle) const

G4bool	isProjectile (const G4CascadParticle &cparticle) const

G4bool	worthToPropagate (const G4CascadParticle &cparticle) const

G4InuclElementaryParticle	generateNucleon (G4int type, G4int zone) const

G4LorentzVector	generateNucleonMomentum (G4int type, G4int zone) const

G4double	absorptionCrossSection (G4double e, G4int type) const

G4double	totalCrossSection (G4double ke, G4int rtype) const

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

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

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

G4bool	passTrailing (const G4ThreeVector &hit_position)

void	boundaryTransition (G4CascadParticle &cparticle)

void	choosePointAlongTraj (G4CascadParticle &cparticle)

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

void	generateInteractionPartners (G4CascadParticle &cparticle)

void	fillBindingEnergies ()

void	fillZoneRadii (G4double nuclearRadius)

G4double	fillZoneVolumes (G4double nuclearRadius)

void	fillPotentials (G4int type, G4double tot_vol)

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

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

G4double	getRatio (G4int ip) const

G4double	getCurrentDensity (G4int ip, G4int izone) const

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

G4double	generateInteractionLength (const G4CascadParticle &cparticle, G4double path, G4double invmfp) 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	verboseLevel

G4LorentzConvertor	dummy_convertor

G4CollisionOutput	EPCoutput

std::vector< G4InuclElementaryParticle >	qdeutrons

std::vector< G4double >	acsecs

std::vector< G4ThreeVector >	coordinates

std::vector< G4LorentzVector >	momentums

std::vector< G4InuclElementaryParticle >	raw_particles

std::vector< G4ThreeVector >	collisionPts

G4double	ur [7]

G4double	v [6]

G4double	v1 [6]

std::vector< G4double >	rod

std::vector< G4double >	pf

std::vector< G4double >	vz

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

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

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

std::vector< G4double >	zone_radii

std::vector< G4double >	zone_volumes

std::vector< G4double >	binding_energies

G4double	nuclei_radius

G4double	nuclei_volume

G4int	number_of_zones

G4int	A

G4int	Z

G4InuclNuclei *	theNucleus

G4int	neutronNumber

G4int	protonNumber

G4int	neutronNumberCurrent

G4int	protonNumberCurrent

G4int	current_nucl1

G4int	current_nucl2

G4CascadeInterpolator< 30 >	gammaQDinterp

const G4double	crossSectionUnits

const G4double	radiusUnits

const G4double	skinDepth

const G4double	radiusScale

const G4double	radiusScale2

const G4double	radiusForSmall

const G4double	radScaleAlpha

const G4double	fermiMomentum

const G4double	R_nucleon

const G4double	gammaQDscale

const G4InuclElementaryParticle	neutronEP

const G4InuclElementaryParticle	protonEP

Static Private Attributes
static const G4double	small = 1.0e-9

static const G4double	large = 1000.

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

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	pion_vp = 0.007

static const G4double	pion_vp_small = 0.007

static const G4double	kaon_vp = 0.015

static const G4double	hyperon_vp = 0.030

Detailed Description

Definition at line 93 of file G4NucleiModel.hh.

Member Typedef Documentation

◆ modelLists

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

Definition at line 163 of file G4NucleiModel.hh.

◆ partner

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

protected

Definition at line 205 of file G4NucleiModel.hh.

Member Enumeration Documentation

◆ PotentialType

enum G4NucleiModel::PotentialType

private

Enumerator
WoodsSaxon
Gaussian

Definition at line 297 of file G4NucleiModel.hh.

297 { WoodsSaxon=0, Gaussian=1 };

G4NucleiModel::WoodsSaxon

Definition: G4NucleiModel.hh:297

G4NucleiModel::Gaussian

Definition: G4NucleiModel.hh:297

Constructor & Destructor Documentation

◆ G4NucleiModel() [1/3]

G4NucleiModel::G4NucleiModel ( )

Definition at line 227 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()),
     neutronEP(neutron), protonEP(proton) {}

◆ G4NucleiModel() [2/3]

G4NucleiModel::G4NucleiModel	(	G4int	a,
		G4int	z
	)

Definition at line 244 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()),
     neutronEP(neutron), protonEP(proton) {
   generateModel(a,z);
 }

Here is the call graph for this function:

◆ G4NucleiModel() [3/3]

G4NucleiModel::G4NucleiModel ( G4InuclNuclei * nuclei )

explicit

Definition at line 263 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()),
     neutronEP(neutron), protonEP(proton) {
   generateModel(nuclei);
 }

Here is the call graph for this function:

◆ ~G4NucleiModel()

G4NucleiModel::~G4NucleiModel ( )

virtual

Definition at line 282 of file G4NucleiModel.cc.

                               {
   delete theNucleus;
   theNucleus = 0;
 }

Member Function Documentation

◆ absorptionCrossSection()

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

Definition at line 1885 of file G4NucleiModel.cc.

                                                                             {
   if (!useQuasiDeuteron(type)) {
     G4cerr << " absorptionCrossSection only valid for incident pions" << 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);     
   }
 
   // Photon cross-section is binned from parametrization by Mi. Kossov
   // See below for initialization of gammaQDinterp, gammaQDxsec
   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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ boundaryTransition()

void G4NucleiModel::boundaryTransition ( G4CascadParticle & cparticle )

protected

Definition at line 1107 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 pr = pos.dot(mom.vect()) / r;
   
   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 p1r = 0.;
   
   if (verboseLevel > 3) {
     G4cout << " type " << type << " zone " << zone << " next " << next_zone
        << " qv " << qv << " dv " << dv << G4endl;
   }
 
   if (qv <= 0.0) {  // reflection
     if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
     p1r = -pr;
     cparticle.incrementReflectionCounter();
   } else {      // transition
     if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
     p1r = std::sqrt(qv);
     if (pr < 0.0) p1r = -p1r;
 
     cparticle.updateZone(next_zone);
     cparticle.resetReflection();
   }
   
   G4double prr = (p1r - pr) / 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;
   }
 
   mom.setVect(mom.vect() + pos*prr);
   cparticle.updateParticleMomentum(mom);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ choosePointAlongTraj()

void G4NucleiModel::choosePointAlongTraj ( G4CascadParticle & cparticle )

protected

Definition at line 1175 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::bind2nd(std::divides<G4double>(), 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;
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ empty()

G4bool G4NucleiModel::empty ( ) const

inline

Definition at line 152 of file G4NucleiModel.hh.

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

Here is the caller graph for this function:

◆ fillBindingEnergies()

void G4NucleiModel::fillBindingEnergies ( )

protected

Definition at line 382 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);
 }

Here is the caller graph for this function:

◆ fillPotentials()

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

protected

Definition at line 471 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ fillZoneRadii()

void G4NucleiModel::fillZoneRadii ( G4double nuclearRadius )

protected

Definition at line 396 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;
     }
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ fillZoneVolumes()

G4double G4NucleiModel::fillZoneVolumes ( G4double nuclearRadius )

protected

Definition at line 436 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
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ forceFirst()

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

Definition at line 1301 of file G4NucleiModel.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateInteractionLength()

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

protected

Definition at line 1853 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateInteractionPartners()

void G4NucleiModel::generateInteractionPartners ( G4CascadParticle & cparticle )

protected

Definition at line 686 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));
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateModel() [1/2]

void G4NucleiModel::generateModel ( G4InuclNuclei * nuclei )

Definition at line 303 of file G4NucleiModel.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateModel() [2/2]

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

Definition at line 307 of file G4NucleiModel.cc.

                                                   {
   if (verboseLevel) {
     G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
        << G4endl;
   }
 
   // If model already built, just return; otherwise intialize 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();
 }

Here is the call graph for this function:

◆ generateNucleon()

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

Definition at line 649 of file G4NucleiModel.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateNucleonMomentum()

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

Definition at line 640 of file G4NucleiModel.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateParticleFate()

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

Definition at line 857 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ generateQuasiDeuteron()

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

protected

Definition at line 660 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getCurrentDensity()

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

protected

Definition at line 1360 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getDensity()

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

inline

Definition at line 112 of file G4NucleiModel.hh.

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

Here is the caller graph for this function:

◆ getFermiKinetic()

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

Definition at line 626 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getFermiMomentum()

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

inline

Definition at line 116 of file G4NucleiModel.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getNumberOfNeutrons()

G4int G4NucleiModel::getNumberOfNeutrons ( ) const

inline

Definition at line 149 of file G4NucleiModel.hh.

149 { return neutronNumberCurrent; }

G4NucleiModel::neutronNumberCurrent

G4int neutronNumberCurrent

Definition: G4NucleiModel.hh:288

Here is the caller graph for this function:

◆ getNumberOfProtons()

G4int G4NucleiModel::getNumberOfProtons ( ) const

inline

Definition at line 150 of file G4NucleiModel.hh.

150 { return protonNumberCurrent; }

G4NucleiModel::protonNumberCurrent

G4int protonNumberCurrent

Definition: G4NucleiModel.hh:289

Here is the caller graph for this function:

◆ getNumberOfZones()

G4int G4NucleiModel::getNumberOfZones ( ) const

inline

Definition at line 143 of file G4NucleiModel.hh.

143 { return number_of_zones; }

G4NucleiModel::number_of_zones

G4int number_of_zones

Definition: G4NucleiModel.hh:279

◆ getPotential()

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

inline

Definition at line 122 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;
   }

Here is the caller graph for this function:

◆ getRadius() [1/2]

G4double G4NucleiModel::getRadius ( ) const

inline

Definition at line 133 of file G4NucleiModel.hh.

133 { return nuclei_radius; }

G4NucleiModel::nuclei_radius

G4double nuclei_radius

Definition: G4NucleiModel.hh:277

◆ getRadius() [2/2]

G4double G4NucleiModel::getRadius ( G4int izone ) const

inline

Definition at line 134 of file G4NucleiModel.hh.

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

◆ getRadiusUnits()

G4double G4NucleiModel::getRadiusUnits ( ) const

inline

Definition at line 131 of file G4NucleiModel.hh.

131 { return radiusUnits*CLHEP::fermi; }

CLHEP::fermi

static const double fermi

Definition: SystemOfUnits.h:82

G4NucleiModel::radiusUnits

const G4double radiusUnits

Definition: G4NucleiModel.hh:301

Here is the caller graph for this function:

◆ getRatio()

G4double G4NucleiModel::getRatio ( G4int ip ) const

protected

Definition at line 1343 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.;
 }

Here is the caller graph for this function:

◆ getTypesOfNucleonsInvolved()

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

inline

Definition at line 169 of file G4NucleiModel.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getVolume()

G4double G4NucleiModel::getVolume ( G4int izone ) const

inline

Definition at line 138 of file G4NucleiModel.hh.

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

Here is the caller graph for this function:

◆ getZone()

G4int G4NucleiModel::getZone ( G4double r ) const

inline

Definition at line 144 of file G4NucleiModel.hh.

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

Here is the caller graph for this function:

◆ initializeCascad() [1/2]

G4CascadParticle G4NucleiModel::initializeCascad ( G4InuclElementaryParticle * particle )

Definition at line 1388 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ initializeCascad() [2/2]

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

Definition at line 1413 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
 }

Here is the call graph for this function:

◆ inverseMeanFreePath()

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

protected

Definition at line 1814 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ isProjectile()

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

Definition at line 1308 of file G4NucleiModel.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ passFermi()

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

protected

Definition at line 1061 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; 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ passTrailing()

G4bool G4NucleiModel::passTrailing ( const G4ThreeVector & hit_position )

protected

Definition at line 1090 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
 }

Here is the caller graph for this function:

◆ printModel()

void G4NucleiModel::printModel ( ) const

Definition at line 604 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ reset()

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

Definition at line 290 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;
 }

Here is the caller graph for this function:

◆ setVerboseLevel()

void G4NucleiModel::setVerboseLevel ( G4int verbose )

inline

Definition at line 101 of file G4NucleiModel.hh.

101 { verboseLevel = verbose; }

G4NucleiModel::verboseLevel

G4int verboseLevel

Definition: G4NucleiModel.hh:247

Here is the call graph for this function:

Here is the caller graph for this function:

◆ sortPartners()

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

inlinestaticprotected

Definition at line 211 of file G4NucleiModel.hh.

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

Here is the call graph for this function:

Here is the caller graph for this function:

◆ stillInside()

G4bool G4NucleiModel::stillInside ( const G4CascadParticle & cparticle )

inline

Definition at line 156 of file G4NucleiModel.hh.

                                                         {
     return cparticle.getCurrentZone() < number_of_zones;
   }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ totalCrossSection()

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

Definition at line 1918 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));
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ useQuasiDeuteron()

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

static

Definition at line 1049 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
 }

Here is the caller graph for this function:

◆ worthToPropagate()

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

Definition at line 1312 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ zoneIntegralGaussian()

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

protected

Definition at line 556 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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ zoneIntegralWoodsSaxon()

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

protected

Definition at line 503 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))));
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

◆ A

G4int G4NucleiModel::A

private

Definition at line 281 of file G4NucleiModel.hh.

◆ acsecs

std::vector<G4double> G4NucleiModel::acsecs

private

Definition at line 254 of file G4NucleiModel.hh.

◆ alfa3

const G4double G4NucleiModel::alfa3 = { 0.7, 0.3, 0.01 }

staticprivate

Definition at line 316 of file G4NucleiModel.hh.

◆ alfa6

const G4double G4NucleiModel::alfa6 = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 }

staticprivate

Definition at line 316 of file G4NucleiModel.hh.

◆ binding_energies

std::vector<G4double> G4NucleiModel::binding_energies

private

Definition at line 276 of file G4NucleiModel.hh.

◆ collisionPts

std::vector<G4ThreeVector> G4NucleiModel::collisionPts

private

Definition at line 260 of file G4NucleiModel.hh.

◆ coordinates

std::vector<G4ThreeVector> G4NucleiModel::coordinates

private

Definition at line 256 of file G4NucleiModel.hh.

◆ crossSectionUnits

const G4double G4NucleiModel::crossSectionUnits

private

Definition at line 300 of file G4NucleiModel.hh.

◆ current_nucl1

G4int G4NucleiModel::current_nucl1

private

Definition at line 291 of file G4NucleiModel.hh.

◆ current_nucl2

G4int G4NucleiModel::current_nucl2

private

Definition at line 292 of file G4NucleiModel.hh.

◆ dummy_convertor

G4LorentzConvertor G4NucleiModel::dummy_convertor

private

Definition at line 250 of file G4NucleiModel.hh.

◆ EPCoutput

G4CollisionOutput G4NucleiModel::EPCoutput

private

Definition at line 251 of file G4NucleiModel.hh.

◆ fermi_momenta

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

private

Definition at line 273 of file G4NucleiModel.hh.

◆ fermiMomentum

const G4double G4NucleiModel::fermiMomentum

private

Definition at line 307 of file G4NucleiModel.hh.

◆ gammaQDinterp

G4CascadeInterpolator<30> G4NucleiModel::gammaQDinterp

private

Definition at line 294 of file G4NucleiModel.hh.

◆ gammaQDscale

const G4double G4NucleiModel::gammaQDscale

private

Definition at line 309 of file G4NucleiModel.hh.

◆ hyperon_vp

const G4double G4NucleiModel::hyperon_vp = 0.030

staticprivate

Definition at line 320 of file G4NucleiModel.hh.

◆ kaon_vp

const G4double G4NucleiModel::kaon_vp = 0.015

staticprivate

Definition at line 319 of file G4NucleiModel.hh.

◆ large

const G4double G4NucleiModel::large = 1000.

staticprivate

Definition at line 313 of file G4NucleiModel.hh.

◆ momentums

std::vector<G4LorentzVector> G4NucleiModel::momentums

private

Definition at line 257 of file G4NucleiModel.hh.

◆ neutronEP

const G4InuclElementaryParticle G4NucleiModel::neutronEP

private

Definition at line 323 of file G4NucleiModel.hh.

◆ neutronNumber

G4int G4NucleiModel::neutronNumber

private

Definition at line 285 of file G4NucleiModel.hh.

◆ neutronNumberCurrent

G4int G4NucleiModel::neutronNumberCurrent

private

Definition at line 288 of file G4NucleiModel.hh.

◆ nuclei_radius

G4double G4NucleiModel::nuclei_radius

private

Definition at line 277 of file G4NucleiModel.hh.

◆ nuclei_volume

G4double G4NucleiModel::nuclei_volume

private

Definition at line 278 of file G4NucleiModel.hh.

◆ nucleon_densities

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

private

Definition at line 271 of file G4NucleiModel.hh.

◆ number_of_zones

G4int G4NucleiModel::number_of_zones

private

Definition at line 279 of file G4NucleiModel.hh.

◆ pf

std::vector<G4double> G4NucleiModel::pf

private

Definition at line 267 of file G4NucleiModel.hh.

◆ pion_vp

const G4double G4NucleiModel::pion_vp = 0.007

staticprivate

Definition at line 317 of file G4NucleiModel.hh.

◆ pion_vp_small

const G4double G4NucleiModel::pion_vp_small = 0.007

staticprivate

Definition at line 318 of file G4NucleiModel.hh.

◆ piTimes4thirds

const G4double G4NucleiModel::piTimes4thirds = pi*4./3.

staticprivate

Definition at line 314 of file G4NucleiModel.hh.

◆ protonEP

const G4InuclElementaryParticle G4NucleiModel::protonEP

private

Definition at line 324 of file G4NucleiModel.hh.

◆ protonNumber

G4int G4NucleiModel::protonNumber

private

Definition at line 286 of file G4NucleiModel.hh.

◆ protonNumberCurrent

G4int G4NucleiModel::protonNumberCurrent

private

Definition at line 289 of file G4NucleiModel.hh.

◆ qdeutrons

std::vector<G4InuclElementaryParticle> G4NucleiModel::qdeutrons

private

Definition at line 253 of file G4NucleiModel.hh.

◆ R_nucleon

const G4double G4NucleiModel::R_nucleon

private

Definition at line 308 of file G4NucleiModel.hh.

◆ radiusForSmall

const G4double G4NucleiModel::radiusForSmall

private

Definition at line 305 of file G4NucleiModel.hh.

◆ radiusScale

const G4double G4NucleiModel::radiusScale

private

Definition at line 303 of file G4NucleiModel.hh.

◆ radiusScale2

const G4double G4NucleiModel::radiusScale2

private

Definition at line 304 of file G4NucleiModel.hh.

◆ radiusUnits

const G4double G4NucleiModel::radiusUnits

private

Definition at line 301 of file G4NucleiModel.hh.

◆ radScaleAlpha

const G4double G4NucleiModel::radScaleAlpha

private

Definition at line 306 of file G4NucleiModel.hh.

◆ raw_particles

std::vector<G4InuclElementaryParticle> G4NucleiModel::raw_particles

private

Definition at line 258 of file G4NucleiModel.hh.

◆ rod

std::vector<G4double> G4NucleiModel::rod

private

Definition at line 266 of file G4NucleiModel.hh.

◆ skinDepth

const G4double G4NucleiModel::skinDepth

private

Definition at line 302 of file G4NucleiModel.hh.

◆ small

const G4double G4NucleiModel::small = 1.0e-9

staticprivate

Definition at line 312 of file G4NucleiModel.hh.

◆ theNucleus

G4InuclNuclei* G4NucleiModel::theNucleus

private

Definition at line 283 of file G4NucleiModel.hh.

◆ thePartners

std::vector<partner> G4NucleiModel::thePartners

protected

Definition at line 207 of file G4NucleiModel.hh.

◆ ur

G4double G4NucleiModel::ur[7]

private

Definition at line 263 of file G4NucleiModel.hh.

◆ v

G4double G4NucleiModel::v[6]

private

Definition at line 264 of file G4NucleiModel.hh.

◆ v1

G4double G4NucleiModel::v1[6]

private

Definition at line 265 of file G4NucleiModel.hh.

◆ verboseLevel

G4int G4NucleiModel::verboseLevel

private

Definition at line 247 of file G4NucleiModel.hh.

◆ vz

std::vector<G4double> G4NucleiModel::vz

private

Definition at line 268 of file G4NucleiModel.hh.

◆ Z

G4int G4NucleiModel::Z

private

Definition at line 282 of file G4NucleiModel.hh.

◆ zone_potentials

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

private

Definition at line 272 of file G4NucleiModel.hh.

◆ zone_radii

std::vector<G4double> G4NucleiModel::zone_radii

private

Definition at line 274 of file G4NucleiModel.hh.

◆ zone_volumes

std::vector<G4double> G4NucleiModel::zone_volumes

private

Definition at line 275 of file G4NucleiModel.hh.

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

Public Types

Public Member Functions

Static Public Member Functions

Protected Types

Protected Member Functions

Static Protected Member Functions

Protected Attributes

Private Types

Private Attributes

Static Private Attributes

Detailed Description

Member Typedef Documentation

◆ modelLists

◆ partner

Member Enumeration Documentation

◆ PotentialType

Constructor & Destructor Documentation

◆ G4NucleiModel() [1/3]

◆ G4NucleiModel() [2/3]

◆ G4NucleiModel() [3/3]

◆ ~G4NucleiModel()

Member Function Documentation

◆ absorptionCrossSection()

◆ boundaryTransition()

◆ choosePointAlongTraj()

◆ empty()

◆ fillBindingEnergies()

◆ fillPotentials()

◆ fillZoneRadii()

◆ fillZoneVolumes()

◆ forceFirst()

◆ generateInteractionLength()

◆ generateInteractionPartners()

◆ generateModel() [1/2]

◆ generateModel() [2/2]

◆ generateNucleon()

◆ generateNucleonMomentum()

◆ generateParticleFate()

◆ generateQuasiDeuteron()

◆ getCurrentDensity()

◆ getDensity()

◆ getFermiKinetic()

◆ getFermiMomentum()

◆ getNumberOfNeutrons()

◆ getNumberOfProtons()

◆ getNumberOfZones()

◆ getPotential()

◆ getRadius() [1/2]

◆ getRadius() [2/2]

◆ getRadiusUnits()

◆ getRatio()

◆ getTypesOfNucleonsInvolved()

◆ getVolume()

◆ getZone()

◆ initializeCascad() [1/2]

◆ initializeCascad() [2/2]

◆ inverseMeanFreePath()

◆ isProjectile()

◆ passFermi()

◆ passTrailing()

◆ printModel()

◆ reset()

◆ setVerboseLevel()

◆ sortPartners()

◆ stillInside()

◆ totalCrossSection()

◆ useQuasiDeuteron()

◆ worthToPropagate()

◆ zoneIntegralGaussian()

◆ zoneIntegralWoodsSaxon()

Member Data Documentation

◆ A

◆ acsecs

◆ alfa3

◆ alfa6

◆ binding_energies

◆ collisionPts

◆ coordinates

◆ crossSectionUnits

◆ current_nucl1