G4NucleiModel Class Reference

#include <G4NucleiModel.hh>

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

Detailed Description

Definition at line 93 of file G4NucleiModel.hh.

Member Typedef Documentation

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

Definition at line 163 of file G4NucleiModel.hh.

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

protected

Definition at line 205 of file G4NucleiModel.hh.

Constructor & Destructor Documentation

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::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::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 ( )

virtual

Definition at line 282 of file G4NucleiModel.cc.

                              {
  delete theNucleus;
  theNucleus = 0;
}

Member Function Documentation

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:

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:

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:

G4bool G4NucleiModel::empty ( ) const

inline

Definition at line 152 of file G4NucleiModel.hh.

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

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 call graph for this function:

Here is the caller graph for this function:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

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 caller graph for this function:

G4int G4NucleiModel::getNumberOfNeutrons ( ) const

inline

Definition at line 149 of file G4NucleiModel.hh.

{ return neutronNumberCurrent; }

G4int G4NucleiModel::getNumberOfProtons ( ) const

inline

Definition at line 150 of file G4NucleiModel.hh.

{ return protonNumberCurrent; }

G4int G4NucleiModel::getNumberOfZones ( ) const

inline

Definition at line 143 of file G4NucleiModel.hh.

{ return number_of_zones; }

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:

G4double G4NucleiModel::getRadius ( ) const

inline

Definition at line 133 of file G4NucleiModel.hh.

{ return nuclei_radius; }

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);
  }

G4double G4NucleiModel::getRadiusUnits ( ) const

inline

Definition at line 131 of file G4NucleiModel.hh.

{ return radiusUnits*CLHEP::fermi; }

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:

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);
  }

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:

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:

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:

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:

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:

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:

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:

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:

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:

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:

void G4NucleiModel::setVerboseLevel ( G4int  verbose )

inline

Definition at line 101 of file G4NucleiModel.hh.

{ verboseLevel = verbose; }

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 caller graph for this function:

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:

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:

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:

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:

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:

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

std::vector<partner> G4NucleiModel::thePartners

protected

Definition at line 207 of file G4NucleiModel.hh.

---

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

• geant4.10.03.p01/source/processes/hadronic/models/cascade/cascade/include/G4NucleiModel.hh
• geant4.10.03.p01/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc