Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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, G4intgetTypesOfNucleonsInvolved () 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< partnerthePartners
 

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.

Definition at line 205 of file G4NucleiModel.hh.

Constructor & Destructor Documentation

G4NucleiModel::G4NucleiModel ( )

Definition at line 227 of file G4NucleiModel.cc.

228  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
229  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
230  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
231  current_nucl2(0), gammaQDinterp(kebins),
232  crossSectionUnits(G4CascadeParameters::xsecScale()),
233  radiusUnits(G4CascadeParameters::radiusScale()),
234  skinDepth(0.611207*radiusUnits),
235  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
236  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
237  radiusForSmall(G4CascadeParameters::radiusSmall()),
238  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
239  fermiMomentum(G4CascadeParameters::fermiScale()),
241  gammaQDscale(G4CascadeParameters::gammaQDScale()),
242  neutronEP(neutron), protonEP(proton) {}
static G4double xsecScale()
static G4double radiusTrailing()
static G4double radiusScale()
static G4double fermiScale()
static G4double radiusAlpha()
static G4bool useTwoParam()
static G4double radiusSmall()
static G4double gammaQDScale()
G4NucleiModel::G4NucleiModel ( G4int  a,
G4int  z 
)

Definition at line 244 of file G4NucleiModel.cc.

245  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
246  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
247  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
248  current_nucl2(0), gammaQDinterp(kebins),
249  crossSectionUnits(G4CascadeParameters::xsecScale()),
250  radiusUnits(G4CascadeParameters::radiusScale()),
251  skinDepth(0.611207*radiusUnits),
252  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
253  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
254  radiusForSmall(G4CascadeParameters::radiusSmall()),
255  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
256  fermiMomentum(G4CascadeParameters::fermiScale()),
258  gammaQDscale(G4CascadeParameters::gammaQDScale()),
259  neutronEP(neutron), protonEP(proton) {
260  generateModel(a,z);
261 }
static G4double xsecScale()
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static G4double radiusTrailing()
static G4double radiusScale()
void generateModel(G4InuclNuclei *nuclei)
static G4double fermiScale()
static G4double radiusAlpha()
static G4bool useTwoParam()
tuple z
Definition: test.py:28
static G4double radiusSmall()
static G4double gammaQDScale()

Here is the call graph for this function:

G4NucleiModel::G4NucleiModel ( G4InuclNuclei nuclei)
explicit

Definition at line 263 of file G4NucleiModel.cc.

264  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
265  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
266  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
267  current_nucl2(0), gammaQDinterp(kebins),
268  crossSectionUnits(G4CascadeParameters::xsecScale()),
269  radiusUnits(G4CascadeParameters::radiusScale()),
270  skinDepth(0.611207*radiusUnits),
271  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
272  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
273  radiusForSmall(G4CascadeParameters::radiusSmall()),
274  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
275  fermiMomentum(G4CascadeParameters::fermiScale()),
277  gammaQDscale(G4CascadeParameters::gammaQDScale()),
278  neutronEP(neutron), protonEP(proton) {
279  generateModel(nuclei);
280 }
static G4double xsecScale()
static G4double radiusTrailing()
static G4double radiusScale()
void generateModel(G4InuclNuclei *nuclei)
static G4double fermiScale()
static G4double radiusAlpha()
static G4bool useTwoParam()
static G4double radiusSmall()
static G4double gammaQDScale()

Here is the call graph for this function:

G4NucleiModel::~G4NucleiModel ( )
virtual

Definition at line 282 of file G4NucleiModel.cc.

282  {
283  delete theNucleus;
284  theNucleus = 0;
285 }

Member Function Documentation

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

Definition at line 1885 of file G4NucleiModel.cc.

1885  {
1886  if (!useQuasiDeuteron(type)) {
1887  G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
1888  return 0.;
1889  }
1890 
1891  G4double csec = 0.;
1892 
1893  // Pion absorption is parametrized for low vs. medium energy
1894  // ... use for muon capture as well
1895  if (type == pionPlus || type == pionMinus || type == pionZero ||
1896  type == muonMinus) {
1897  if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
1898  + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
1899  else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);
1900  }
1901 
1902  // Photon cross-section is binned from parametrization by Mi. Kossov
1903  // See below for initialization of gammaQDinterp, gammaQDxsec
1904  if (type == photon) {
1905  csec = gammaQDinterp.interpolate(ke, gammaQDxsec) * gammaQDscale;
1906  }
1907 
1908  if (csec < 0.0) csec = 0.0;
1909 
1910  if (verboseLevel > 2) {
1911  G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;
1912  }
1913 
1914  return crossSectionUnits * csec;
1915 }
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
G4GLOB_DLL std::ostream G4cout
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr

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.

1107  {
1108  if (verboseLevel > 1) {
1109  G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
1110  }
1111 
1112  G4int zone = cparticle.getCurrentZone();
1113 
1114  if (cparticle.movingInsideNuclei() && zone == 0) {
1115  if (verboseLevel) G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
1116  return;
1117  }
1118 
1119  G4LorentzVector mom = cparticle.getMomentum();
1120  G4ThreeVector pos = cparticle.getPosition();
1121 
1122  G4int type = cparticle.getParticle().type();
1123 
1124  G4double r = pos.mag();
1125  G4double pr = pos.dot(mom.vect()) / r;
1126 
1127  G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
1128 
1129  // NOTE: dv is the height of the wall seen by the particle
1130  G4double dv = getPotential(type,next_zone) - getPotential(type,zone);
1131  if (verboseLevel > 3) {
1132  G4cout << "Potentials for type " << type << " = "
1133  << getPotential(type,zone) << " , "
1134  << getPotential(type,next_zone) << G4endl;
1135  }
1136 
1137  G4double qv = dv * dv + 2.0 * dv * mom.e() + pr * pr;
1138 
1139  G4double p1r = 0.;
1140 
1141  if (verboseLevel > 3) {
1142  G4cout << " type " << type << " zone " << zone << " next " << next_zone
1143  << " qv " << qv << " dv " << dv << G4endl;
1144  }
1145 
1146  if (qv <= 0.0) { // reflection
1147  if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
1148  p1r = -pr;
1149  cparticle.incrementReflectionCounter();
1150  } else { // transition
1151  if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
1152  p1r = std::sqrt(qv);
1153  if (pr < 0.0) p1r = -p1r;
1154 
1155  cparticle.updateZone(next_zone);
1156  cparticle.resetReflection();
1157  }
1158 
1159  G4double prr = (p1r - pr) / r;
1160 
1161  if (verboseLevel > 3) {
1162  G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
1163  << prr*pos.y() << " pz " << prr*pos.z() << " mag "
1164  << std::fabs(prr*r) << G4endl;
1165  }
1166 
1167  mom.setVect(mom.vect() + pos*prr);
1168  cparticle.updateParticleMomentum(mom);
1169 }
void updateParticleMomentum(const G4LorentzVector &mom)
G4LorentzVector getMomentum() const
void incrementReflectionCounter()
G4bool movingInsideNuclei() const
double x() const
double dot(const Hep3Vector &) const
int G4int
Definition: G4Types.hh:78
const G4InuclElementaryParticle & getParticle() const
double z() const
void updateZone(G4int izone)
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
double y() const
G4int getCurrentZone() const
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & getPosition() const
void setVect(const Hep3Vector &)
double G4double
Definition: G4Types.hh:76
double mag() const
static const G4double pos
G4GLOB_DLL std::ostream G4cerr
G4double getPotential(G4int ip, G4int izone) const

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.

1175  {
1176  if (verboseLevel > 1)
1177  G4cout << " >>> G4NucleiModel::choosePointAlongTraj" << G4endl;
1178 
1179  // Get trajectory through nucleus by computing exit point of line,
1180  // assuming that current position is on surface
1181 
1182  // FIXME: These need to be reusable (static) buffers
1183  G4ThreeVector pos = cparticle.getPosition();
1184  G4ThreeVector rhat = pos.unit();
1185 
1186  G4ThreeVector phat = cparticle.getMomentum().vect().unit();
1187  if (cparticle.getMomentum().vect().mag() < small) phat.set(0.,0.,1.);
1188 
1189  if (verboseLevel > 3)
1190  G4cout << " pos " << pos << " phat " << phat << " rhat " << rhat << G4endl;
1191 
1192  G4ThreeVector posout = pos;
1193  G4double prang = rhat.angle(-phat);
1194 
1195  if (prang < 1e-6) posout = -pos; // Check for radial incidence
1196  else {
1197  G4double posrot = 2.*prang - pi;
1198  posout.rotate(posrot, phat.cross(rhat));
1199  if (verboseLevel > 3) G4cout << " posrot " << posrot/deg << " deg";
1200  }
1201 
1202  if (verboseLevel > 3) G4cout << " posout " << posout << G4endl;
1203 
1204  // Get list of zone crossings along trajectory
1205  G4ThreeVector posmid = (pos+posout)/2.; // Midpoint of trajectory
1206  G4double r2mid = posmid.mag2();
1207  G4double lenmid = (posout-pos).mag()/2.; // Half-length of trajectory
1208 
1209  G4int zoneout = number_of_zones-1;
1210  G4int zonemid = getZone(posmid.mag()); // Middle zone traversed
1211 
1212  // Every zone is entered then exited, so preallocate vector
1213  G4int ncross = (number_of_zones-zonemid)*2;
1214 
1215  if (verboseLevel > 3) {
1216  G4cout << " posmid " << posmid << " lenmid " << lenmid
1217  << " zoneout " << zoneout << " zonemid " << zonemid
1218  << " ncross " << ncross << G4endl;
1219  }
1220 
1221  // FIXME: These need to be reusable (static) buffers
1222  std::vector<G4double> wtlen(ncross,0.); // CDF from entry point
1223  std::vector<G4double> len(ncross,0.); // Distance from entry point
1224 
1225  // Work from outside in, to accumulate CDF steps properly
1226  G4int i; // Loop variable, used multiple times
1227  for (i=0; i<ncross/2; i++) {
1228  G4int iz = zoneout-i;
1229  G4double ds = std::sqrt(zone_radii[iz]*zone_radii[iz]-r2mid);
1230 
1231  len[i] = lenmid - ds; // Distance from entry to crossing
1232  len[ncross-1-i] = lenmid + ds; // Distance to outbound crossing
1233 
1234  if (verboseLevel > 3) {
1235  G4cout << " i " << i << " iz " << iz << " ds " << ds
1236  << " len " << len[i] << G4endl;
1237  }
1238  }
1239 
1240  // Compute weights for each zone along trajectory and integrate
1241  for (i=1; i<ncross; i++) {
1242  G4int iz = (i<ncross/2) ? zoneout-i+1 : zoneout-ncross+i+1;
1243 
1244  G4double dlen = len[i]-len[i-1]; // Distance from previous crossing
1245 
1246  // Uniform probability across entire zone
1247  //*** G4double wt = dlen*(getDensity(1,iz)+getDensity(2,iz));
1248 
1249  // Probability based on interaction length through zone
1250  G4double invmfp = (inverseMeanFreePath(cparticle, neutronEP, iz)
1251  + inverseMeanFreePath(cparticle, protonEP, iz));
1252 
1253  // Integral of exp(-len/mfp) from start of zone to end
1254  G4double wt = (G4Exp(-len[i-1]*invmfp)-G4Exp(-len[i]*invmfp)) / invmfp;
1255 
1256  wtlen[i] = wtlen[i-1] + wt;
1257 
1258  if (verboseLevel > 3) {
1259  G4cout << " i " << i << " iz " << iz << " avg.mfp " << 1./invmfp
1260  << " dlen " << dlen << " wt " << wt << " wtlen " << wtlen[i]
1261  << G4endl;
1262  }
1263  }
1264 
1265  // Normalize CDF to unit integral
1266  std::transform(wtlen.begin(), wtlen.end(), wtlen.begin(),
1267  std::bind2nd(std::divides<G4double>(), wtlen.back()));
1268 
1269  if (verboseLevel > 3) {
1270  G4cout << " weights";
1271  for (i=0; i<ncross; i++) G4cout << " " << wtlen[i];
1272  G4cout << G4endl;
1273  }
1274 
1275  // Choose random point along trajectory, weighted by density
1276  G4double rand = G4UniformRand();
1277  G4int ir = std::upper_bound(wtlen.begin(),wtlen.end(),rand) - wtlen.begin();
1278 
1279  G4double frac = (rand-wtlen[ir-1]) / (wtlen[ir]-wtlen[ir-1]);
1280  G4double drand = (1.-frac)*len[ir-1] + frac*len[ir];
1281 
1282  if (verboseLevel > 3) {
1283  G4cout << " rand " << rand << " ir " << ir << " frac " << frac
1284  << " drand " << drand << G4endl;
1285  }
1286 
1287  pos += drand * phat; // Distance from entry point along trajectory
1288 
1289  // Update input particle with new location
1290  cparticle.updatePosition(pos);
1291  cparticle.updateZone(getZone(pos.mag()));
1292 
1293  if (verboseLevel > 2) {
1294  G4cout << " moved particle to zone " << cparticle.getCurrentZone()
1295  << " @ " << pos << G4endl;
1296  }
1297 }
const XML_Char int len
Definition: expat.h:262
void set(double x, double y, double z)
G4LorentzVector getMomentum() const
double angle(const Hep3Vector &) const
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
void updatePosition(const G4ThreeVector &pos)
int G4int
Definition: G4Types.hh:78
void updateZone(G4int izone)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
Hep3Vector unit() const
G4int getCurrentZone() const
double mag2() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
Hep3Vector cross(const Hep3Vector &) const
const G4ThreeVector & getPosition() const
double G4double
Definition: G4Types.hh:76
static constexpr double deg
Definition: G4SIunits.hh:152
double mag() const
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4int getZone(G4double r) const
static const G4double pos

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.

152  {
153  return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
154  }
void G4NucleiModel::fillBindingEnergies ( )
protected

Definition at line 382 of file G4NucleiModel.cc.

382  {
383  if (verboseLevel > 1)
384  G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
385 
386  G4double dm = bindingEnergy(A,Z);
387 
388  // Binding energy differences for proton and neutron loss, respectively
389  // FIXME: Why is fabs() used here instead of the signed difference?
390  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
391  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
392 }
G4GLOB_DLL std::ostream G4cout
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double bindingEnergy(G4int A, G4int Z)

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.

471  {
472  if (verboseLevel > 1)
473  G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
474 
475  if (type != proton && type != neutron) return;
476 
478 
479  // FIXME: This is the fabs() binding energy difference, not signed
480  const G4double dm = binding_energies[type-1];
481 
482  rod.clear(); rod.reserve(number_of_zones);
483  pf.clear(); pf.reserve(number_of_zones);
484  vz.clear(); vz.reserve(number_of_zones);
485 
486  G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
487  G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
488 
489  for (G4int i = 0; i < number_of_zones; i++) {
490  G4double rd = dd0 * v[i] / v1[i];
491  rod.push_back(rd);
492  G4double pff = fermiMomentum * G4cbrt(rd);
493  pf.push_back(pff);
494  vz.push_back(0.5 * pff * pff / mass + dm);
495  }
496 
497  nucleon_densities.push_back(rod);
498  fermi_momenta.push_back(pf);
499  zone_potentials.push_back(vz);
500 }
static G4double getParticleMass(G4int type)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

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.

396  {
397  if (verboseLevel > 1)
398  G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
399 
400  G4double skinRatio = nuclearRadius/skinDepth;
401  G4double skinDecay = G4Exp(-skinRatio);
402 
403  if (A < 5) { // Light ions treated as simple balls
404  zone_radii.push_back(nuclearRadius);
405  ur[0] = 0.;
406  ur[1] = 1.;
407  } else if (A < 12) { // Small nuclei have Gaussian potential
408  G4double rSq = nuclearRadius * nuclearRadius;
409  G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
410 
411  ur[0] = 0.0;
412  for (G4int i = 0; i < number_of_zones; i++) {
413  G4double y = std::sqrt(-G4Log(alfa3[i]));
414  zone_radii.push_back(gaussRadius * y);
415  ur[i+1] = y;
416  }
417  } else if (A < 100) { // Intermediate nuclei have three-zone W-S
418  ur[0] = -skinRatio;
419  for (G4int i = 0; i < number_of_zones; i++) {
420  G4double y = G4Log((1.0 + skinDecay)/alfa3[i] - 1.0);
421  zone_radii.push_back(nuclearRadius + skinDepth * y);
422  ur[i+1] = y;
423  }
424  } else { // Heavy nuclei have six-zone W-S
425  ur[0] = -skinRatio;
426  for (G4int i = 0; i < number_of_zones; i++) {
427  G4double y = G4Log((1.0 + skinDecay)/alfa6[i] - 1.0);
428  zone_radii.push_back(nuclearRadius + skinDepth * y);
429  ur[i+1] = y;
430  }
431  }
432 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

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.

436  {
437  if (verboseLevel > 1)
438  G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
439 
440  G4double tot_vol = 0.; // Return value omitting 4pi/3 for sphere!
441 
442  if (A < 5) { // Light ions are treated as simple balls
443  v[0] = v1[0] = 1.;
444  tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
445  zone_volumes.push_back(tot_vol*piTimes4thirds); // True volume of sphere
446 
447  return tot_vol;
448  }
449 
450  PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
451 
452  for (G4int i = 0; i < number_of_zones; i++) {
453  if (usePotential == WoodsSaxon) {
454  v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
455  } else {
456  v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
457  }
458 
459  tot_vol += v[i];
460 
461  v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
462  if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
463 
464  zone_volumes.push_back(v1[i]*piTimes4thirds); // True volume of shell
465  }
466 
467  return tot_vol; // Sum of zone integrals, not geometric volume
468 }
G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2, G4double nuclearRadius) const
G4double zoneIntegralGaussian(G4double ur1, G4double ur2, G4double nuclearRadius) const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

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.

1301  {
1302  return (isProjectile(cparticle) &&
1303  (cparticle.getParticle().isPhoton() ||
1304  cparticle.getParticle().isMuon())
1305  );
1306 }
const G4InuclElementaryParticle & getParticle() const
G4bool isProjectile(const G4CascadParticle &cparticle) const

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.

1854  {
1855  // Delay interactions of newly formed secondaries (minimum int. length)
1856  const G4double young_cut = std::sqrt(10.0) * 0.25;
1857  const G4double huge_num = 50.0; // Argument to exponential
1858 
1859  G4double spath = large; // Buffer for return value
1860 
1861  if (invmfp < small) return spath; // No interaction, avoid unnecessary work
1862 
1863  G4double pw = -path * invmfp; // Ratio of path in zone to MFP
1864  if (pw < -huge_num) pw = -huge_num;
1865  pw = 1.0 - G4Exp(pw);
1866 
1867  if (verboseLevel > 2)
1868  G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
1869 
1870  // Primary particle(s) should always interact at least once
1871  if (forceFirst(cparticle) || (inuclRndm() < pw)) {
1872  spath = -G4Log(1.0 - pw * inuclRndm()) / invmfp;
1873  if (cparticle.young(young_cut, spath)) spath = large;
1874 
1875  if (verboseLevel > 2)
1876  G4cout << " spath " << spath << " path " << path << G4endl;
1877  }
1878 
1879  return spath;
1880 }
G4bool young(G4double young_path_cut, G4double cpath) const
G4GLOB_DLL std::ostream G4cout
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4bool forceFirst(const G4CascadParticle &cparticle) const

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.

686  {
687  if (verboseLevel > 1) {
688  G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
689  }
690 
691  thePartners.clear(); // Reset buffer for next cycle
692 
693  G4int ptype = cparticle.getParticle().type();
694  G4int zone = cparticle.getCurrentZone();
695 
696  G4double r_in; // Destination radius if inbound
697  G4double r_out; // Destination radius if outbound
698 
699  if (zone == number_of_zones) {
700  r_in = nuclei_radius;
701  r_out = 0.0;
702  } else if (zone == 0) { // particle is outside core
703  r_in = 0.0;
704  r_out = zone_radii[0];
705  } else {
706  r_in = zone_radii[zone - 1];
707  r_out = zone_radii[zone];
708  }
709 
710  G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
711 
712  if (verboseLevel > 2) {
713  if (isProjectile(cparticle)) G4cout << " incident particle: ";
714  G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
715  << G4endl;
716  }
717 
718  if (path < -small) { // something wrong
719  if (verboseLevel)
720  G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
721  return;
722  }
723 
724  if (std::fabs(path) < small) { // Not moving, or just at boundary
725  if (cparticle.getMomentum().vect().mag() > small) {
726  if (verboseLevel > 3)
727  G4cout << " generateInteractionPartners-> zero path" << G4endl;
728 
729  thePartners.push_back(partner()); // Dummy list terminator with zero path
730  return;
731  }
732 
733  if (zone >= number_of_zones) // Place captured-at-rest in nucleus
734  zone = number_of_zones-1;
735  }
736 
737  G4double invmfp = 0.; // Buffers for interaction probability
738  G4double spath = 0.;
739  for (G4int ip = 1; ip < 3; ip++) {
740  // Only process nucleons which remain active in target
741  if (ip==proton && protonNumberCurrent < 1) continue;
742  if (ip==neutron && neutronNumberCurrent < 1) continue;
743  if (ip==neutron && ptype==muonMinus) continue; // mu-/n forbidden
744 
745  // All nucleons are assumed to be at rest when colliding
746  G4InuclElementaryParticle particle = generateNucleon(ip, zone);
747  invmfp = inverseMeanFreePath(cparticle, particle);
748  spath = generateInteractionLength(cparticle, path, invmfp);
749 
750  if (path<small || spath < path) {
751  if (verboseLevel > 3) {
752  G4cout << " adding partner[" << thePartners.size() << "]: "
753  << particle << G4endl;
754  }
755  thePartners.push_back(partner(particle, spath));
756  }
757  } // for (G4int ip...
758 
759  if (verboseLevel > 2) {
760  G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
761  }
762 
763  // Absorption possible for pions or photons interacting with dibaryons
764  if (useQuasiDeuteron(cparticle.getParticle().type())) {
765  if (verboseLevel > 2) {
766  G4cout << " trying quasi-deuterons with bullet: "
767  << cparticle.getParticle() << G4endl;
768  }
769 
770  // Initialize buffers for quasi-deuteron results
771  qdeutrons.clear();
772  acsecs.clear();
773 
774  G4double tot_invmfp = 0.0; // Total inv. mean-free-path for all QDs
775 
776  // Proton-proton state interacts with pi-, mu- or neutrals
777  if (protonNumberCurrent >= 2 && ptype != pip) {
779  if (verboseLevel > 2)
780  G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
781 
782  invmfp = inverseMeanFreePath(cparticle, ppd);
783  tot_invmfp += invmfp;
784  acsecs.push_back(invmfp);
785  qdeutrons.push_back(ppd);
786  }
787 
788  // Proton-neutron state interacts with any pion type or photon
789  if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
791  if (verboseLevel > 2)
792  G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
793 
794  invmfp = inverseMeanFreePath(cparticle, npd);
795  tot_invmfp += invmfp;
796  acsecs.push_back(invmfp);
797  qdeutrons.push_back(npd);
798  }
799 
800  // Neutron-neutron state interacts with pi+ or neutrals
801  if (neutronNumberCurrent >= 2 && ptype != pim && ptype != mum) {
803  if (verboseLevel > 2)
804  G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
805 
806  invmfp = inverseMeanFreePath(cparticle, nnd);
807  tot_invmfp += invmfp;
808  acsecs.push_back(invmfp);
809  qdeutrons.push_back(nnd);
810  }
811 
812  // Select quasideuteron interaction from non-zero cross-section choices
813  if (verboseLevel > 2) {
814  for (size_t i=0; i<qdeutrons.size(); i++) {
815  G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
816  << "] " << acsecs[i];
817  }
818  G4cout << G4endl;
819  }
820 
821  if (tot_invmfp > small) { // Must have absorption cross-section
822  G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
823 
824  if (path<small || apath < path) { // choose the qdeutron
825  G4double sl = inuclRndm() * tot_invmfp;
826  G4double as = 0.0;
827 
828  for (size_t i = 0; i < qdeutrons.size(); i++) {
829  as += acsecs[i];
830  if (sl < as) {
831  if (verboseLevel > 2)
832  G4cout << " deut type " << qdeutrons[i] << G4endl;
833 
834  thePartners.push_back(partner(qdeutrons[i], apath));
835  break;
836  }
837  } // for (qdeutrons...
838  } // if (apath < path)
839  } // if (tot_invmfp > small)
840  } // if (useQuasiDeuteron) [pion, muon or photon]
841 
842  if (verboseLevel > 2) {
843  G4cout << " after deuterons " << thePartners.size() << " partners"
844  << G4endl;
845  }
846 
847  if (thePartners.size() > 1) { // Sort list by path length
848  std::sort(thePartners.begin(), thePartners.end(), sortPartners);
849  }
850 
851  G4InuclElementaryParticle particle; // Total path at end of list
852  thePartners.push_back(partner(particle, path));
853 }
G4InuclElementaryParticle generateQuasiDeuteron(G4int type1, G4int type2, G4int zone) const
G4LorentzVector getMomentum() const
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
G4double getPathToTheNextZone(G4double rz_in, G4double rz_out)
std::vector< partner > thePartners
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
int G4int
Definition: G4Types.hh:78
const G4InuclElementaryParticle & getParticle() const
G4bool isProjectile(const G4CascadParticle &cparticle) const
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
G4double generateInteractionLength(const G4CascadParticle &cparticle, G4double path, G4double invmfp) const
std::pair< G4InuclElementaryParticle, G4double > partner
static G4bool sortPartners(const partner &p1, const partner &p2)
G4int getCurrentZone() const
#define G4endl
Definition: G4ios.hh:61
G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const
double G4double
Definition: G4Types.hh:76
double mag() const
G4GLOB_DLL std::ostream G4cerr

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.

303  {
304  generateModel(nuclei->getA(), nuclei->getZ());
305 }
G4int getZ() const
G4int getA() const
void generateModel(G4InuclNuclei *nuclei)

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.

307  {
308  if (verboseLevel) {
309  G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
310  << G4endl;
311  }
312 
313  // If model already built, just return; otherwise intialize everything
314  if (a == A && z == Z) {
315  if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
316  reset();
317  return;
318  }
319 
320  A = a;
321  Z = z;
322  delete theNucleus;
323  theNucleus = new G4InuclNuclei(A,Z); // For conservation checking
324 
325  neutronNumber = A - Z;
326  protonNumber = Z;
327  reset();
328 
329  if (verboseLevel > 3) {
330  G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl
331  << " radiusUnits = " << radiusUnits << G4endl
332  << " skinDepth = " << skinDepth << G4endl
333  << " radiusScale = " << radiusScale << G4endl
334  << " radiusScale2 = " << radiusScale2 << G4endl
335  << " radiusForSmall = " << radiusForSmall << G4endl
336  << " radScaleAlpha = " << radScaleAlpha << G4endl
337  << " fermiMomentum = " << fermiMomentum << G4endl
338  << " piTimes4thirds = " << piTimes4thirds << G4endl;
339  }
340 
341  G4double nuclearRadius; // Nuclear radius computed from A
342  if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
343  else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
344 
345  // This will be used to pre-allocate lots of arrays below
346  number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
347 
348  // Clear all parameters arrays for reloading
349  binding_energies.clear();
350  nucleon_densities.clear();
351  zone_potentials.clear();
352  fermi_momenta.clear();
353  zone_radii.clear();
354  zone_volumes.clear();
355 
357  fillZoneRadii(nuclearRadius);
358 
359  G4double tot_vol = fillZoneVolumes(nuclearRadius); // Woods-Saxon integral
360 
361  fillPotentials(proton, tot_vol); // Protons
362  fillPotentials(neutron, tot_vol); // Neutrons
363 
364  // Additional flat zone potentials for other hadrons
365  const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
366  const std::vector<G4double> kp(number_of_zones, kaon_vp);
367  const std::vector<G4double> hp(number_of_zones, hyperon_vp);
368 
369  zone_potentials.push_back(vp);
370  zone_potentials.push_back(kp);
371  zone_potentials.push_back(hp);
372 
373  nuclei_radius = zone_radii.back();
374  nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
375 
376  if (verboseLevel > 3) printModel();
377 }
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void fillBindingEnergies()
void fillZoneRadii(G4double nuclearRadius)
G4GLOB_DLL std::ostream G4cout
void fillPotentials(G4int type, G4double tot_vol)
G4double fillZoneVolumes(G4double nuclearRadius)
tuple z
Definition: test.py:28
void printModel() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Definition at line 649 of file G4NucleiModel.cc.

649  {
650  if (verboseLevel > 1) {
651  G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
652  }
653 
654  G4LorentzVector mom = generateNucleonMomentum(type, zone);
655  return G4InuclElementaryParticle(mom, type);
656 }
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

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.

640  {
641  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
643 
644  return generateWithRandomAngles(pmod, mass);
645 }
static G4double getParticleMass(G4int type)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double getFermiMomentum(G4int ip, G4int izone) const
double G4double
Definition: G4Types.hh:76

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.

859  {
860  if (verboseLevel > 1)
861  G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
862 
863  if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
864 
865  // Create four-vector checking
866 #ifdef G4CASCADE_CHECK_ECONS
867  G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV
868  balance.setVerboseLevel(verboseLevel);
869 #endif
870 
871  outgoing_cparticles.clear(); // Clear return buffer for this event
872  generateInteractionPartners(cparticle); // Fills "thePartners" data
873 
874  if (thePartners.empty()) { // smth. is wrong -> needs special treatment
875  if (verboseLevel)
876  G4cerr << " generateParticleFate-> got empty interaction-partners list "
877  << G4endl;
878  return;
879  }
880 
881  G4int npart = thePartners.size(); // Last item is a total-path placeholder
882 
883  if (npart == 1) { // cparticle is on the next zone entry
884  if (verboseLevel > 1)
885  G4cout << " no interactions; moving to next zone" << G4endl;
886 
888  cparticle.incrementCurrentPath(thePartners[0].second);
889  boundaryTransition(cparticle);
890  outgoing_cparticles.push_back(cparticle);
891 
892  if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
893 
894  //A.R. 19-Jun-2013: Fixed rare cases of non-reproducibility.
895  current_nucl1 = 0;
896  current_nucl2 = 0;
897 
898  return;
899  } // if (npart == 1)
900 
901  if (verboseLevel > 1)
902  G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
903 
904  G4ThreeVector old_position = cparticle.getPosition();
905  G4InuclElementaryParticle& bullet = cparticle.getParticle();
906  G4bool no_interaction = true;
907  G4int zone = cparticle.getCurrentZone();
908 
909  for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder
910  if (i > 0) cparticle.updatePosition(old_position);
911 
913 
914  if (verboseLevel > 3) {
915  if (target.quasi_deutron()) G4cout << " try absorption: ";
916  G4cout << " target " << target.type() << " bullet " << bullet.type()
917  << G4endl;
918  }
919 
920  EPCoutput.reset();
921  // Pass current (A,Z) configuration for possible recoils
922  G4int massNumberCurrent = protonNumberCurrent+neutronNumberCurrent;
923  theEPCollider->setNucleusState(massNumberCurrent, protonNumberCurrent);
924  theEPCollider->collide(&bullet, &target, EPCoutput);
925 
926  // If collision failed, exit loop over partners
927  if (EPCoutput.numberOfOutgoingParticles() == 0) break;
928 
929  if (verboseLevel > 2) {
930  EPCoutput.printCollisionOutput();
931 #ifdef G4CASCADE_CHECK_ECONS
932  balance.collide(&bullet, &target, EPCoutput);
933  balance.okay(); // Do checks, but ignore result
934 #endif
935  }
936 
937  // Get list of outgoing particles for evaluation
938  std::vector<G4InuclElementaryParticle>& outgoing_particles =
939  EPCoutput.getOutgoingParticles();
940 
941  if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
942 
943  // Trailing effect: reject interaction at previously hit nucleon
944  cparticle.propagateAlongThePath(thePartners[i].second);
945  const G4ThreeVector& new_position = cparticle.getPosition();
946 
947  if (!passTrailing(new_position)) continue;
948  collisionPts.push_back(new_position);
949 
950  // Sort particles according to beta (fastest first)
951  std::sort(outgoing_particles.begin(), outgoing_particles.end(),
953 
954  if (verboseLevel > 2)
955  G4cout << " adding " << outgoing_particles.size()
956  << " output particles" << G4endl;
957 
958  // NOTE: Embedded temporary is optimized away (no copying gets done)
959  G4int nextGen = cparticle.getGeneration()+1;
960  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) {
961  outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
962  new_position, zone,
963  0.0, nextGen));
964  }
965 
966  no_interaction = false;
967  current_nucl1 = 0;
968  current_nucl2 = 0;
969 
970 #ifdef G4CASCADE_DEBUG_CHARGE
971  {
972  G4double out_charge = 0.0;
973 
974  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
975  out_charge += outgoing_particles[ip].getCharge();
976 
977  G4cout << " multiplicity " << outgoing_particles.size()
978  << " bul type " << bullet.type()
979  << " targ type " << target.type()
980  << "\n initial charge "
981  << bullet.getCharge() + target.getCharge()
982  << " out charge " << out_charge << G4endl;
983  }
984 #endif
985 
986  if (verboseLevel > 2)
987  G4cout << " partner type " << target.type() << G4endl;
988 
989  if (target.nucleon()) {
990  current_nucl1 = target.type();
991  } else {
992  if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
993  current_nucl1 = (target.type() - 100) / 10;
994  current_nucl2 = target.type() - 100 - 10 * current_nucl1;
995  }
996 
997  if (current_nucl1 == 1) {
998  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
999  protonNumberCurrent--;
1000  } else {
1001  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
1002  neutronNumberCurrent--;
1003  }
1004 
1005  if (current_nucl2 == 1) {
1006  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
1007  protonNumberCurrent--;
1008  } else if (current_nucl2 == 2) {
1009  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
1010  neutronNumberCurrent--;
1011  }
1012 
1013  break;
1014  } // loop over partners
1015 
1016  if (no_interaction) { // still no interactions
1017  if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
1018 
1019  // For conservation checking (below), get particle before updating
1020  static G4ThreadLocal G4InuclElementaryParticle *prescatCP_G4MT_TLS_ = 0;
1021  if (!prescatCP_G4MT_TLS_) {
1022  prescatCP_G4MT_TLS_ = new G4InuclElementaryParticle;
1023  G4AutoDelete::Register(prescatCP_G4MT_TLS_);
1024  }
1025  G4InuclElementaryParticle &prescatCP = *prescatCP_G4MT_TLS_; // Avoid memory churn
1026  prescatCP = cparticle.getParticle();
1027 
1028  // Last "partner" is just a total-path placeholder
1029  cparticle.updatePosition(old_position);
1030  cparticle.propagateAlongThePath(thePartners[npart-1].second);
1031  cparticle.incrementCurrentPath(thePartners[npart-1].second);
1032  boundaryTransition(cparticle);
1033  outgoing_cparticles.push_back(cparticle);
1034 
1035  // Check conservation for simple scattering (ignore target nucleus!)
1036 #ifdef G4CASCADE_CHECK_ECONS
1037  if (verboseLevel > 2) {
1038  balance.collide(&prescatCP, 0, outgoing_cparticles);
1039  balance.okay(); // Report violations, but don't act on them
1040  }
1041 #endif
1042  } // if (no_interaction)
1043 
1044  return;
1045 }
void incrementCurrentPath(G4double npath)
const XML_Char * target
Definition: expat.h:268
void generateInteractionPartners(G4CascadParticle &cparticle)
G4int getGeneration() const
void printCollisionOutput(std::ostream &os=G4cout) const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4bool passFermi(const std::vector< G4InuclElementaryParticle > &particles, G4int zone)
void propagateAlongThePath(G4double path)
std::vector< partner > thePartners
void updatePosition(const G4ThreeVector &pos)
static constexpr double second
Definition: G4SIunits.hh:157
void boundaryTransition(G4CascadParticle &cparticle)
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
const G4InuclElementaryParticle & getParticle() const
void Register(T *inst)
Definition: G4AutoDelete.hh:65
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4int numberOfOutgoingParticles() const
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
G4int getCurrentZone() const
G4double getCharge() const
#define G4endl
Definition: G4ios.hh:61
const G4ThreeVector & getPosition() const
double G4double
Definition: G4Types.hh:76
G4bool passTrailing(const G4ThreeVector &hit_position)
G4GLOB_DLL std::ostream G4cerr

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.

661  {
662  if (verboseLevel > 1) {
663  G4cout << " >>> G4NucleiModel::generateQuasiDeuteron" << G4endl;
664  }
665 
666  // Quasideuteron consists of an unbound but associated nucleon pair
667 
668  // FIXME: Why generate two separate nucleon momenta (randomly!) and
669  // add them, instead of just throwing a net momentum for the
670  // dinulceon state? And why do I have to capture the two
671  // return values into local variables?
672  G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
673  G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
674  G4LorentzVector dmom = mom1+mom2;
675 
676  G4int dtype = 0;
677  if (type1*type2 == pro*pro) dtype = 111;
678  else if (type1*type2 == pro*neu) dtype = 112;
679  else if (type1*type2 == neu*neu) dtype = 122;
680 
681  return G4InuclElementaryParticle(dmom, dtype);
682 }
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

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.

1360  {
1361  const G4double pn_spec = 1.0; // Scale factor for pn vs. pp/nn
1362  //const G4double pn_spec = 0.5;
1363 
1364  G4double dens = 0.;
1365 
1366  if (ip < 100) dens = getDensity(ip,izone);
1367  else { // For dibaryons, remove extra 1/volume term in density product
1368  switch (ip) {
1369  case diproton:
1370  dens = getDensity(proton,izone) * getDensity(proton,izone);
1371  break;
1372  case unboundPN:
1373  dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
1374  break;
1375  case dineutron:
1376  dens = getDensity(neutron,izone) * getDensity(neutron,izone);
1377  break;
1378  default: dens = 0.;
1379  }
1380  dens *= getVolume(izone);
1381  }
1382 
1383  return getRatio(ip) * dens;
1384 }
G4double getVolume(G4int izone) const
G4double getRatio(G4int ip) const
double G4double
Definition: G4Types.hh:76
G4double getDensity(G4int ip, G4int izone) const

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.

112  {
113  return nucleon_densities[ip - 1][izone];
114  }

Here is the caller graph for this function:

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

Definition at line 626 of file G4NucleiModel.cc.

626  {
627  G4double ekin = 0.0;
628 
629  if (ip < 3 && izone < number_of_zones) { // ip for proton/neutron only
630  G4double pfermi = fermi_momenta[ip - 1][izone];
632 
633  ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
634  }
635  return ekin;
636 }
static G4double getParticleMass(G4int type)
double G4double
Definition: G4Types.hh:76

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.

116  {
117  return fermi_momenta[ip - 1][izone];
118  }

Here is the caller graph for this function:

G4int G4NucleiModel::getNumberOfNeutrons ( ) const
inline

Definition at line 149 of file G4NucleiModel.hh.

149 { return neutronNumberCurrent; }
G4int G4NucleiModel::getNumberOfProtons ( ) const
inline

Definition at line 150 of file G4NucleiModel.hh.

150 { return protonNumberCurrent; }
G4int G4NucleiModel::getNumberOfZones ( ) const
inline

Definition at line 143 of file G4NucleiModel.hh.

143 { return number_of_zones; }
G4double G4NucleiModel::getPotential ( G4int  ip,
G4int  izone 
) const
inline

Definition at line 122 of file G4NucleiModel.hh.

122  {
123  if (ip == 9 || ip < 0) return 0.0; // Photons and leptons
124  G4int ip0 = ip < 3 ? ip - 1 : 2;
125  if (ip > 10 && ip < 18) ip0 = 3;
126  if (ip > 20) ip0 = 4;
127  return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
128  }
int G4int
Definition: G4Types.hh:78

Here is the caller graph for this function:

G4double G4NucleiModel::getRadius ( ) const
inline

Definition at line 133 of file G4NucleiModel.hh.

133 { return nuclei_radius; }
G4double G4NucleiModel::getRadius ( G4int  izone) const
inline

Definition at line 134 of file G4NucleiModel.hh.

134  {
135  return ( (izone<0) ? 0.
136  : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
137  }
G4double G4NucleiModel::getRadiusUnits ( ) const
inline

Definition at line 131 of file G4NucleiModel.hh.

131 { return radiusUnits*CLHEP::fermi; }
static constexpr double fermi
Definition: SystemOfUnits.h:83
G4double G4NucleiModel::getRatio ( G4int  ip) const
protected

Definition at line 1343 of file G4NucleiModel.cc.

1343  {
1344  if (verboseLevel > 4) {
1345  G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
1346  }
1347 
1348  switch (ip) {
1349  case proton: return G4double(protonNumberCurrent)/G4double(protonNumber);
1350  case neutron: return G4double(neutronNumberCurrent)/G4double(neutronNumber);
1351  case diproton: return getRatio(proton)*getRatio(proton);
1352  case unboundPN: return getRatio(proton)*getRatio(neutron);
1353  case dineutron: return getRatio(neutron)*getRatio(neutron);
1354  default: return 0.;
1355  }
1356 
1357  return 0.;
1358 }
G4GLOB_DLL std::ostream G4cout
G4double getRatio(G4int ip) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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

Definition at line 169 of file G4NucleiModel.hh.

169  {
170  return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
171  }
G4double G4NucleiModel::getVolume ( G4int  izone) const
inline

Definition at line 138 of file G4NucleiModel.hh.

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

Here is the caller graph for this function:

G4int G4NucleiModel::getZone ( G4double  r) const
inline

Definition at line 144 of file G4NucleiModel.hh.

144  {
145  for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
146  return number_of_zones;
147  }
int G4int
Definition: G4Types.hh:78

Here is the caller graph for this function:

G4CascadParticle G4NucleiModel::initializeCascad ( G4InuclElementaryParticle particle)

Definition at line 1388 of file G4NucleiModel.cc.

1388  {
1389  if (verboseLevel > 1) {
1390  G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
1391  }
1392 
1393  // FIXME: Previous version generated random sin(theta), then used -cos(theta)
1394  // Using generateWithRandomAngles changes result!
1395  // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
1396  G4double costh = std::sqrt(1.0 - inuclRndm());
1397  G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
1398 
1399  // Start particle outside nucleus, unless capture-at-rest
1400  G4int zone = number_of_zones;
1401  if (particle->getKineticEnergy() < small) zone--;
1402 
1403  G4CascadParticle cpart(*particle, pos, zone, large, 0);
1404 
1405  // SPECIAL CASE: Inbound photons are emplanted along through-path
1406  if (forceFirst(cpart)) choosePointAlongTraj(cpart);
1407 
1408  if (verboseLevel > 2) G4cout << cpart << G4endl;
1409 
1410  return cpart;
1411 }
int G4int
Definition: G4Types.hh:78
G4double getKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
void choosePointAlongTraj(G4CascadParticle &cparticle)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static const G4double pos
G4bool forceFirst(const G4CascadParticle &cparticle) const

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.

1415  {
1416  if (verboseLevel) {
1417  G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
1418  << G4endl;
1419  }
1420 
1421  const G4double max_a_for_cascad = 5.0;
1422  const G4double ekin_cut = 2.0;
1423  const G4double small_ekin = 1.0e-6;
1424  const G4double r_large2for3 = 62.0;
1425  const G4double r0forAeq3 = 3.92;
1426  const G4double s3max = 6.5;
1427  const G4double r_large2for4 = 69.14;
1428  const G4double r0forAeq4 = 4.16;
1429  const G4double s4max = 7.0;
1430  const G4int itry_max = 100;
1431 
1432  // Convenient references to input buffer contents
1433  std::vector<G4CascadParticle>& casparticles = output.first;
1434  std::vector<G4InuclElementaryParticle>& particles = output.second;
1435 
1436  casparticles.clear(); // Reset input buffer to fill this event
1437  particles.clear();
1438 
1439  // first decide whether it will be cascad or compound final nuclei
1440  G4int ab = bullet->getA();
1441  G4int zb = bullet->getZ();
1442  G4int at = target->getA();
1443  G4int zt = target->getZ();
1444 
1445  G4double massb = bullet->getMass(); // For creating LorentzVectors below
1446 
1447  if (ab < max_a_for_cascad) {
1448 
1449  G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
1450  G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
1451  G4double ben = benb < bent ? bent : benb;
1452 
1453  if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
1454  G4int itryg = 0;
1455 
1456  /* Loop checking 08.06.2015 MHK */
1457  while (casparticles.size() == 0 && itryg < itry_max) {
1458  itryg++;
1459  particles.clear();
1460 
1461  // nucleons coordinates and momenta in nuclei rest frame
1462  coordinates.clear();
1463  momentums.clear();
1464 
1465  if (ab < 3) { // deuteron, simplest case
1466  G4double r = 2.214 - 3.4208 * G4Log(1.0 - 0.981 * inuclRndm());
1468  coordinates.push_back(coord1);
1469  coordinates.push_back(-coord1);
1470 
1471  G4double p = 0.0;
1472  G4bool bad = true;
1473  G4int itry = 0;
1474 
1475  while (bad && itry < itry_max) { /* Loop checking 08.06.2015 MHK */
1476  itry++;
1477  p = 456.0 * inuclRndm();
1478 
1479  if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1480  p * r > 312.0) bad = false;
1481  }
1482 
1483  if (itry == itry_max)
1484  if (verboseLevel > 2){
1485  G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;
1486  }
1487 
1488  p = 0.0005 * p;
1489 
1490  if (verboseLevel > 2){
1491  G4cout << " p nuc " << p << G4endl;
1492  }
1493 
1494  G4LorentzVector mom = generateWithRandomAngles(p, massb);
1495 
1496  momentums.push_back(mom);
1497  mom.setVect(-mom.vect());
1498  momentums.push_back(-mom);
1499  } else {
1500  G4ThreeVector coord1;
1501 
1502  G4bool badco = true;
1503 
1504  G4int itry = 0;
1505 
1506  if (ab == 3) {
1507  while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
1508  if (itry > 0) coordinates.clear();
1509  itry++;
1510  G4int i(0);
1511 
1512  for (i = 0; i < 2; i++) {
1513  G4int itry1 = 0;
1514  G4double ss, u, rho;
1515  G4double fmax = G4Exp(-0.5) / std::sqrt(0.5);
1516 
1517  while (itry1 < itry_max) { /* Loop checking 08.06.2015 MHK */
1518  itry1++;
1519  ss = -G4Log(inuclRndm());
1520  u = fmax * inuclRndm();
1521  rho = std::sqrt(ss) * G4Exp(-ss);
1522 
1523  if (rho > u && ss < s3max) {
1524  ss = r0forAeq3 * std::sqrt(ss);
1525  coord1 = generateWithRandomAngles(ss).vect();
1526  coordinates.push_back(coord1);
1527 
1528  if (verboseLevel > 2){
1529  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1530  }
1531  break;
1532  }
1533  }
1534 
1535  if (itry1 == itry_max) { // bad case
1536  coord1.set(10000.,10000.,10000.);
1537  coordinates.push_back(coord1);
1538  break;
1539  }
1540  }
1541 
1542  coord1 = -coordinates[0] - coordinates[1];
1543  if (verboseLevel > 2) {
1544  G4cout << " 3 r " << coord1.mag() << G4endl;
1545  }
1546 
1547  coordinates.push_back(coord1);
1548 
1549  G4bool large_dist = false;
1550 
1551  for (i = 0; i < 2; i++) {
1552  for (G4int j = i+1; j < 3; j++) {
1553  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1554 
1555  if (verboseLevel > 2) {
1556  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1557  }
1558 
1559  if (r2 > r_large2for3) {
1560  large_dist = true;
1561 
1562  break;
1563  }
1564  }
1565 
1566  if (large_dist) break;
1567  }
1568 
1569  if(!large_dist) badco = false;
1570 
1571  }
1572 
1573  } else { // a >= 4
1574  G4double b = 3./(ab - 2.0);
1575  G4double b1 = 1.0 - b / 2.0;
1576  G4double u = b1 + std::sqrt(b1 * b1 + b);
1577  G4double fmax = (1.0 + u/b) * u * G4Exp(-u);
1578 
1579  while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
1580 
1581  if (itry > 0) coordinates.clear();
1582  itry++;
1583  G4int i(0);
1584 
1585  for (i = 0; i < ab-1; i++) {
1586  G4int itry1 = 0;
1587  G4double ss;
1588 
1589  while (itry1 < itry_max) { /* Loop checking 08.06.2015 MHK */
1590  itry1++;
1591  ss = -G4Log(inuclRndm());
1592  u = fmax * inuclRndm();
1593 
1594  if (std::sqrt(ss) * G4Exp(-ss) * (1.0 + ss/b) > u
1595  && ss < s4max) {
1596  ss = r0forAeq4 * std::sqrt(ss);
1597  coord1 = generateWithRandomAngles(ss).vect();
1598  coordinates.push_back(coord1);
1599 
1600  if (verboseLevel > 2) {
1601  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1602  }
1603 
1604  break;
1605  }
1606  }
1607 
1608  if (itry1 == itry_max) { // bad case
1609  coord1.set(10000.,10000.,10000.);
1610  coordinates.push_back(coord1);
1611  break;
1612  }
1613  }
1614 
1615  coord1 *= 0.0; // Cheap way to reset
1616  for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
1617 
1618  coordinates.push_back(coord1);
1619 
1620  if (verboseLevel > 2){
1621  G4cout << " last r " << coord1.mag() << G4endl;
1622  }
1623 
1624  G4bool large_dist = false;
1625 
1626  for (i = 0; i < ab-1; i++) {
1627  for (G4int j = i+1; j < ab; j++) {
1628 
1629  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1630 
1631  if (verboseLevel > 2){
1632  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1633  }
1634 
1635  if (r2 > r_large2for4) {
1636  large_dist = true;
1637 
1638  break;
1639  }
1640  }
1641 
1642  if (large_dist) break;
1643  }
1644 
1645  if (!large_dist) badco = false;
1646  }
1647  }
1648 
1649  if(badco) {
1650  G4cout << " can not generate the nucleons coordinates for a "
1651  << ab << G4endl;
1652 
1653  casparticles.clear(); // Return empty buffer on error
1654  particles.clear();
1655  return;
1656 
1657  } else { // momentums
1658  G4double p, u, x;
1659  G4LorentzVector mom;
1660  //G4bool badp = True;
1661 
1662  for (G4int i = 0; i < ab - 1; i++) {
1663  G4int itry2 = 0;
1664 
1665  while(itry2 < itry_max) { /* Loop checking 08.06.2015 MHK */
1666  itry2++;
1667  u = -G4Log(0.879853 - 0.8798502 * inuclRndm());
1668  x = u * G4Exp(-u);
1669 
1670  if(x > inuclRndm()) {
1671  p = std::sqrt(0.01953 * u);
1672  mom = generateWithRandomAngles(p, massb);
1673  momentums.push_back(mom);
1674 
1675  break;
1676  }
1677  } // while (itry2 < itry_max)
1678 
1679  if(itry2 == itry_max) {
1680  G4cout << " can not generate proper momentum for a "
1681  << ab << G4endl;
1682 
1683  casparticles.clear(); // Return empty buffer on error
1684  particles.clear();
1685  return;
1686  }
1687  } // for (i=0 ...
1688  // last momentum
1689 
1690  mom *= 0.; // Cheap way to reset
1691  mom.setE(bullet->getEnergy()+target->getEnergy());
1692 
1693  for(G4int j=0; j< ab-1; j++) mom -= momentums[j];
1694 
1695  momentums.push_back(mom);
1696  }
1697  }
1698 
1699  // Coordinates and momenta at rest are generated, now back to the lab
1700  G4double rb = 0.0;
1701  G4int i(0);
1702 
1703  for(i = 0; i < G4int(coordinates.size()); i++) {
1704  G4double rp = coordinates[i].mag2();
1705 
1706  if(rp > rb) rb = rp;
1707  }
1708 
1709  // nuclei i.p. as a whole
1710  G4double s1 = std::sqrt(inuclRndm());
1711  G4double phi = randomPHI();
1712  G4double rz = (nuclei_radius + rb) * s1;
1713  G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
1714  -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
1715 
1716  for (i = 0; i < G4int(coordinates.size()); i++) {
1717  coordinates[i] += global_pos;
1718  }
1719 
1720  // all nucleons at rest
1721  raw_particles.clear();
1722 
1723  for (G4int ipa = 0; ipa < ab; ipa++) {
1724  G4int knd = ipa < zb ? 1 : 2;
1725  raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1726  }
1727 
1728  G4InuclElementaryParticle dummy(small_ekin, 1);
1729  G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
1730  toTheBulletRestFrame.toTheTargetRestFrame();
1731 
1732  particleIterator ipart;
1733 
1734  for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1735  ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum()));
1736  }
1737 
1738  // fill cascad particles and outgoing particles
1739 
1740  for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1741  G4LorentzVector mom = raw_particles[ip].getMomentum();
1742  G4double pmod = mom.rho();
1743  G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
1744  G4double det = t0 * t0 + nuclei_radius * nuclei_radius
1745  - coordinates[ip].mag2();
1746  G4double tr = -1.0;
1747 
1748  if(det > 0.0) {
1749  G4double t1 = t0 + std::sqrt(det);
1750  G4double t2 = t0 - std::sqrt(det);
1751 
1752  if(std::fabs(t1) <= std::fabs(t2)) {
1753  if(t1 > 0.0) {
1754  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1755  }
1756 
1757  if(tr < 0.0 && t2 > 0.0) {
1758 
1759  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1760  }
1761 
1762  } else {
1763  if(t2 > 0.0) {
1764 
1765  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1766  }
1767 
1768  if(tr < 0.0 && t1 > 0.0) {
1769  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1770  }
1771  }
1772 
1773  }
1774 
1775  if(tr >= 0.0) { // cascad particle
1776  coordinates[ip] += mom.vect()*tr / pmod;
1777  casparticles.push_back(G4CascadParticle(raw_particles[ip],
1778  coordinates[ip],
1779  number_of_zones, large, 0));
1780 
1781  } else {
1782  particles.push_back(raw_particles[ip]);
1783  }
1784  }
1785  }
1786 
1787  if(casparticles.size() == 0) {
1788  particles.clear();
1789 
1790  G4cout << " can not generate proper distribution for " << itry_max
1791  << " steps " << G4endl;
1792  }
1793  }
1794  }
1795 
1796  if(verboseLevel > 2){
1797  G4int ip(0);
1798 
1799  G4cout << " cascad particles: " << casparticles.size() << G4endl;
1800  for(ip = 0; ip < G4int(casparticles.size()); ip++)
1801  G4cout << casparticles[ip] << G4endl;
1802 
1803  G4cout << " outgoing particles: " << particles.size() << G4endl;
1804  for(ip = 0; ip < G4int(particles.size()); ip++)
1805  G4cout << particles[ip] << G4endl;
1806  }
1807 
1808  return; // Buffer has been filled
1809 }
void set(double x, double y, double z)
G4int getZ() const
double dot(const Hep3Vector &) const
const char * p
Definition: xmltok.h:285
G4double getEnergy() const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double getKineticEnergy() const
tuple b
Definition: test.py:12
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
G4int getA() const
bool G4bool
Definition: G4Types.hh:79
double rho() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const G4double ab
tuple t1
Definition: plottest35.py:33
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
#define G4endl
Definition: G4ios.hh:61
void setVect(const Hep3Vector &)
double G4double
Definition: G4Types.hh:76
G4double bindingEnergy(G4int A, G4int Z)
double mag() const
G4double getMass() const

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.

1816  {
1817  G4int ptype = cparticle.getParticle().type();
1818  G4int ip = target.type();
1819 
1820  // Ensure that zone specified is within nucleus, for array lookups
1821  if (zone<0) zone = cparticle.getCurrentZone();
1822  if (zone>=number_of_zones) zone = number_of_zones-1;
1823 
1824  // Special cases: neutrinos, and muon-on-neutron, have infinite path
1825  if (cparticle.getParticle().isNeutrino()) return 0.;
1826  if (ptype == muonMinus && ip == neutron) return 0.;
1827 
1828  // Configure frame transformation to get kinetic energy for lookups
1829  dummy_convertor.setBullet(cparticle.getParticle());
1830  dummy_convertor.setTarget(&target);
1831  dummy_convertor.toTheCenterOfMass(); // Fill internal kinematics
1832  G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
1833 
1834  // Get cross section for interacting with target (dibaryons are absorptive)
1835  G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
1836  : absorptionCrossSection(ekin, ptype);
1837 
1838  if (verboseLevel > 2) {
1839  G4cout << " ip " << ip << " zone " << zone << " ekin " << ekin
1840  << " dens " << getCurrentDensity(ip, zone)
1841  << " csec " << csec << G4endl;
1842  }
1843 
1844  if (csec <= 0.) return 0.; // No interaction, avoid unnecessary work
1845 
1846  // Get nuclear density and compute mean free path
1847  return csec * getCurrentDensity(ip, zone);
1848 }
void setBullet(const G4InuclParticle *bullet)
int G4int
Definition: G4Types.hh:78
G4double absorptionCrossSection(G4double e, G4int type) const
const G4InuclElementaryParticle & getParticle() const
G4double getCurrentDensity(G4int ip, G4int izone) const
G4GLOB_DLL std::ostream G4cout
G4double getKinEnergyInTheTRS() const
G4double totalCrossSection(G4double ke, G4int rtype) const
G4int getCurrentZone() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void setTarget(const G4InuclParticle *target)

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.

1308  {
1309  return (cparticle.getGeneration() == 0); // Only initial state particles
1310 }
G4int getGeneration() const

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.

1062  {
1063  if (verboseLevel > 1) {
1064  G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
1065  }
1066 
1067  // Only check Fermi momenta for nucleons
1068  for (G4int i = 0; i < G4int(particles.size()); i++) {
1069  if (!particles[i].nucleon()) continue;
1070 
1071  G4int type = particles[i].type();
1072  G4double mom = particles[i].getMomModule();
1073  G4double pfermi = fermi_momenta[type-1][zone];
1074 
1075  if (verboseLevel > 2)
1076  G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
1077 
1078  if (mom < pfermi) {
1079  if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
1080  return false;
1081  }
1082  }
1083  return true;
1084 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4bool nucleon(G4int ityp)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

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.

1090  {
1091  if (verboseLevel > 1)
1092  G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
1093 
1094  G4double dist;
1095  for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
1096  dist = (collisionPts[i] - hit_position).mag();
1097  if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
1098  if (dist < R_nucleon) {
1099  if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
1100  return false;
1101  }
1102  }
1103  return true; // New point far enough away to be used
1104 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

void G4NucleiModel::printModel ( ) const

Definition at line 604 of file G4NucleiModel.cc.

604  {
605  if (verboseLevel > 1) {
606  G4cout << " >>> G4NucleiModel::printModel" << G4endl;
607  }
608 
609  G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
610  << " proton binding energy " << binding_energies[0]
611  << " neutron binding energy " << binding_energies[1] << G4endl
612  << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
613  << " number of zones " << number_of_zones << G4endl;
614 
615  for (G4int i = 0; i < number_of_zones; i++)
616  G4cout << " zone " << i+1 << " radius " << zone_radii[i]
617  << " volume " << zone_volumes[i] << G4endl
618  << " protons: density " << getDensity(1,i) << " PF " <<
619  getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
620  << " neutrons: density " << getDensity(2,i) << " PF " <<
621  getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
622  << " pions: VP " << getPotential(3,i) << G4endl;
623 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double getFermiMomentum(G4int ip, G4int izone) const
#define G4endl
Definition: G4ios.hh:61
G4double getDensity(G4int ip, G4int izone) const
G4double getPotential(G4int ip, G4int izone) const

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.

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

Here is the caller graph for this function:

void G4NucleiModel::setVerboseLevel ( G4int  verbose)
inline

Definition at line 101 of file G4NucleiModel.hh.

101 { verboseLevel = verbose; }
static G4bool G4NucleiModel::sortPartners ( const partner p1,
const partner p2 
)
inlinestaticprotected

Definition at line 211 of file G4NucleiModel.hh.

211  {
212  return (p2.second > p1.second);
213  }

Here is the caller graph for this function:

G4bool G4NucleiModel::stillInside ( const G4CascadParticle cparticle)
inline

Definition at line 156 of file G4NucleiModel.hh.

156  {
157  return cparticle.getCurrentZone() < number_of_zones;
158  }
G4int getCurrentZone() const

Here is the call graph for this function:

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

Definition at line 1918 of file G4NucleiModel.cc.

1919 {
1920  // All scattering cross-sections are available from tables
1921  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
1922  if (!xsecTable) {
1923  G4cerr << " unknown collison type = " << rtype << G4endl;
1924  return 0.;
1925  }
1926 
1927  return (crossSectionUnits * xsecTable->getCrossSection(ke));
1928 }
static const G4CascadeChannel * GetTable(G4int initialState)
virtual G4double getCrossSection(double ke) const =0
#define G4endl
Definition: G4ios.hh:61
G4GLOB_DLL std::ostream G4cerr

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.

1049  {
1050  if (qdtype==pn || qdtype==0) // All absorptive particles
1051  return (ptype==pi0 || ptype==pip || ptype==pim || ptype==gam || ptype==mum);
1052  else if (qdtype==pp) // Negative or neutral only
1053  return (ptype==pi0 || ptype==pim || ptype==gam || ptype==mum);
1054  else if (qdtype==nn) // Positive or neutral only
1055  return (ptype==pi0 || ptype==pip || ptype==gam);
1056 
1057  return false; // Not a quasideuteron target
1058 }

Here is the caller graph for this function:

G4bool G4NucleiModel::worthToPropagate ( const G4CascadParticle cparticle) const

Definition at line 1312 of file G4NucleiModel.cc.

1312  {
1313  if (verboseLevel > 1) {
1314  G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
1315  }
1316 
1317  const G4double ekin_scale = 2.0;
1318 
1319  G4bool worth = true;
1320 
1321  if (cparticle.reflectedNow()) { // Just reflected -- keep going?
1322  G4int zone = cparticle.getCurrentZone();
1323  G4int ip = cparticle.getParticle().type();
1324 
1325  // NOTE: Temporarily backing out use of potential for non-nucleons
1326  G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
1327  getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
1328 
1329  worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
1330 
1331  if (verboseLevel > 3) {
1332  G4cout << " type=" << ip
1333  << " ekin=" << cparticle.getParticle().getKineticEnergy()
1334  << " potential=" << ekin_cut
1335  << " : worth? " << worth << G4endl;
1336  }
1337  }
1338 
1339  return worth;
1340 }
G4bool reflectedNow() const
G4double getFermiKinetic(G4int ip, G4int izone) const
int G4int
Definition: G4Types.hh:78
const G4InuclElementaryParticle & getParticle() const
G4double getKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4int getCurrentZone() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

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.

557  {
558  if (verboseLevel > 1) {
559  G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
560  }
561 
562  G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
563 
564  const G4double epsilon = 1.0e-3;
565  const G4int itry_max = 1000;
566 
567  G4double dr = r2 - r1;
568  G4double fr1 = r1 * r1 * G4Exp(-r1 * r1);
569  G4double fr2 = r2 * r2 * G4Exp(-r2 * r2);
570  G4double fi = (fr1 + fr2) / 2.;
571  G4double fun1 = fi * dr;
572  G4double fun;
573  G4int jc = 1;
574  G4double dr1 = dr;
575  G4int itry = 0;
576 
577  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
578  dr /= 2.;
579  itry++;
580  G4double r = r1 - dr;
581  fi = 0.0;
582 
583  for (G4int i = 0; i < jc; i++) {
584  r += dr1;
585  fi += r * r * G4Exp(-r * r);
586  }
587 
588  fun = 0.5 * fun1 + fi * dr;
589 
590  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
591 
592  jc *= 2;
593  dr1 = dr;
594  fun1 = fun;
595  } // while (itry < itry_max)
596 
597  if (verboseLevel > 2 && itry == itry_max)
598  G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
599 
600  return gaussRadius*gaussRadius*gaussRadius * fun;
601 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
double epsilon(double density, double temperature)
G4GLOB_DLL std::ostream G4cerr

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.

504  {
505  if (verboseLevel > 1) {
506  G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
507  }
508 
509  const G4double epsilon = 1.0e-3;
510  const G4int itry_max = 1000;
511 
512  G4double skinRatio = nuclearRadius / skinDepth;
513 
514  G4double d2 = 2.0 * skinRatio;
515  G4double dr = r2 - r1;
516  G4double fr1 = r1 * (r1 + d2) / (1.0 + G4Exp(r1));
517  G4double fr2 = r2 * (r2 + d2) / (1.0 + G4Exp(r2));
518  G4double fi = (fr1 + fr2) / 2.;
519  G4double fun1 = fi * dr;
520  G4double fun;
521  G4int jc = 1;
522  G4double dr1 = dr;
523  G4int itry = 0;
524 
525  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
526  dr /= 2.;
527  itry++;
528 
529  G4double r = r1 - dr;
530  fi = 0.0;
531 
532  for (G4int i = 0; i < jc; i++) {
533  r += dr1;
534  fi += r * (r + d2) / (1.0 + G4Exp(r));
535  }
536 
537  fun = 0.5 * fun1 + fi * dr;
538 
539  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
540 
541  jc *= 2;
542  dr1 = dr;
543  fun1 = fun;
544  } // while (itry < itry_max)
545 
546  if (verboseLevel > 2 && itry == itry_max)
547  G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
548 
549  G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
550 
551  return skinDepth3 * (fun + skinRatio*skinRatio*G4Log((1.0 + G4Exp(-r1)) / (1.0 + G4Exp(-r2))));
552 }
static const G4double d2
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
double epsilon(double density, double temperature)

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: