Geant4_10
G4NucleiModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4NucleiModel.cc 71989 2013-07-02 17:12:22Z mkelsey $
27 //
28 // 20100112 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100114 M. Kelsey -- Use G4ThreeVector for position
30 // 20100309 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
31 // eliminate some unnecessary std::pow(), move ctor here
32 // 20100407 M. Kelsey -- Replace std::vector<>::resize(0) with ::clear().
33 // Move output vectors from generateParticleFate() and
34 // ::generateInteractionPartners() to be data members; return via
35 // const-ref instead of by value.
36 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
37 // 20100418 M. Kelsey -- Reference output particle lists via const-ref
38 // 20100421 M. Kelsey -- Replace hardwired p/n masses with G4PartDef's
39 // 20100517 M. Kelsey -- Use G4CascadeINterpolator for cross-section
40 // calculations. use G4Cascade*Channel for total xsec in pi-N
41 // N-N channels. Move absorptionCrossSection() from SpecialFunc.
42 // 20100610 M. Kelsey -- Replace another random-angle code block; add some
43 // diagnostic output for partner-list production.
44 // 20100617 M. Kelsey -- Replace preprocessor flag CHC_CHECK with
45 // G4CASCADE_DEBUG_CHARGE
46 // 20100620 M. Kelsey -- Improve error message on empty partners list, add
47 // four-momentum checking after EPCollider
48 // 20100621 M. Kelsey -- In boundaryTransition() account for momentum transfer
49 // to secondary by boosting into recoil nucleus "rest" frame.
50 // Don't need temporaries to create dummy "partners" for list.
51 // 20100622 M. Kelsey -- Restore use of "bindingEnergy()" function name, which
52 // is now a wrapper for G4NucleiProperties::GetBindingEnergy().
53 // 20100623 M. Kelsey -- Eliminate some temporaries terminating partner-list,
54 // discard recoil boost for now. Initialize all data
55 // members in ctors. Allow generateModel() to be called
56 // mutliple times, by clearing vectors each time through;
57 // avoid extra work by returning if A and Z are same as
58 // before.
59 // 20100628 M. Kelsey -- Two momentum-recoil bugs; don't subtract energies!
60 // 20100715 M. Kelsey -- Make G4InuclNuclei in generateModel(), use for
61 // balance checking (only verbose>2) in generateParticleFate().
62 // 20100721 M. Kelsey -- Use new G4CASCADE_CHECK_ECONS for conservation checks
63 // 20100723 M. Kelsey -- Move G4CollisionOutput buffer to .hh for reuse
64 // 20100726 M. Kelsey -- Preallocate arrays with number_of_zones dimension.
65 // 20100902 M. Kelsey -- Remove resize(3) directives from qdeutron/acsecs
66 // 20100907 M. Kelsey -- Limit interaction targets based on current nucleon
67 // counts (e.g., no pp if protonNumberCurrent < 2!)
68 // 20100914 M. Kelsey -- Migrate to integer A and Z
69 // 20100928 M. Kelsey -- worthToPropagate() use nuclear potential for hadrons.
70 // 20101005 M. Kelsey -- Annotate hardwired numbers for upcoming rationalizing,
71 // move hardwired numbers out to static data members, factor out
72 // all the model construction pieces to separate functions, fix
73 // wrong value for 4/3 pi.
74 // 20101019 M. Kelsey -- CoVerity reports: unitialized constructor, dtor leak
75 // 20101020 M. Kelsey -- Bug fixes to refactoring changes (5 Oct). Back out
76 // worthToPropagate() changes for better regression testing.
77 // 20101020 M. Kelsey -- Re-activate worthToPropagate() changes.
78 // 20101119 M. Kelsey -- Hide "negative path" and "no partners" messages in
79 // verbosity.
80 // 20110218 M. Kelsey -- Add crossSectionUnits and radiusUnits scale factors,
81 // use "theoretical" numbers for radii etc., multipled by scale
82 // factor; set scale factors using environment variables
83 // 20110303 M. Kelsey -- Add comments why using fabs() with B.E. differences?
84 // 20110321 M. Kelsey -- Replace strtof() with strtod() for envvar conversion
85 // 20110321 M. Kelsey -- Use fm and fm^2 as default units, Per D. Wright
86 // (NOTE: Restored from original 20110318 commit)
87 // 20110324 D. Wright -- Implement trailing effect
88 // 20110324 M. Kelsey -- Move ::reset() here, as it has more code.
89 // 20110519 M. Kelsey -- Used "rho" after assignment, instead of recomputing
90 // 20110525 M. Kelsey -- Revert scale factor changes (undo 20110321 changes)
91 // 20110617 M. Kelsey -- Apply scale factor to trailing-effect radius, make
92 // latter runtime adjustable (G4NUCMODEL_RAD_TRAILING)
93 // 20110720 M. Kelsey -- Follow interface change for cross-section tables,
94 // eliminating switch blocks.
95 // 20110806 M. Kelsey -- Reduce memory churn by pre-allocating buffers
96 // 20110823 M. Kelsey -- Remove local cross-section tables entirely
97 // 20110825 M. Kelsey -- Add comments regarding Fermi momentum scale, set of
98 // "best guess" parameter values
99 // 20110831 M. Kelsey -- Make "best guess" parameters the defaults
100 // 20110922 M. Kelsey -- Follow migrations G4InuclParticle::print(ostream&)
101 // and G4CascadParticle::print(ostream&)
102 // 20111018 M. Kelsey -- Correct kaon potential to be positive, not negative
103 // 20111107 M. Kelsey -- *** REVERT TO OLD NON-PHYSICAL PARAMETERS FOR 9.5 ***
104 // 20120306 M. Kelsey -- For incident projectile (CParticle incoming to
105 // nucleus from outside) don't use pw cut; force interaction.
106 // 20120307 M. Kelsey -- Compute zone volumes during setup, not during each
107 // interaction. Include photons in absorption by dibaryons.
108 // Discard unused code in totalCrossSection(). Encapsulate
109 // interaction path calculations in generateInteractionLength().
110 // 20120308 M. Kelsey -- Add binned photon absorption cross-section.
111 // 20120320 M. Kelsey -- Bug fix: fill zone_volumes for single-zone case
112 // 20120321 M. Kelsey -- Add check in isProjectile() for single-zone case.
113 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
114 // 20120822 M. Kelsey -- Move envvars to G4CascadeParameters.
115 // 20121205 M. Kelsey -- Bug fix in generateParticleFate(): daughters should
116 // have generation count incremented from parent. This allows
117 // isProjectile() to simply check generation number == 0.
118 // 20130221 M. Kelsey -- For incident photons, move to random point along
119 // trajectory through nucleus for first (forced) interaction.
120 // 20130223 M. Kelsey -- Fix bugs in weighting and selection of traj points
121 // 20130302 M. Kelsey -- Replace "isProjectile()" with "forceFirst()" in
122 // path-length calculation.
123 // 20130314 M. Kelsey -- Restore null initializer and if-block for _TLS_.
124 // 20130508 D. Wright -- Support muon-nuclear interactions.
125 // 20130510 M. Kelsey -- Check for zero-momentum in choosePointAlongTraj().
126 // 20130511 M. Kelsey -- Move choosePointAlongTraj() to initializeCascad(),
127 // instead of after generateIP(). Change "spath<path" to <= to
128 // handle at-rest particles. Skip mu-/neutron interactions.
129 // Hide "zone 0 transition" error message behind verbosity.
130 // 20130523 M. Kelsey -- Bug fix: Initialize dummy_convertor in inverseMFP.
131 // For capture-at-rest, replace exterior nucleus with outer zone.
132 // 20130524 M. Kelsey -- Move "large" and "small" cutoffs to static const.
133 // Check zone argument to inverseMFP() to ensure within range.
134 // 20130611 M. Kelsey -- Undo "spath<=path" change (20130511), replace with
135 // explicit check on spath==0.
136 // 20130619 A. Ribon -- Fixed reproducibility problem in the method
137 // generateParticleFate
138 // 20130627 M. Kelsey -- Check "path==0.", rather than spath.
139 // 20130628 M. Kelsey -- Print deuteron type code, rather than array index,
140 // Extend useQuasiDeuteron() to check good absorption
141 // 20130701 M. Kelsey -- Don't average 1/MFP for total interaction; just sum;
142 // inverseMFP() returns "large" value instead of input path.
143 // 20130806 M. Kelsey -- Move static G4InuclEP's to file scope to avoid
144 // thread collisions
145 // 20130924 M. Kelsey -- Use G4Log, G4Exp, G4Pow for CPU speedup
146 // 20130925 M. Kelsey -- Eliminate unnecessary use of pow() in integrals
147 // 20131001 M. Kelsey -- Move QDinterp object to data member, thread local
148 
149 #include "G4NucleiModel.hh"
150 #include "G4PhysicalConstants.hh"
151 #include "G4SystemOfUnits.hh"
152 #include "G4CascadeChannel.hh"
153 #include "G4CascadeChannelTables.hh"
154 #include "G4CascadeCheckBalance.hh"
155 #include "G4CascadeInterpolator.hh"
156 #include "G4CascadeParameters.hh"
157 #include "G4CollisionOutput.hh"
159 #include "G4HadTmpUtil.hh"
160 #include "G4InuclNuclei.hh"
161 #include "G4InuclParticleNames.hh"
163 #include "G4LorentzConvertor.hh"
164 #include "G4Neutron.hh"
165 #include "G4ParticleLargerBeta.hh"
166 #include "G4ParticleDefinition.hh"
167 #include "G4Proton.hh"
168 #include "G4Log.hh"
169 #include "G4Exp.hh"
170 #include <vector>
171 #include <algorithm>
172 #include <functional>
173 #include <numeric>
174 
175 using namespace G4InuclParticleNames;
176 using namespace G4InuclSpecialFunctions;
177 
178 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
179 
180 // For the best approximation to a physical-units model, set the following:
181 // setenv G4NUCMODEL_XSEC_SCALE 0.1
182 // setenv G4NUCMODEL_RAD_SCALE 1.0
183 // setenv G4NUCMODEL_RAD_2PAR 1
184 // setenv G4NUCMODEL_RAD_SMALL 1.992
185 // setenv G4NUCMODEL_RAD_ALPHA 0.84
186 // setenv G4NUCMODEL_FERMI_SCALE 0.685
187 // setenv G4NUCMODEL_RAD_TRAILING 0.70
188 
189 // Scaling factors for radii and cross-sections, currently different!
190 const G4double G4NucleiModel::crossSectionUnits = G4CascadeParameters::xsecScale();
191 const G4double G4NucleiModel::radiusUnits = G4CascadeParameters::radiusScale();
192 const G4double G4NucleiModel::skinDepth = 0.611207*radiusUnits;
193 
194 // One- vs. two-parameter nuclear radius based on envvar
195 // ==> radius = radiusScale*cbrt(A) + radiusScale2/cbrt(A)
196 
197 const G4double G4NucleiModel::radiusScale =
198  (G4CascadeParameters::useTwoParam() ? 1.16 : 1.2) * radiusUnits;
199 const G4double G4NucleiModel::radiusScale2 =
200  (G4CascadeParameters::useTwoParam() ? -1.3456 : 0.) * radiusUnits;
201 
202 // NOTE: Old code used R_small = 8.0 (~2.83*units), and R_alpha = 0.7*R_small
203 // Published data suggests R_small ~ 1.992 fm, R_alpha = 0.84*R_small
204 
205 const G4double G4NucleiModel::radiusForSmall = G4CascadeParameters::radiusSmall();
206 
207 const G4double G4NucleiModel::radScaleAlpha = G4CascadeParameters::radiusAlpha();
208 
209 // Scale factor relating Fermi momentum to density of states, units GeV.fm
210 // NOTE: Old code has 0.685*units GeV.fm, literature suggests 0.470 GeV.fm,
211 // but this gives too small momentum; old value gives <P_F> ~ 270 MeV
212 const G4double G4NucleiModel::fermiMomentum = G4CascadeParameters::fermiScale();
213 
214 // Effective radius (0.87 to 1.2 fm) of nucleon, for trailing effect
215 const G4double G4NucleiModel::R_nucleon = G4CascadeParameters::radiusTrailing();
216 
217 // Zone boundaries as fraction of nuclear radius (from outside in)
218 const G4double G4NucleiModel::alfa3[3] = { 0.7, 0.3, 0.01 };
219 const G4double G4NucleiModel::alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
220 
221 // Flat nuclear potentials for mesons and hyperons (GeV)
222 const G4double G4NucleiModel::pion_vp = 0.007;
223 const G4double G4NucleiModel::pion_vp_small = 0.007;
224 const G4double G4NucleiModel::kaon_vp = 0.015;
225 const G4double G4NucleiModel::hyperon_vp = 0.030;
226 
227 // Cross-section scaling parameter for gamma-quasideuteron scattering
228 const G4double G4NucleiModel::gammaQDscale = G4CascadeParameters::gammaQDScale();
229 
230 // For convenience in computing densities
231 const G4double G4NucleiModel::piTimes4thirds = pi*4./3.;
232 
233 // Cut-off limits for extreme values in calculations
234 const G4double G4NucleiModel::small = 1.0e-9;
235 const G4double G4NucleiModel::large = 1000.;
236 
237 namespace {
238  const G4double kebins[] =
239  { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
240  0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
241  2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
242 
243  const G4double gammaQDxsec[30] =
244  { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098,
245  0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005,
246  0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
247  0.0001, 0.0001, 0.0001 };
248 }
249 
250 
251 // Constructors and destructor
253  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
254  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
255  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
256  current_nucl2(0), gammaQDinterp(kebins),
257  neutronEP(neutron), protonEP(proton) {}
258 
260  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
261  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
262  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
263  current_nucl2(0), gammaQDinterp(kebins),
264  neutronEP(neutron), protonEP(proton) {
265  generateModel(a,z);
266 }
267 
269  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
270  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
271  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
272  current_nucl2(0), gammaQDinterp(kebins),
273  neutronEP(neutron), protonEP(proton) {
274  generateModel(nuclei);
275 }
276 
278  delete theNucleus;
279  theNucleus = 0;
280 }
281 
282 
283 // Initialize model state for new cascade
284 
285 void G4NucleiModel::reset(G4int nHitNeutrons, G4int nHitProtons,
286  const std::vector<G4ThreeVector>* hitPoints) {
287  neutronNumberCurrent = neutronNumber - nHitNeutrons;
288  protonNumberCurrent = protonNumber - nHitProtons;
289 
290  // zero or copy collision point array for trailing effect
291  if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
292  else collisionPts = *hitPoints;
293 }
294 
295 
296 // Generate nuclear model parameters for given nucleus
297 
299  generateModel(nuclei->getA(), nuclei->getZ());
300 }
301 
303  if (verboseLevel) {
304  G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
305  << G4endl;
306  }
307 
308  // If model already built, just return; otherwise intialize everything
309  if (a == A && z == Z) {
310  if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
311  reset();
312  return;
313  }
314 
315  A = a;
316  Z = z;
317  delete theNucleus;
318  theNucleus = new G4InuclNuclei(A,Z); // For conservation checking
319 
320  neutronNumber = A - Z;
321  protonNumber = Z;
322  reset();
323 
324  if (verboseLevel > 3) {
325  G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl
326  << " radiusUnits = " << radiusUnits << G4endl
327  << " skinDepth = " << skinDepth << G4endl
328  << " radiusScale = " << radiusScale << G4endl
329  << " radiusScale2 = " << radiusScale2 << G4endl
330  << " radiusForSmall = " << radiusForSmall << G4endl
331  << " radScaleAlpha = " << radScaleAlpha << G4endl
332  << " fermiMomentum = " << fermiMomentum << G4endl
333  << " piTimes4thirds = " << piTimes4thirds << G4endl;
334  }
335 
336  G4double nuclearRadius; // Nuclear radius computed from A
337  if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
338  else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
339 
340  // This will be used to pre-allocate lots of arrays below
341  number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
342 
343  // Clear all parameters arrays for reloading
344  binding_energies.clear();
345  nucleon_densities.clear();
346  zone_potentials.clear();
347  fermi_momenta.clear();
348  zone_radii.clear();
349  zone_volumes.clear();
350 
352  fillZoneRadii(nuclearRadius);
353 
354  G4double tot_vol = fillZoneVolumes(nuclearRadius); // Woods-Saxon integral
355 
356  fillPotentials(proton, tot_vol); // Protons
357  fillPotentials(neutron, tot_vol); // Neutrons
358 
359  // Additional flat zone potentials for other hadrons
360  const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
361  const std::vector<G4double> kp(number_of_zones, kaon_vp);
362  const std::vector<G4double> hp(number_of_zones, hyperon_vp);
363 
364  zone_potentials.push_back(vp);
365  zone_potentials.push_back(kp);
366  zone_potentials.push_back(hp);
367 
368  nuclei_radius = zone_radii.back();
369  nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
370 
371  if (verboseLevel > 3) printModel();
372 }
373 
374 
375 // Load binding energy array for current nucleus
376 
378  if (verboseLevel > 1)
379  G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
380 
381  G4double dm = bindingEnergy(A,Z);
382 
383  // Binding energy differences for proton and neutron loss, respectively
384  // FIXME: Why is fabs() used here instead of the signed difference?
385  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
386  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
387 }
388 
389 // Load zone boundary radius array for current nucleus
390 
392  if (verboseLevel > 1)
393  G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
394 
395  G4double skinRatio = nuclearRadius/skinDepth;
396  G4double skinDecay = G4Exp(-skinRatio);
397 
398  if (A < 5) { // Light ions treated as simple balls
399  zone_radii.push_back(nuclearRadius);
400  ur[0] = 0.;
401  ur[1] = 1.;
402  } else if (A < 12) { // Small nuclei have Gaussian potential
403  G4double rSq = nuclearRadius * nuclearRadius;
404  G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
405 
406  ur[0] = 0.0;
407  for (G4int i = 0; i < number_of_zones; i++) {
408  G4double y = std::sqrt(-G4Log(alfa3[i]));
409  zone_radii.push_back(gaussRadius * y);
410  ur[i+1] = y;
411  }
412  } else if (A < 100) { // Intermediate nuclei have three-zone W-S
413  ur[0] = -skinRatio;
414  for (G4int i = 0; i < number_of_zones; i++) {
415  G4double y = G4Log((1.0 + skinDecay)/alfa3[i] - 1.0);
416  zone_radii.push_back(nuclearRadius + skinDepth * y);
417  ur[i+1] = y;
418  }
419  } else { // Heavy nuclei have six-zone W-S
420  ur[0] = -skinRatio;
421  for (G4int i = 0; i < number_of_zones; i++) {
422  G4double y = G4Log((1.0 + skinDecay)/alfa6[i] - 1.0);
423  zone_radii.push_back(nuclearRadius + skinDepth * y);
424  ur[i+1] = y;
425  }
426  }
427 }
428 
429 // Compute zone-by-zone density integrals
430 
432  if (verboseLevel > 1)
433  G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
434 
435  G4double tot_vol = 0.; // Return value omitting 4pi/3 for sphere!
436 
437  if (A < 5) { // Light ions are treated as simple balls
438  v[0] = v1[0] = 1.;
439  tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
440  zone_volumes.push_back(tot_vol*piTimes4thirds); // True volume of sphere
441 
442  return tot_vol;
443  }
444 
445  PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
446 
447  for (G4int i = 0; i < number_of_zones; i++) {
448  if (usePotential == WoodsSaxon) {
449  v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
450  } else {
451  v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
452  }
453 
454  tot_vol += v[i];
455 
456  v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
457  if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
458 
459  zone_volumes.push_back(v1[i]*piTimes4thirds); // True volume of shell
460  }
461 
462  return tot_vol; // Sum of zone integrals, not geometric volume
463 }
464 
465 // Load potentials for nucleons (using G4InuclParticle codes)
467  if (verboseLevel > 1)
468  G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
469 
470  if (type != proton && type != neutron) return;
471 
473 
474  // FIXME: This is the fabs() binding energy difference, not signed
475  const G4double dm = binding_energies[type-1];
476 
477  rod.clear(); rod.reserve(number_of_zones);
478  pf.clear(); pf.reserve(number_of_zones);
479  vz.clear(); vz.reserve(number_of_zones);
480 
481  G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
482  G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
483 
484  for (G4int i = 0; i < number_of_zones; i++) {
485  G4double rd = dd0 * v[i] / v1[i];
486  rod.push_back(rd);
487  G4double pff = fermiMomentum * G4cbrt(rd);
488  pf.push_back(pff);
489  vz.push_back(0.5 * pff * pff / mass + dm);
490  }
491 
492  nucleon_densities.push_back(rod);
493  fermi_momenta.push_back(pf);
494  zone_potentials.push_back(vz);
495 }
496 
497 // Zone integral of Woods-Saxon density function
499  G4double nuclearRadius) const {
500  if (verboseLevel > 1) {
501  G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
502  }
503 
504  const G4double epsilon = 1.0e-3;
505  const G4int itry_max = 1000;
506 
507  G4double skinRatio = nuclearRadius / skinDepth;
508 
509  G4double d2 = 2.0 * skinRatio;
510  G4double dr = r2 - r1;
511  G4double fr1 = r1 * (r1 + d2) / (1.0 + G4Exp(r1));
512  G4double fr2 = r2 * (r2 + d2) / (1.0 + G4Exp(r2));
513  G4double fi = (fr1 + fr2) / 2.;
514  G4double fun1 = fi * dr;
515  G4double fun;
516  G4int jc = 1;
517  G4double dr1 = dr;
518  G4int itry = 0;
519 
520  while (itry < itry_max) {
521  dr /= 2.;
522  itry++;
523 
524  G4double r = r1 - dr;
525  fi = 0.0;
526 
527  for (G4int i = 0; i < jc; i++) {
528  r += dr1;
529  fi += r * (r + d2) / (1.0 + G4Exp(r));
530  }
531 
532  fun = 0.5 * fun1 + fi * dr;
533 
534  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
535 
536  jc *= 2;
537  dr1 = dr;
538  fun1 = fun;
539  } // while (itry < itry_max)
540 
541  if (verboseLevel > 2 && itry == itry_max)
542  G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
543 
544  G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
545 
546  return skinDepth3 * (fun + skinRatio*skinRatio*G4Log((1.0 + G4Exp(-r1)) / (1.0 + G4Exp(-r2))));
547 }
548 
549 
550 // Zone integral of Gaussian density function
552  G4double nucRad) const {
553  if (verboseLevel > 1) {
554  G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
555  }
556 
557  G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
558 
559  const G4double epsilon = 1.0e-3;
560  const G4int itry_max = 1000;
561 
562  G4double dr = r2 - r1;
563  G4double fr1 = r1 * r1 * G4Exp(-r1 * r1);
564  G4double fr2 = r2 * r2 * G4Exp(-r2 * r2);
565  G4double fi = (fr1 + fr2) / 2.;
566  G4double fun1 = fi * dr;
567  G4double fun;
568  G4int jc = 1;
569  G4double dr1 = dr;
570  G4int itry = 0;
571 
572  while (itry < itry_max) {
573  dr /= 2.;
574  itry++;
575  G4double r = r1 - dr;
576  fi = 0.0;
577 
578  for (G4int i = 0; i < jc; i++) {
579  r += dr1;
580  fi += r * r * G4Exp(-r * r);
581  }
582 
583  fun = 0.5 * fun1 + fi * dr;
584 
585  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
586 
587  jc *= 2;
588  dr1 = dr;
589  fun1 = fun;
590  } // while (itry < itry_max)
591 
592  if (verboseLevel > 2 && itry == itry_max)
593  G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
594 
595  return gaussRadius*gaussRadius*gaussRadius * fun;
596 }
597 
598 
600  if (verboseLevel > 1) {
601  G4cout << " >>> G4NucleiModel::printModel" << G4endl;
602  }
603 
604  G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
605  << " proton binding energy " << binding_energies[0]
606  << " neutron binding energy " << binding_energies[1] << G4endl
607  << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
608  << " number of zones " << number_of_zones << G4endl;
609 
610  for (G4int i = 0; i < number_of_zones; i++)
611  G4cout << " zone " << i+1 << " radius " << zone_radii[i]
612  << " volume " << zone_volumes[i] << G4endl
613  << " protons: density " << getDensity(1,i) << " PF " <<
614  getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
615  << " neutrons: density " << getDensity(2,i) << " PF " <<
616  getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
617  << " pions: VP " << getPotential(3,i) << G4endl;
618 }
619 
620 
622  G4double ekin = 0.0;
623 
624  if (ip < 3 && izone < number_of_zones) { // ip for proton/neutron only
625  G4double pfermi = fermi_momenta[ip - 1][izone];
627 
628  ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
629  }
630  return ekin;
631 }
632 
633 
636  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
638 
639  return generateWithRandomAngles(pmod, mass);
640 }
641 
642 
645  if (verboseLevel > 1) {
646  G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
647  }
648 
649  G4LorentzVector mom = generateNucleonMomentum(type, zone);
650  return G4InuclElementaryParticle(mom, type);
651 }
652 
653 
656  G4int zone) const {
657  if (verboseLevel > 1) {
658  G4cout << " >>> G4NucleiModel::generateQuasiDeuteron" << G4endl;
659  }
660 
661  // Quasideuteron consists of an unbound but associated nucleon pair
662 
663  // FIXME: Why generate two separate nucleon momenta (randomly!) and
664  // add them, instead of just throwing a net momentum for the
665  // dinulceon state? And why do I have to capture the two
666  // return values into local variables?
667  G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
668  G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
669  G4LorentzVector dmom = mom1+mom2;
670 
671  G4int dtype = 0;
672  if (type1*type2 == pro*pro) dtype = 111;
673  else if (type1*type2 == pro*neu) dtype = 112;
674  else if (type1*type2 == neu*neu) dtype = 122;
675 
676  return G4InuclElementaryParticle(dmom, dtype);
677 }
678 
679 
680 void
682  if (verboseLevel > 1) {
683  G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
684  }
685 
686  thePartners.clear(); // Reset buffer for next cycle
687 
688  G4int ptype = cparticle.getParticle().type();
689  G4int zone = cparticle.getCurrentZone();
690 
691  G4double r_in; // Destination radius if inbound
692  G4double r_out; // Destination radius if outbound
693 
694  if (zone == number_of_zones) {
695  r_in = nuclei_radius;
696  r_out = 0.0;
697  } else if (zone == 0) { // particle is outside core
698  r_in = 0.0;
699  r_out = zone_radii[0];
700  } else {
701  r_in = zone_radii[zone - 1];
702  r_out = zone_radii[zone];
703  }
704 
705  G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
706 
707  if (verboseLevel > 2) {
708  if (isProjectile(cparticle)) G4cout << " incident particle: ";
709  G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
710  << G4endl;
711  }
712 
713  if (path < -small) { // something wrong
714  if (verboseLevel)
715  G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
716  return;
717  }
718 
719  if (std::fabs(path) < small) { // Not moving, or just at boundary
720  if (cparticle.getMomentum().vect().mag() > small) {
721  if (verboseLevel > 3)
722  G4cout << " generateInteractionPartners-> zero path" << G4endl;
723 
724  thePartners.push_back(partner()); // Dummy list terminator with zero path
725  return;
726  }
727 
728  if (zone >= number_of_zones) // Place captured-at-rest in nucleus
729  zone = number_of_zones-1;
730  }
731 
732  G4double invmfp = 0.; // Buffers for interaction probability
733  G4double spath = 0.;
734  for (G4int ip = 1; ip < 3; ip++) {
735  // Only process nucleons which remain active in target
736  if (ip==proton && protonNumberCurrent < 1) continue;
737  if (ip==neutron && neutronNumberCurrent < 1) continue;
738  if (ip==neutron && ptype==muonMinus) continue; // mu-/n forbidden
739 
740  // All nucleons are assumed to be at rest when colliding
741  G4InuclElementaryParticle particle = generateNucleon(ip, zone);
742  invmfp = inverseMeanFreePath(cparticle, particle);
743  spath = generateInteractionLength(cparticle, path, invmfp);
744 
745  if (path<small || spath < path) {
746  if (verboseLevel > 3) {
747  G4cout << " adding partner[" << thePartners.size() << "]: "
748  << particle << G4endl;
749  }
750  thePartners.push_back(partner(particle, spath));
751  }
752  } // for (G4int ip...
753 
754  if (verboseLevel > 2) {
755  G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
756  }
757 
758  // Absorption possible for pions or photons interacting with dibaryons
759  if (useQuasiDeuteron(cparticle.getParticle().type())) {
760  if (verboseLevel > 2) {
761  G4cout << " trying quasi-deuterons with bullet: "
762  << cparticle.getParticle() << G4endl;
763  }
764 
765  // Initialize buffers for quasi-deuteron results
766  qdeutrons.clear();
767  acsecs.clear();
768 
769  G4double tot_invmfp = 0.0; // Total inv. mean-free-path for all QDs
770 
771  // Proton-proton state interacts with pi-, mu- or neutrals
772  if (protonNumberCurrent >= 2 && ptype != pip) {
774  if (verboseLevel > 2)
775  G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
776 
777  invmfp = inverseMeanFreePath(cparticle, ppd);
778  tot_invmfp += invmfp;
779  acsecs.push_back(invmfp);
780  qdeutrons.push_back(ppd);
781  }
782 
783  // Proton-neutron state interacts with any pion type or photon
784  if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
786  if (verboseLevel > 2)
787  G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
788 
789  invmfp = inverseMeanFreePath(cparticle, npd);
790  tot_invmfp += invmfp;
791  acsecs.push_back(invmfp);
792  qdeutrons.push_back(npd);
793  }
794 
795  // Neutron-neutron state interacts with pi+ or neutrals
796  if (neutronNumberCurrent >= 2 && ptype != pim && ptype != mum) {
798  if (verboseLevel > 2)
799  G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
800 
801  invmfp = inverseMeanFreePath(cparticle, nnd);
802  tot_invmfp += invmfp;
803  acsecs.push_back(invmfp);
804  qdeutrons.push_back(nnd);
805  }
806 
807  // Select quasideuteron interaction from non-zero cross-section choices
808  if (verboseLevel > 2) {
809  for (size_t i=0; i<qdeutrons.size(); i++) {
810  G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
811  << "] " << acsecs[i];
812  }
813  G4cout << G4endl;
814  }
815 
816  if (tot_invmfp > small) { // Must have absorption cross-section
817  G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
818 
819  if (path<small || apath < path) { // choose the qdeutron
820  G4double sl = inuclRndm() * tot_invmfp;
821  G4double as = 0.0;
822 
823  for (size_t i = 0; i < qdeutrons.size(); i++) {
824  as += acsecs[i];
825  if (sl < as) {
826  if (verboseLevel > 2)
827  G4cout << " deut type " << qdeutrons[i] << G4endl;
828 
829  thePartners.push_back(partner(qdeutrons[i], apath));
830  break;
831  }
832  } // for (qdeutrons...
833  } // if (apath < path)
834  } // if (tot_invmfp > small)
835  } // if (useQuasiDeuteron) [pion, muon or photon]
836 
837  if (verboseLevel > 2) {
838  G4cout << " after deuterons " << thePartners.size() << " partners"
839  << G4endl;
840  }
841 
842  if (thePartners.size() > 1) { // Sort list by path length
844  }
845 
846  G4InuclElementaryParticle particle; // Total path at end of list
847  thePartners.push_back(partner(particle, path));
848 }
849 
850 
851 void G4NucleiModel::
853  G4ElementaryParticleCollider* theEPCollider,
854  std::vector<G4CascadParticle>& outgoing_cparticles) {
855  if (verboseLevel > 1)
856  G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
857 
858  if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
859 
860  // Create four-vector checking
861 #ifdef G4CASCADE_CHECK_ECONS
862  G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV
863  balance.setVerboseLevel(verboseLevel);
864 #endif
865 
866  outgoing_cparticles.clear(); // Clear return buffer for this event
867  generateInteractionPartners(cparticle); // Fills "thePartners" data
868 
869  if (thePartners.empty()) { // smth. is wrong -> needs special treatment
870  if (verboseLevel)
871  G4cerr << " generateParticleFate-> got empty interaction-partners list "
872  << G4endl;
873  return;
874  }
875 
876  G4int npart = thePartners.size(); // Last item is a total-path placeholder
877 
878  if (npart == 1) { // cparticle is on the next zone entry
879  if (verboseLevel > 1)
880  G4cout << " no interactions; moving to next zone" << G4endl;
881 
883  cparticle.incrementCurrentPath(thePartners[0].second);
884  boundaryTransition(cparticle);
885  outgoing_cparticles.push_back(cparticle);
886 
887  if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
888 
889  //A.R. 19-Jun-2013: Fixed rare cases of non-reproducibility.
890  current_nucl1 = 0;
891  current_nucl2 = 0;
892 
893  return;
894  } // if (npart == 1)
895 
896  if (verboseLevel > 1)
897  G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
898 
899  G4ThreeVector old_position = cparticle.getPosition();
900  G4InuclElementaryParticle& bullet = cparticle.getParticle();
901  G4bool no_interaction = true;
902  G4int zone = cparticle.getCurrentZone();
903 
904  for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder
905  if (i > 0) cparticle.updatePosition(old_position);
906 
908 
909  if (verboseLevel > 3) {
910  if (target.quasi_deutron()) G4cout << " try absorption: ";
911  G4cout << " target " << target.type() << " bullet " << bullet.type()
912  << G4endl;
913  }
914 
915  EPCoutput.reset();
916  theEPCollider->collide(&bullet, &target, EPCoutput);
917 
918  // If collision failed, exit loop over partners
919  if (EPCoutput.numberOfOutgoingParticles() == 0) break;
920 
921  if (verboseLevel > 2) {
922  EPCoutput.printCollisionOutput();
923 #ifdef G4CASCADE_CHECK_ECONS
924  balance.collide(&bullet, &target, EPCoutput);
925  balance.okay(); // Do checks, but ignore result
926 #endif
927  }
928 
929  // Get list of outgoing particles for evaluation
930  std::vector<G4InuclElementaryParticle>& outgoing_particles =
931  EPCoutput.getOutgoingParticles();
932 
933  if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
934 
935  // Trailing effect: reject interaction at previously hit nucleon
937  const G4ThreeVector& new_position = cparticle.getPosition();
938 
939  if (!passTrailing(new_position)) continue;
940  collisionPts.push_back(new_position);
941 
942  // Sort particles according to beta (fastest first)
943  std::sort(outgoing_particles.begin(), outgoing_particles.end(),
945 
946  if (verboseLevel > 2)
947  G4cout << " adding " << outgoing_particles.size()
948  << " output particles" << G4endl;
949 
950  // NOTE: Embedded temporary is optimized away (no copying gets done)
951  G4int nextGen = cparticle.getGeneration()+1;
952  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) {
953  outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
954  new_position, zone,
955  0.0, nextGen));
956  }
957 
958  no_interaction = false;
959  current_nucl1 = 0;
960  current_nucl2 = 0;
961 
962 #ifdef G4CASCADE_DEBUG_CHARGE
963  {
964  G4double out_charge = 0.0;
965 
966  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
967  out_charge += outgoing_particles[ip].getCharge();
968 
969  G4cout << " multiplicity " << outgoing_particles.size()
970  << " bul type " << bullet.type()
971  << " targ type " << target.type()
972  << "\n initial charge "
973  << bullet.getCharge() + target.getCharge()
974  << " out charge " << out_charge << G4endl;
975  }
976 #endif
977 
978  if (verboseLevel > 2)
979  G4cout << " partner type " << target.type() << G4endl;
980 
981  if (target.nucleon()) {
982  current_nucl1 = target.type();
983  } else {
984  if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
985  current_nucl1 = (target.type() - 100) / 10;
986  current_nucl2 = target.type() - 100 - 10 * current_nucl1;
987  }
988 
989  if (current_nucl1 == 1) {
990  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
991  protonNumberCurrent--;
992  } else {
993  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
994  neutronNumberCurrent--;
995  }
996 
997  if (current_nucl2 == 1) {
998  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
999  protonNumberCurrent--;
1000  } else if (current_nucl2 == 2) {
1001  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
1002  neutronNumberCurrent--;
1003  }
1004 
1005  break;
1006  } // loop over partners
1007 
1008  if (no_interaction) { // still no interactions
1009  if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
1010 
1011  // For conservation checking (below), get particle before updating
1012  static G4ThreadLocal G4InuclElementaryParticle *prescatCP_G4MT_TLS_ = 0;
1013  if (!prescatCP_G4MT_TLS_)
1014  prescatCP_G4MT_TLS_ = new G4InuclElementaryParticle;
1015  G4InuclElementaryParticle &prescatCP = *prescatCP_G4MT_TLS_; // Avoid memory churn
1016  prescatCP = cparticle.getParticle();
1017 
1018  // Last "partner" is just a total-path placeholder
1019  cparticle.updatePosition(old_position);
1020  cparticle.propagateAlongThePath(thePartners[npart-1].second);
1021  cparticle.incrementCurrentPath(thePartners[npart-1].second);
1022  boundaryTransition(cparticle);
1023  outgoing_cparticles.push_back(cparticle);
1024 
1025  // Check conservation for simple scattering (ignore target nucleus!)
1026 #ifdef G4CASCADE_CHECK_ECONS
1027  if (verboseLevel > 2) {
1028  balance.collide(&prescatCP, 0, outgoing_cparticles);
1029  balance.okay(); // Report violations, but don't act on them
1030  }
1031 #endif
1032  } // if (no_interaction)
1033 
1034  return;
1035 }
1036 
1037 // Test if particle is suitable for absorption with dibaryons
1038 
1040  if (qdtype==pn || qdtype==0) // All absorptive particles
1041  return (ptype==pi0 || ptype==pip || ptype==pim || ptype==gam || ptype==mum);
1042  else if (qdtype==pp) // Negative or neutral only
1043  return (ptype==pi0 || ptype==pim || ptype==gam || ptype==mum);
1044  else if (qdtype==nn) // Positive or neutral only
1045  return (ptype==pi0 || ptype==pip || ptype==gam);
1046 
1047  return false; // Not a quasideuteron target
1048 }
1049 
1050 
1051 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles,
1052  G4int zone) {
1053  if (verboseLevel > 1) {
1054  G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
1055  }
1056 
1057  // Only check Fermi momenta for nucleons
1058  for (G4int i = 0; i < G4int(particles.size()); i++) {
1059  if (!particles[i].nucleon()) continue;
1060 
1061  G4int type = particles[i].type();
1062  G4double mom = particles[i].getMomModule();
1063  G4double pfermi = fermi_momenta[type-1][zone];
1064 
1065  if (verboseLevel > 2)
1066  G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
1067 
1068  if (mom < pfermi) {
1069  if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
1070  return false;
1071  }
1072  }
1073  return true;
1074 }
1075 
1076 
1077 // Test here for trailing effect: loop over all previous collision
1078 // locations and test for d > R_nucleon
1079 
1081  if (verboseLevel > 1)
1082  G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
1083 
1084  G4double dist;
1085  for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
1086  dist = (collisionPts[i] - hit_position).mag();
1087  if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
1088  if (dist < R_nucleon) {
1089  if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
1090  return false;
1091  }
1092  }
1093  return true; // New point far enough away to be used
1094 }
1095 
1096 
1098  if (verboseLevel > 1) {
1099  G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
1100  }
1101 
1102  G4int zone = cparticle.getCurrentZone();
1103 
1104  if (cparticle.movingInsideNuclei() && zone == 0) {
1105  if (verboseLevel) G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
1106  return;
1107  }
1108 
1109  G4LorentzVector mom = cparticle.getMomentum();
1110  G4ThreeVector pos = cparticle.getPosition();
1111 
1112  G4int type = cparticle.getParticle().type();
1113 
1114  G4double r = pos.mag();
1115  G4double pr = pos.dot(mom.vect()) / r;
1116 
1117  G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
1118 
1119  G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
1120  if (verboseLevel > 3) {
1121  G4cout << "Potentials for type " << type << " = "
1122  << getPotential(type,zone) << " , "
1123  << getPotential(type,next_zone) << G4endl;
1124  }
1125 
1126  G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
1127 
1128  G4double p1r;
1129 
1130  if (verboseLevel > 3) {
1131  G4cout << " type " << type << " zone " << zone << " next " << next_zone
1132  << " qv " << qv << " dv " << dv << G4endl;
1133  }
1134 
1135  if (qv <= 0.0) { // reflection
1136  if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
1137  p1r = -pr;
1138  cparticle.incrementReflectionCounter();
1139  } else { // transition
1140  if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
1141  p1r = std::sqrt(qv);
1142  if(pr < 0.0) p1r = -p1r;
1143  cparticle.updateZone(next_zone);
1144  cparticle.resetReflection();
1145  }
1146 
1147  G4double prr = (p1r - pr) / r;
1148 
1149  if (verboseLevel > 3) {
1150  G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
1151  << prr*pos.y() << " pz " << prr*pos.z() << " mag "
1152  << std::fabs(prr*r) << G4endl;
1153  }
1154 
1155  mom.setVect(mom.vect() + pos*prr);
1156  cparticle.updateParticleMomentum(mom);
1157 }
1158 
1159 
1160 // Select random point along full trajectory through nucleus
1161 // NOTE: Intended for projectile photons for initial interaction
1162 
1164  if (verboseLevel > 1)
1165  G4cout << " >>> G4NucleiModel::choosePointAlongTraj" << G4endl;
1166 
1167  // Get trajectory through nucleus by computing exit point of line,
1168  // assuming that current position is on surface
1169 
1170  // FIXME: These need to be reusable (static) buffers
1171  G4ThreeVector pos = cparticle.getPosition();
1172  G4ThreeVector rhat = pos.unit();
1173 
1174  G4ThreeVector phat = cparticle.getMomentum().vect().unit();
1175  if (cparticle.getMomentum().vect().mag() < small) phat.set(0.,0.,1.);
1176 
1177  if (verboseLevel > 3)
1178  G4cout << " pos " << pos << " phat " << phat << " rhat " << rhat << G4endl;
1179 
1180  G4ThreeVector posout = pos;
1181  G4double prang = rhat.angle(-phat);
1182 
1183  if (prang < 1e-6) posout = -pos; // Check for radial incidence
1184  else {
1185  G4double posrot = 2.*prang - pi;
1186  posout.rotate(posrot, phat.cross(rhat));
1187  if (verboseLevel > 3) G4cout << " posrot " << posrot/deg << " deg";
1188  }
1189 
1190  if (verboseLevel > 3) G4cout << " posout " << posout << G4endl;
1191 
1192  // Get list of zone crossings along trajectory
1193  G4ThreeVector posmid = (pos+posout)/2.; // Midpoint of trajectory
1194  G4double r2mid = posmid.mag2();
1195  G4double lenmid = (posout-pos).mag()/2.; // Half-length of trajectory
1196 
1197  G4int zoneout = number_of_zones-1;
1198  G4int zonemid = getZone(posmid.mag()); // Middle zone traversed
1199 
1200  // Every zone is entered then exited, so preallocate vector
1201  G4int ncross = (number_of_zones-zonemid)*2;
1202 
1203  if (verboseLevel > 3) {
1204  G4cout << " posmid " << posmid << " lenmid " << lenmid
1205  << " zoneout " << zoneout << " zonemid " << zonemid
1206  << " ncross " << ncross << G4endl;
1207  }
1208 
1209  // FIXME: These need to be reusable (static) buffers
1210  std::vector<G4double> wtlen(ncross,0.); // CDF from entry point
1211  std::vector<G4double> len(ncross,0.); // Distance from entry point
1212 
1213  // Work from outside in, to accumulate CDF steps properly
1214  G4int i; // Loop variable, used multiple times
1215  for (i=0; i<ncross/2; i++) {
1216  G4int iz = zoneout-i;
1217  G4double ds = std::sqrt(zone_radii[iz]*zone_radii[iz]-r2mid);
1218 
1219  len[i] = lenmid - ds; // Distance from entry to crossing
1220  len[ncross-1-i] = lenmid + ds; // Distance to outbound crossing
1221 
1222  if (verboseLevel > 3) {
1223  G4cout << " i " << i << " iz " << iz << " ds " << ds
1224  << " len " << len[i] << G4endl;
1225  }
1226  }
1227 
1228  // Compute weights for each zone along trajectory and integrate
1229  for (i=1; i<ncross; i++) {
1230  G4int iz = (i<ncross/2) ? zoneout-i+1 : zoneout-ncross+i+1;
1231 
1232  G4double dlen = len[i]-len[i-1]; // Distance from previous crossing
1233 
1234  // Uniform probability across entire zone
1235  //*** G4double wt = dlen*(getDensity(1,iz)+getDensity(2,iz));
1236 
1237  // Probability based on interaction length through zone
1238  G4double invmfp = (inverseMeanFreePath(cparticle, neutronEP, iz)
1239  + inverseMeanFreePath(cparticle, protonEP, iz));
1240 
1241  // Integral of exp(-len/mfp) from start of zone to end
1242  G4double wt = (G4Exp(-len[i-1]*invmfp)-G4Exp(-len[i]*invmfp)) / invmfp;
1243 
1244  wtlen[i] = wtlen[i-1] + wt;
1245 
1246  if (verboseLevel > 3) {
1247  G4cout << " i " << i << " iz " << iz << " avg.mfp " << 1./invmfp
1248  << " dlen " << dlen << " wt " << wt << " wtlen " << wtlen[i]
1249  << G4endl;
1250  }
1251  }
1252 
1253  // Normalize CDF to unit integral
1254  std::transform(wtlen.begin(), wtlen.end(), wtlen.begin(),
1255  std::bind2nd(std::divides<G4double>(), wtlen.back()));
1256 
1257  if (verboseLevel > 3) {
1258  G4cout << " weights";
1259  for (i=0; i<ncross; i++) G4cout << " " << wtlen[i];
1260  G4cout << G4endl;
1261  }
1262 
1263  // Choose random point along trajectory, weighted by density
1264  G4double rand = G4UniformRand();
1265  G4int ir = std::upper_bound(wtlen.begin(),wtlen.end(),rand) - wtlen.begin();
1266 
1267  G4double frac = (rand-wtlen[ir-1]) / (wtlen[ir]-wtlen[ir-1]);
1268  G4double drand = (1.-frac)*len[ir-1] + frac*len[ir];
1269 
1270  if (verboseLevel > 3) {
1271  G4cout << " rand " << rand << " ir " << ir << " frac " << frac
1272  << " drand " << drand << G4endl;
1273  }
1274 
1275  pos += drand * phat; // Distance from entry point along trajectory
1276 
1277  // Update input particle with new location
1278  cparticle.updatePosition(pos);
1279  cparticle.updateZone(getZone(pos.mag()));
1280 
1281  if (verboseLevel > 2) {
1282  G4cout << " moved particle to zone " << cparticle.getCurrentZone()
1283  << " @ " << pos << G4endl;
1284  }
1285 }
1286 
1287 
1288 // Returns true if particle should interact immediately
1290  return (isProjectile(cparticle) &&
1291  (cparticle.getParticle().isPhoton() ||
1292  cparticle.getParticle().isMuon())
1293  );
1294 }
1295 
1297  return (cparticle.getGeneration() == 0); // Only initial state particles
1298 }
1299 
1301  if (verboseLevel > 1) {
1302  G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
1303  }
1304 
1305  const G4double ekin_scale = 2.0;
1306 
1307  G4bool worth = true;
1308 
1309  if (cparticle.reflectedNow()) { // Just reflected -- keep going?
1310  G4int zone = cparticle.getCurrentZone();
1311  G4int ip = cparticle.getParticle().type();
1312 
1313  // NOTE: Temporarily backing out use of potential for non-nucleons
1314  G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
1315  getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
1316 
1317  worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
1318 
1319  if (verboseLevel > 3) {
1320  G4cout << " type=" << ip
1321  << " ekin=" << cparticle.getParticle().getKineticEnergy()
1322  << " potential=" << ekin_cut
1323  << " : worth? " << worth << G4endl;
1324  }
1325  }
1326 
1327  return worth;
1328 }
1329 
1330 
1332  if (verboseLevel > 4) {
1333  G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
1334  }
1335 
1336  switch (ip) {
1337  case proton: return G4double(protonNumberCurrent)/G4double(protonNumber);
1338  case neutron: return G4double(neutronNumberCurrent)/G4double(neutronNumber);
1339  case diproton: return getRatio(proton)*getRatio(proton);
1340  case unboundPN: return getRatio(proton)*getRatio(neutron);
1341  case dineutron: return getRatio(neutron)*getRatio(neutron);
1342  default: return 0.;
1343  }
1344 
1345  return 0.;
1346 }
1347 
1349  const G4double pn_spec = 1.0; // Scale factor for pn vs. pp/nn
1350  //const G4double pn_spec = 0.5;
1351 
1352  G4double dens = 0.;
1353 
1354  if (ip < 100) dens = getDensity(ip,izone);
1355  else { // For dibaryons, remove extra 1/volume term in density product
1356  switch (ip) {
1357  case diproton:
1358  dens = getDensity(proton,izone) * getDensity(proton,izone);
1359  break;
1360  case unboundPN:
1361  dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
1362  break;
1363  case dineutron:
1364  dens = getDensity(neutron,izone) * getDensity(neutron,izone);
1365  break;
1366  default: dens = 0.;
1367  }
1368  dens *= getVolume(izone);
1369  }
1370 
1371  return getRatio(ip) * dens;
1372 }
1373 
1374 
1377  if (verboseLevel > 1) {
1378  G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
1379  }
1380 
1381  // FIXME: Previous version generated random sin(theta), then used -cos(theta)
1382  // Using generateWithRandomAngles changes result!
1383  // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
1384  G4double costh = std::sqrt(1.0 - inuclRndm());
1385  G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
1386 
1387  // Start particle outside nucleus, unless capture-at-rest
1388  G4int zone = number_of_zones;
1389  if (particle->getKineticEnergy() < small) zone--;
1390 
1391  G4CascadParticle cpart(*particle, pos, zone, large, 0);
1392 
1393  // SPECIAL CASE: Inbound photons are emplanted along through-path
1394  if (forceFirst(cpart)) choosePointAlongTraj(cpart);
1395 
1396  if (verboseLevel > 2) G4cout << cpart << G4endl;
1397 
1398  return cpart;
1399 }
1400 
1403  modelLists& output) {
1404  if (verboseLevel) {
1405  G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
1406  << G4endl;
1407  }
1408 
1409  const G4double max_a_for_cascad = 5.0;
1410  const G4double ekin_cut = 2.0;
1411  const G4double small_ekin = 1.0e-6;
1412  const G4double r_large2for3 = 62.0;
1413  const G4double r0forAeq3 = 3.92;
1414  const G4double s3max = 6.5;
1415  const G4double r_large2for4 = 69.14;
1416  const G4double r0forAeq4 = 4.16;
1417  const G4double s4max = 7.0;
1418  const G4int itry_max = 100;
1419 
1420  // Convenient references to input buffer contents
1421  std::vector<G4CascadParticle>& casparticles = output.first;
1422  std::vector<G4InuclElementaryParticle>& particles = output.second;
1423 
1424  casparticles.clear(); // Reset input buffer to fill this event
1425  particles.clear();
1426 
1427  // first decide whether it will be cascad or compound final nuclei
1428  G4int ab = bullet->getA();
1429  G4int zb = bullet->getZ();
1430  G4int at = target->getA();
1431  G4int zt = target->getZ();
1432 
1433  G4double massb = bullet->getMass(); // For creating LorentzVectors below
1434 
1435  if (ab < max_a_for_cascad) {
1436 
1437  G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
1438  G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
1439  G4double ben = benb < bent ? bent : benb;
1440 
1441  if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
1442  G4int itryg = 0;
1443 
1444  while (casparticles.size() == 0 && itryg < itry_max) {
1445  itryg++;
1446  particles.clear();
1447 
1448  // nucleons coordinates and momenta in nuclei rest frame
1449  coordinates.clear();
1450  momentums.clear();
1451 
1452  if (ab < 3) { // deuteron, simplest case
1453  G4double r = 2.214 - 3.4208 * G4Log(1.0 - 0.981 * inuclRndm());
1455  coordinates.push_back(coord1);
1456  coordinates.push_back(-coord1);
1457 
1458  G4double p = 0.0;
1459  G4bool bad = true;
1460  G4int itry = 0;
1461 
1462  while (bad && itry < itry_max) {
1463  itry++;
1464  p = 456.0 * inuclRndm();
1465 
1466  if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1467  p * r > 312.0) bad = false;
1468  }
1469 
1470  if (itry == itry_max)
1471  if (verboseLevel > 2){
1472  G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;
1473  }
1474 
1475  p = 0.0005 * p;
1476 
1477  if (verboseLevel > 2){
1478  G4cout << " p nuc " << p << G4endl;
1479  }
1480 
1481  G4LorentzVector mom = generateWithRandomAngles(p, massb);
1482 
1483  momentums.push_back(mom);
1484  mom.setVect(-mom.vect());
1485  momentums.push_back(-mom);
1486  } else {
1487  G4ThreeVector coord1;
1488 
1489  G4bool badco = true;
1490 
1491  G4int itry = 0;
1492 
1493  if (ab == 3) {
1494  while (badco && itry < itry_max) {
1495  if (itry > 0) coordinates.clear();
1496  itry++;
1497  G4int i(0);
1498 
1499  for (i = 0; i < 2; i++) {
1500  G4int itry1 = 0;
1501  G4double ss, u, rho;
1502  G4double fmax = G4Exp(-0.5) / std::sqrt(0.5);
1503 
1504  while (itry1 < itry_max) {
1505  itry1++;
1506  ss = -G4Log(inuclRndm());
1507  u = fmax * inuclRndm();
1508  rho = std::sqrt(ss) * G4Exp(-ss);
1509 
1510  if (rho > u && ss < s3max) {
1511  ss = r0forAeq3 * std::sqrt(ss);
1512  coord1 = generateWithRandomAngles(ss).vect();
1513  coordinates.push_back(coord1);
1514 
1515  if (verboseLevel > 2){
1516  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1517  }
1518  break;
1519  }
1520  }
1521 
1522  if (itry1 == itry_max) { // bad case
1523  coord1.set(10000.,10000.,10000.);
1524  coordinates.push_back(coord1);
1525  break;
1526  }
1527  }
1528 
1529  coord1 = -coordinates[0] - coordinates[1];
1530  if (verboseLevel > 2) {
1531  G4cout << " 3 r " << coord1.mag() << G4endl;
1532  }
1533 
1534  coordinates.push_back(coord1);
1535 
1536  G4bool large_dist = false;
1537 
1538  for (i = 0; i < 2; i++) {
1539  for (G4int j = i+1; j < 3; j++) {
1540  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1541 
1542  if (verboseLevel > 2) {
1543  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1544  }
1545 
1546  if (r2 > r_large2for3) {
1547  large_dist = true;
1548 
1549  break;
1550  }
1551  }
1552 
1553  if (large_dist) break;
1554  }
1555 
1556  if(!large_dist) badco = false;
1557 
1558  }
1559 
1560  } else { // a >= 4
1561  G4double b = 3./(ab - 2.0);
1562  G4double b1 = 1.0 - b / 2.0;
1563  G4double u = b1 + std::sqrt(b1 * b1 + b);
1564  G4double fmax = (1.0 + u/b) * u * G4Exp(-u);
1565 
1566  while (badco && itry < itry_max) {
1567 
1568  if (itry > 0) coordinates.clear();
1569  itry++;
1570  G4int i(0);
1571 
1572  for (i = 0; i < ab-1; i++) {
1573  G4int itry1 = 0;
1574  G4double ss;
1575 
1576  while (itry1 < itry_max) {
1577  itry1++;
1578  ss = -G4Log(inuclRndm());
1579  u = fmax * inuclRndm();
1580 
1581  if (std::sqrt(ss) * G4Exp(-ss) * (1.0 + ss/b) > u
1582  && ss < s4max) {
1583  ss = r0forAeq4 * std::sqrt(ss);
1584  coord1 = generateWithRandomAngles(ss).vect();
1585  coordinates.push_back(coord1);
1586 
1587  if (verboseLevel > 2) {
1588  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1589  }
1590 
1591  break;
1592  }
1593  }
1594 
1595  if (itry1 == itry_max) { // bad case
1596  coord1.set(10000.,10000.,10000.);
1597  coordinates.push_back(coord1);
1598  break;
1599  }
1600  }
1601 
1602  coord1 *= 0.0; // Cheap way to reset
1603  for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
1604 
1605  coordinates.push_back(coord1);
1606 
1607  if (verboseLevel > 2){
1608  G4cout << " last r " << coord1.mag() << G4endl;
1609  }
1610 
1611  G4bool large_dist = false;
1612 
1613  for (i = 0; i < ab-1; i++) {
1614  for (G4int j = i+1; j < ab; j++) {
1615 
1616  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1617 
1618  if (verboseLevel > 2){
1619  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1620  }
1621 
1622  if (r2 > r_large2for4) {
1623  large_dist = true;
1624 
1625  break;
1626  }
1627  }
1628 
1629  if (large_dist) break;
1630  }
1631 
1632  if (!large_dist) badco = false;
1633  }
1634  }
1635 
1636  if(badco) {
1637  G4cout << " can not generate the nucleons coordinates for a "
1638  << ab << G4endl;
1639 
1640  casparticles.clear(); // Return empty buffer on error
1641  particles.clear();
1642  return;
1643 
1644  } else { // momentums
1645  G4double p, u, x;
1646  G4LorentzVector mom;
1647  //G4bool badp = True;
1648 
1649  for (G4int i = 0; i < ab - 1; i++) {
1650  G4int itry2 = 0;
1651 
1652  while(itry2 < itry_max) {
1653  itry2++;
1654  u = -G4Log(0.879853 - 0.8798502 * inuclRndm());
1655  x = u * G4Exp(-u);
1656 
1657  if(x > inuclRndm()) {
1658  p = std::sqrt(0.01953 * u);
1659  mom = generateWithRandomAngles(p, massb);
1660  momentums.push_back(mom);
1661 
1662  break;
1663  }
1664  } // while (itry2 < itry_max)
1665 
1666  if(itry2 == itry_max) {
1667  G4cout << " can not generate proper momentum for a "
1668  << ab << G4endl;
1669 
1670  casparticles.clear(); // Return empty buffer on error
1671  particles.clear();
1672  return;
1673  }
1674  } // for (i=0 ...
1675  // last momentum
1676 
1677  mom *= 0.; // Cheap way to reset
1678  mom.setE(bullet->getEnergy()+target->getEnergy());
1679 
1680  for(G4int j=0; j< ab-1; j++) mom -= momentums[j];
1681 
1682  momentums.push_back(mom);
1683  }
1684  }
1685 
1686  // Coordinates and momenta at rest are generated, now back to the lab
1687  G4double rb = 0.0;
1688  G4int i(0);
1689 
1690  for(i = 0; i < G4int(coordinates.size()); i++) {
1691  G4double rp = coordinates[i].mag2();
1692 
1693  if(rp > rb) rb = rp;
1694  }
1695 
1696  // nuclei i.p. as a whole
1697  G4double s1 = std::sqrt(inuclRndm());
1698  G4double phi = randomPHI();
1699  G4double rz = (nuclei_radius + rb) * s1;
1700  G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
1701  -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
1702 
1703  for (i = 0; i < G4int(coordinates.size()); i++) {
1704  coordinates[i] += global_pos;
1705  }
1706 
1707  // all nucleons at rest
1708  raw_particles.clear();
1709 
1710  for (G4int ipa = 0; ipa < ab; ipa++) {
1711  G4int knd = ipa < zb ? 1 : 2;
1712  raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1713  }
1714 
1715  G4InuclElementaryParticle dummy(small_ekin, 1);
1716  G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
1717  toTheBulletRestFrame.toTheTargetRestFrame();
1718 
1720 
1721  for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1722  ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum()));
1723  }
1724 
1725  // fill cascad particles and outgoing particles
1726 
1727  for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1728  G4LorentzVector mom = raw_particles[ip].getMomentum();
1729  G4double pmod = mom.rho();
1730  G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
1731  G4double det = t0 * t0 + nuclei_radius * nuclei_radius
1732  - coordinates[ip].mag2();
1733  G4double tr = -1.0;
1734 
1735  if(det > 0.0) {
1736  G4double t1 = t0 + std::sqrt(det);
1737  G4double t2 = t0 - std::sqrt(det);
1738 
1739  if(std::fabs(t1) <= std::fabs(t2)) {
1740  if(t1 > 0.0) {
1741  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1742  }
1743 
1744  if(tr < 0.0 && t2 > 0.0) {
1745 
1746  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1747  }
1748 
1749  } else {
1750  if(t2 > 0.0) {
1751 
1752  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1753  }
1754 
1755  if(tr < 0.0 && t1 > 0.0) {
1756  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1757  }
1758  }
1759 
1760  }
1761 
1762  if(tr >= 0.0) { // cascad particle
1763  coordinates[ip] += mom.vect()*tr / pmod;
1764  casparticles.push_back(G4CascadParticle(raw_particles[ip],
1765  coordinates[ip],
1766  number_of_zones, large, 0));
1767 
1768  } else {
1769  particles.push_back(raw_particles[ip]);
1770  }
1771  }
1772  }
1773 
1774  if(casparticles.size() == 0) {
1775  particles.clear();
1776 
1777  G4cout << " can not generate proper distribution for " << itry_max
1778  << " steps " << G4endl;
1779  }
1780  }
1781  }
1782 
1783  if(verboseLevel > 2){
1784  G4int ip(0);
1785 
1786  G4cout << " cascad particles: " << casparticles.size() << G4endl;
1787  for(ip = 0; ip < G4int(casparticles.size()); ip++)
1788  G4cout << casparticles[ip] << G4endl;
1789 
1790  G4cout << " outgoing particles: " << particles.size() << G4endl;
1791  for(ip = 0; ip < G4int(particles.size()); ip++)
1792  G4cout << particles[ip] << G4endl;
1793  }
1794 
1795  return; // Buffer has been filled
1796 }
1797 
1798 
1799 // Compute mean free path for inclusive interaction of projectile and target
1800 G4double
1803  G4int zone) {
1804  G4int ptype = cparticle.getParticle().type();
1805  G4int ip = target.type();
1806 
1807  // Ensure that zone specified is within nucleus, for array lookups
1808  if (zone<0) zone = cparticle.getCurrentZone();
1809  if (zone>=number_of_zones) zone = number_of_zones-1;
1810 
1811  // Special cases: neutrinos, and muon-on-neutron, have infinite path
1812  if (cparticle.getParticle().isNeutrino()) return 0.;
1813  if (ptype == muonMinus && ip == neutron) return 0.;
1814 
1815  // Configure frame transformation to get kinetic energy for lookups
1816  dummy_convertor.setBullet(cparticle.getParticle());
1817  dummy_convertor.setTarget(&target);
1818  dummy_convertor.toTheCenterOfMass(); // Fill internal kinematics
1819  G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
1820 
1821  // Get cross section for interacting with target (dibaryons are absorptive)
1822  G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
1823  : absorptionCrossSection(ekin, ptype);
1824 
1825  if (verboseLevel > 2) {
1826  G4cout << " ip " << ip << " zone " << zone << " ekin " << ekin
1827  << " dens " << getCurrentDensity(ip, zone)
1828  << " csec " << csec << G4endl;
1829  }
1830 
1831  if (csec <= 0.) return 0.; // No interaction, avoid unnecessary work
1832 
1833  // Get nuclear density and compute mean free path
1834  return csec * getCurrentDensity(ip, zone);
1835 }
1836 
1837 // Throw random distance for interaction of particle using path and MFP
1838 
1839 G4double
1841  G4double path, G4double invmfp) const {
1842  // Delay interactions of newly formed secondaries (minimum int. length)
1843  const G4double young_cut = std::sqrt(10.0) * 0.25;
1844  const G4double huge_num = 50.0; // Argument to exponential
1845 
1846  G4double spath = large; // Buffer for return value
1847 
1848  if (invmfp < small) return spath; // No interaction, avoid unnecessary work
1849 
1850  G4double pw = -path * invmfp; // Ratio of path in zone to MFP
1851  if (pw < -huge_num) pw = -huge_num;
1852  pw = 1.0 - G4Exp(pw);
1853 
1854  if (verboseLevel > 2)
1855  G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
1856 
1857  // Primary particle(s) should always interact at least once
1858  if (forceFirst(cparticle) || (inuclRndm() < pw)) {
1859  spath = -G4Log(1.0 - pw * inuclRndm()) / invmfp;
1860  if (cparticle.young(young_cut, spath)) spath = large;
1861 
1862  if (verboseLevel > 2)
1863  G4cout << " spath " << spath << " path " << path << G4endl;
1864  }
1865 
1866  return spath;
1867 }
1868 
1869 
1870 // Parametrized cross section for pion and photon absorption
1871 
1873  if (!useQuasiDeuteron(type)) {
1874  G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
1875  return 0.;
1876  }
1877 
1878  G4double csec = 0.;
1879 
1880  // Pion absorption is parametrized for low vs. medium energy
1881  // ... use for muon capture as well
1882  if (type == pionPlus || type == pionMinus || type == pionZero ||
1883  type == muonMinus) {
1884  if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
1885  + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
1886  else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);
1887  }
1888 
1889  // Photon cross-section is binned from parametrization by Mi. Kossov
1890  // See below for initialization of gammaQDinterp, gammaQDxsec
1891  if (type == photon) {
1892  csec = gammaQDinterp.interpolate(ke, gammaQDxsec) * gammaQDscale;
1893  }
1894 
1895  if (csec < 0.0) csec = 0.0;
1896 
1897  if (verboseLevel > 2) {
1898  G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;
1899  }
1900 
1901  return crossSectionUnits * csec;
1902 }
1903 
1904 
1906 {
1907  // All scattering cross-sections are available from tables
1908  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
1909  if (!xsecTable) {
1910  G4cerr << " unknown collison type = " << rtype << G4endl;
1911  return 0.;
1912  }
1913 
1914  return (crossSectionUnits * xsecTable->getCrossSection(ke));
1915 }
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2, G4double nuclearRadius) const
void set(double x, double y, double z)
void updateParticleMomentum(const G4LorentzVector &mom)
G4InuclElementaryParticle generateQuasiDeuteron(G4int type1, G4int type2, G4int zone) const
std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
TTree * t1
Definition: plottest35.C:26
G4bool reflectedNow() const
static const G4CascadeChannel * GetTable(G4int initialState)
G4int getZ() const
static G4double xsecScale()
G4LorentzVector getMomentum() const
void incrementReflectionCounter()
void incrementCurrentPath(G4double npath)
double angle(const Hep3Vector &) const
tuple a
Definition: test.py:11
G4bool movingInsideNuclei() const
void generateInteractionPartners(G4CascadParticle &cparticle)
double x() const
double dot(const Hep3Vector &) const
G4int getGeneration() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static G4double radiusTrailing()
Int_t ipart
Definition: Style.C:10
const char * p
Definition: xmltok.h:285
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
static G4double radiusScale()
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4double getFermiKinetic(G4int ip, G4int izone) const
G4bool passFermi(const std::vector< G4InuclElementaryParticle > &particles, G4int zone)
G4double zoneIntegralGaussian(G4double ur1, G4double ur2, G4double nuclearRadius) const
void propagateAlongThePath(G4double path)
G4double getPathToTheNextZone(G4double rz_in, G4double rz_out)
std::vector< partner > thePartners
G4double getEnergy() const
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
void updatePosition(const G4ThreeVector &pos)
tuple x
Definition: test.py:50
void boundaryTransition(G4CascadParticle &cparticle)
subroutine sort(A, N)
Definition: dpm25nuc7.f:4670
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
const XML_Char * target
Definition: expat.h:268
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
#define G4ThreadLocal
Definition: tls.hh:52
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4CascadParticle initializeCascad(G4InuclElementaryParticle *particle)
G4bool young(G4double young_path_cut, G4double cpath) const
G4double getVolume(G4int izone) const
int G4int
Definition: G4Types.hh:78
G4double absorptionCrossSection(G4double e, G4int type) const
const G4InuclElementaryParticle & getParticle() const
double z() const
virtual void setVerboseLevel(G4int verbose=0)
virtual G4double getCrossSection(double ke) const =0
G4bool isProjectile(const G4CascadParticle &cparticle) const
G4double getKineticEnergy() const
void fillBindingEnergies()
Double_t y
Definition: plot.C:279
void updateZone(G4int izone)
tuple b
Definition: test.py:12
G4double getCurrentDensity(G4int ip, G4int izone) const
void fillZoneRadii(G4double nuclearRadius)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4int getA() const
void fillPotentials(G4int type, G4double tot_vol)
Float_t Z
Definition: plot.C:39
G4bool nucleon(G4int ityp)
bool G4bool
Definition: G4Types.hh:79
G4double iz
Definition: TRTMaterials.hh:39
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
G4int numberOfOutgoingParticles() const
double rho() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double generateInteractionLength(const G4CascadParticle &cparticle, G4double path, G4double invmfp) const
G4double fillZoneVolumes(G4double nuclearRadius)
void generateModel(G4InuclNuclei *nuclei)
std::pair< G4InuclElementaryParticle, G4double > partner
G4double getKinEnergyInTheTRS() const
void generateParticleFate(G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
static G4double fermiScale()
G4double getRatio(G4int ip) const
static G4double radiusAlpha()
jump r
Definition: plot.C:36
G4double totalCrossSection(G4double ke, G4int rtype) const
G4double getFermiMomentum(G4int ip, G4int izone) const
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
Hep3Vector unit() const
G4bool worthToPropagate(const G4CascadParticle &cparticle) const
static G4bool sortPartners(const partner &p1, const partner &p2)
double y() const
void choosePointAlongTraj(G4CascadParticle &cparticle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
G4int getCurrentZone() const
double mag2() const
static G4bool useTwoParam()
tuple z
Definition: test.py:28
G4double getCharge() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:64
TTree * t2
Definition: plottest35.C:36
const XML_Char int len
Definition: expat.h:262
void printModel() const
Int_t npart
Definition: Style.C:7
#define G4endl
Definition: G4ios.hh:61
Hep3Vector cross(const Hep3Vector &) const
G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const
const G4ThreeVector & getPosition() const
void setVect(const Hep3Vector &)
double G4double
Definition: G4Types.hh:76
static G4double radiusSmall()
virtual ~G4NucleiModel()
G4double getDensity(G4int ip, G4int izone) const
G4double bindingEnergy(G4int A, G4int Z)
double mag() const
G4bool passTrailing(const G4ThreeVector &hit_position)
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4int getZone(G4double r) const
G4bool forceFirst(const G4CascadParticle &cparticle) const
static G4double gammaQDScale()
G4GLOB_DLL std::ostream G4cerr
G4double getMass() const
void setTarget(const G4InuclParticle *target)
G4double getPotential(G4int ip, G4int izone) const