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 91 of file G4NucleiModel.hh.

Member Typedef Documentation

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

Definition at line 161 of file G4NucleiModel.hh.

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

protected

Definition at line 203 of file G4NucleiModel.hh.

Constructor & Destructor Documentation

G4NucleiModel::G4NucleiModel

(

)

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

G4NucleiModel::G4NucleiModel	(	G4int	a,
		G4int	z
	)

Definition at line 259 of file G4NucleiModel.cc.

References generateModel().

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

G4NucleiModel::G4NucleiModel ( G4InuclNuclei *  nuclei )

explicit

Definition at line 268 of file G4NucleiModel.cc.

References generateModel().

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

G4NucleiModel::~G4NucleiModel ( )

virtual

Definition at line 277 of file G4NucleiModel.cc.

                              {
  delete theNucleus;
  theNucleus = 0;
}

Member Function Documentation

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

Definition at line 1872 of file G4NucleiModel.cc.

References G4cerr, G4cout, G4endl, G4CascadeInterpolator< NBINS >::interpolate(), G4InuclParticleNames::muonMinus, photon, G4InuclParticleNames::pionMinus, G4InuclParticleNames::pionPlus, G4InuclParticleNames::pionZero, and useQuasiDeuteron().

Referenced by inverseMeanFreePath().

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

void G4NucleiModel::boundaryTransition ( G4CascadParticle &  cparticle )

protected

Definition at line 1097 of file G4NucleiModel.cc.

References CLHEP::Hep3Vector::dot(), CLHEP::HepLorentzVector::e(), G4cerr, G4cout, G4endl, G4CascadParticle::getCurrentZone(), G4CascadParticle::getMomentum(), G4CascadParticle::getParticle(), G4CascadParticle::getPosition(), getPotential(), G4CascadParticle::incrementReflectionCounter(), CLHEP::Hep3Vector::mag(), G4CascadParticle::movingInsideNuclei(), G4InuclParticleNames::pos, gammaraytel::pr, G4CascadParticle::resetReflection(), CLHEP::HepLorentzVector::setVect(), G4InuclElementaryParticle::type(), G4CascadParticle::updateParticleMomentum(), G4CascadParticle::updateZone(), CLHEP::HepLorentzVector::vect(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by generateParticleFate().

                                                                  {
  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;
  
  G4double dv = getPotential(type,zone) - getPotential(type, next_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;
  
  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);
}

void G4NucleiModel::choosePointAlongTraj ( G4CascadParticle &  cparticle )

protected

Definition at line 1163 of file G4NucleiModel.cc.

References CLHEP::Hep3Vector::angle(), CLHEP::Hep3Vector::cross(), python.hepunit::deg, G4cout, G4endl, G4Exp(), G4UniformRand, G4CascadParticle::getCurrentZone(), G4CascadParticle::getMomentum(), G4CascadParticle::getPosition(), getZone(), inverseMeanFreePath(), iz, CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), python.hepunit::pi, G4InuclParticleNames::pos, CLHEP::Hep3Vector::rotate(), CLHEP::Hep3Vector::set(), CLHEP::Hep3Vector::unit(), G4CascadParticle::updatePosition(), G4CascadParticle::updateZone(), and CLHEP::HepLorentzVector::vect().

Referenced by initializeCascad().

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

G4bool G4NucleiModel::empty ( ) const

inline

Definition at line 150 of file G4NucleiModel.hh.

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

void G4NucleiModel::fillBindingEnergies ( )

protected

Definition at line 377 of file G4NucleiModel.cc.

References G4InuclSpecialFunctions::bindingEnergy(), G4cout, G4endl, and python.hepunit::GeV.

Referenced by generateModel().

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

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

protected

Definition at line 466 of file G4NucleiModel.cc.

References G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, G4InuclElementaryParticle::getParticleMass(), neutron, and G4InuclParticleNames::proton.

Referenced by generateModel().

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

void G4NucleiModel::fillZoneRadii ( G4double  nuclearRadius )

protected

Definition at line 391 of file G4NucleiModel.cc.

References G4cout, G4endl, G4Exp(), and G4Log().

Referenced by generateModel().

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

G4double G4NucleiModel::fillZoneVolumes ( G4double  nuclearRadius )

protected

Definition at line 431 of file G4NucleiModel.cc.

References G4cout, G4endl, zoneIntegralGaussian(), and zoneIntegralWoodsSaxon().

Referenced by generateModel().

                                                              {
  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
}

G4bool G4NucleiModel::forceFirst

(

const G4CascadParticle &

cparticle

)

const

Definition at line 1289 of file G4NucleiModel.cc.

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

Referenced by generateInteractionLength(), and initializeCascad().

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

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

protected

Definition at line 1840 of file G4NucleiModel.cc.

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

Referenced by generateInteractionPartners().

                                                           {
  // 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;
}

void G4NucleiModel::generateInteractionPartners ( G4CascadParticle &  cparticle )

protected

Definition at line 681 of file G4NucleiModel.cc.

References G4cerr, G4cout, G4endl, generateInteractionLength(), generateNucleon(), generateQuasiDeuteron(), G4CascadParticle::getCurrentZone(), G4CascadParticle::getMomentum(), G4CascadParticle::getParticle(), G4CascadParticle::getPathToTheNextZone(), G4InuclSpecialFunctions::inuclRndm(), inverseMeanFreePath(), isProjectile(), CLHEP::Hep3Vector::mag(), G4InuclParticleNames::mum, G4InuclParticleNames::muonMinus, G4InuclParticleNames::neu, neutron, G4InuclParticleNames::pim, G4InuclParticleNames::pip, G4InuclParticleNames::pro, G4InuclParticleNames::proton, sort(), sortPartners(), thePartners, G4InuclElementaryParticle::type(), useQuasiDeuteron(), and CLHEP::HepLorentzVector::vect().

Referenced by generateParticleFate().

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

void G4NucleiModel::generateModel

(

G4InuclNuclei *

nuclei

)

Definition at line 298 of file G4NucleiModel.cc.

References G4InuclNuclei::getA(), and G4InuclNuclei::getZ().

Referenced by G4NucleiModel().

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

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

Definition at line 302 of file G4NucleiModel.cc.

References test::a, fillBindingEnergies(), fillPotentials(), fillZoneRadii(), fillZoneVolumes(), G4InuclSpecialFunctions::G4cbrt(), G4cout, G4endl, neutron, printModel(), G4InuclParticleNames::proton, reset(), and z.

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

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

Definition at line 644 of file G4NucleiModel.cc.

References G4cout, G4endl, and generateNucleonMomentum().

Referenced by generateInteractionPartners().

                                                           {
  if (verboseLevel > 1) {
    G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
  }

  G4LorentzVector mom = generateNucleonMomentum(type, zone);
  return G4InuclElementaryParticle(mom, type);
}

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

Definition at line 635 of file G4NucleiModel.cc.

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

Referenced by generateNucleon(), and generateQuasiDeuteron().

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

  return generateWithRandomAngles(pmod, mass);
}

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

Definition at line 852 of file G4NucleiModel.cc.

References boundaryTransition(), G4ElementaryParticleCollider::collide(), G4CascadeCheckBalance::collide(), G4cerr, G4cout, G4endl, G4ThreadLocal, generateInteractionPartners(), G4InuclParticle::getCharge(), G4CascadParticle::getCurrentZone(), G4CascadParticle::getGeneration(), G4CollisionOutput::getOutgoingParticles(), G4CascadParticle::getParticle(), G4CascadParticle::getPosition(), G4CascadParticle::incrementCurrentPath(), G4InuclElementaryParticle::nucleon(), G4CollisionOutput::numberOfOutgoingParticles(), G4CascadeCheckBalance::okay(), passFermi(), passTrailing(), G4CollisionOutput::printCollisionOutput(), G4CascadParticle::propagateAlongThePath(), G4InuclElementaryParticle::quasi_deutron(), G4CollisionOutput::reset(), python.hepunit::second, G4VCascadeCollider::setVerboseLevel(), sort(), thePartners, G4InuclElementaryParticle::type(), and G4CascadParticle::updatePosition().

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

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

protected

Definition at line 655 of file G4NucleiModel.cc.

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

Referenced by generateInteractionPartners().

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

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

protected

Definition at line 1348 of file G4NucleiModel.cc.

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

Referenced by inverseMeanFreePath().

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

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

inline

Definition at line 110 of file G4NucleiModel.hh.

Referenced by getCurrentDensity(), and printModel().

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

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

Definition at line 621 of file G4NucleiModel.cc.

References G4InuclElementaryParticle::getParticleMass().

Referenced by worthToPropagate().

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

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

inline

Definition at line 114 of file G4NucleiModel.hh.

Referenced by generateNucleonMomentum(), and printModel().

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

G4int G4NucleiModel::getNumberOfNeutrons ( ) const

inline

Definition at line 147 of file G4NucleiModel.hh.

{ return neutronNumberCurrent; }

G4int G4NucleiModel::getNumberOfProtons ( ) const

inline

Definition at line 148 of file G4NucleiModel.hh.

{ return protonNumberCurrent; }

G4int G4NucleiModel::getNumberOfZones ( ) const

inline

Definition at line 141 of file G4NucleiModel.hh.

{ return number_of_zones; }

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

inline

Definition at line 120 of file G4NucleiModel.hh.

Referenced by boundaryTransition(), and printModel().

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

G4double G4NucleiModel::getRadius ( ) const

inline

Definition at line 131 of file G4NucleiModel.hh.

{ return nuclei_radius; }

G4double G4NucleiModel::getRadius ( G4int  izone ) const

inline

Definition at line 132 of file G4NucleiModel.hh.

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

G4double G4NucleiModel::getRadiusUnits ( ) const

inline

Definition at line 129 of file G4NucleiModel.hh.

{ return radiusUnits*CLHEP::fermi; }

G4double G4NucleiModel::getRatio ( G4int  ip ) const

protected

Definition at line 1331 of file G4NucleiModel.cc.

References G4InuclParticleNames::dineutron, G4InuclParticleNames::diproton, G4cout, G4endl, neutron, G4InuclParticleNames::proton, and G4InuclParticleNames::unboundPN.

Referenced by getCurrentDensity().

                                               {
  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.;
}

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

inline

Definition at line 167 of file G4NucleiModel.hh.

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

G4double G4NucleiModel::getVolume ( G4int  izone ) const

inline

Definition at line 136 of file G4NucleiModel.hh.

Referenced by getCurrentDensity().

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

G4int G4NucleiModel::getZone ( G4double  r ) const

inline

Definition at line 142 of file G4NucleiModel.hh.

References iz.

Referenced by choosePointAlongTraj().

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

G4CascadParticle G4NucleiModel::initializeCascad

(

G4InuclElementaryParticle *

particle

)

Definition at line 1376 of file G4NucleiModel.cc.

References choosePointAlongTraj(), forceFirst(), G4cout, G4endl, G4InuclSpecialFunctions::generateWithFixedTheta(), G4InuclParticle::getKineticEnergy(), G4InuclSpecialFunctions::inuclRndm(), and G4InuclParticleNames::pos.

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

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

Definition at line 1401 of file G4NucleiModel.cc.

References test::b, G4LorentzConvertor::backToTheLab(), G4InuclSpecialFunctions::bindingEnergy(), CLHEP::Hep3Vector::dot(), G4cout, G4endl, G4Exp(), G4Log(), G4InuclSpecialFunctions::generateWithRandomAngles(), G4InuclNuclei::getA(), G4InuclParticle::getEnergy(), G4InuclParticle::getKineticEnergy(), G4InuclParticle::getMass(), G4InuclNuclei::getZ(), python.hepunit::GeV, G4InuclSpecialFunctions::inuclRndm(), CLHEP::Hep3Vector::mag(), G4InuclSpecialFunctions::randomPHI(), CLHEP::HepLorentzVector::rho(), CLHEP::Hep3Vector::set(), CLHEP::HepLorentzVector::setE(), CLHEP::HepLorentzVector::setVect(), plottest35::t1, G4LorentzConvertor::toTheTargetRestFrame(), CLHEP::HepLorentzVector::vect(), test::x, z, and CLHEP::HepLorentzVector::z().

                                         {
  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;

      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) {
        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) {
          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) {
          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) {

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

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

protected

Definition at line 1801 of file G4NucleiModel.cc.

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

Referenced by choosePointAlongTraj(), and generateInteractionPartners().

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

G4bool G4NucleiModel::isProjectile

(

const G4CascadParticle &

cparticle

)

const

Definition at line 1296 of file G4NucleiModel.cc.

References G4CascadParticle::getGeneration().

Referenced by forceFirst(), and generateInteractionPartners().

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

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

protected

Definition at line 1051 of file G4NucleiModel.cc.

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

Referenced by generateParticleFate().

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

G4bool G4NucleiModel::passTrailing ( const G4ThreeVector &  hit_position )

protected

Definition at line 1080 of file G4NucleiModel.cc.

References G4cout, and G4endl.

Referenced by generateParticleFate().

                                                                    {
  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
}

void G4NucleiModel::printModel

(

)

const

Definition at line 599 of file G4NucleiModel.cc.

References G4cout, G4endl, getDensity(), getFermiMomentum(), and getPotential().

Referenced by generateModel().

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

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

Definition at line 285 of file G4NucleiModel.cc.

Referenced by generateModel().

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

void G4NucleiModel::setVerboseLevel ( G4int  verbose )

inline

Definition at line 99 of file G4NucleiModel.hh.

{ verboseLevel = verbose; }

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

inlinestaticprotected

Definition at line 209 of file G4NucleiModel.hh.

Referenced by generateInteractionPartners().

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

G4bool G4NucleiModel::stillInside ( const G4CascadParticle &  cparticle )

inline

Definition at line 154 of file G4NucleiModel.hh.

References G4CascadParticle::getCurrentZone().

                                                        {
    return cparticle.getCurrentZone() < number_of_zones;
  }

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

Definition at line 1905 of file G4NucleiModel.cc.

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

Referenced by inverseMeanFreePath().

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

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

static

Definition at line 1039 of file G4NucleiModel.cc.

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

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

                                                                {
  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
}

G4bool G4NucleiModel::worthToPropagate

(

const G4CascadParticle &

cparticle

)

const

Definition at line 1300 of file G4NucleiModel.cc.

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

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

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

protected

Definition at line 551 of file G4NucleiModel.cc.

References G4cerr, G4cout, G4endl, and G4Exp().

Referenced by fillZoneVolumes().

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

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

protected

Definition at line 498 of file G4NucleiModel.cc.

References G4cout, G4endl, G4Exp(), and G4Log().

Referenced by fillZoneVolumes().

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

Field Documentation

std::vector<partner> G4NucleiModel::thePartners

protected

Definition at line 205 of file G4NucleiModel.hh.

Referenced by generateInteractionPartners(), and generateParticleFate().

---

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

• G4NucleiModel.hh
• G4NucleiModel.cc