Geant4  10.02
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 // For the best approximation to a physical-units model, set the following:
29 // setenv G4NUCMODEL_XSEC_SCALE 0.1
30 // setenv G4NUCMODEL_RAD_SCALE 1.0
31 // setenv G4NUCMODEL_RAD_2PAR 1
32 // setenv G4NUCMODEL_RAD_SMALL 1.992
33 // setenv G4NUCMODEL_RAD_ALPHA 0.84
34 // setenv G4NUCMODEL_FERMI_SCALE 0.685
35 // setenv G4NUCMODEL_RAD_TRAILING 0.70
36 //
37 // 20100112 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
38 // 20100114 M. Kelsey -- Use G4ThreeVector for position
39 // 20100309 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
40 // eliminate some unnecessary std::pow(), move ctor here
41 // 20100407 M. Kelsey -- Replace std::vector<>::resize(0) with ::clear().
42 // Move output vectors from generateParticleFate() and
43 // ::generateInteractionPartners() to be data members; return via
44 // const-ref instead of by value.
45 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
46 // 20100418 M. Kelsey -- Reference output particle lists via const-ref
47 // 20100421 M. Kelsey -- Replace hardwired p/n masses with G4PartDef's
48 // 20100517 M. Kelsey -- Use G4CascadeINterpolator for cross-section
49 // calculations. use G4Cascade*Channel for total xsec in pi-N
50 // N-N channels. Move absorptionCrossSection() from SpecialFunc.
51 // 20100610 M. Kelsey -- Replace another random-angle code block; add some
52 // diagnostic output for partner-list production.
53 // 20100617 M. Kelsey -- Replace preprocessor flag CHC_CHECK with
54 // G4CASCADE_DEBUG_CHARGE
55 // 20100620 M. Kelsey -- Improve error message on empty partners list, add
56 // four-momentum checking after EPCollider
57 // 20100621 M. Kelsey -- In boundaryTransition() account for momentum transfer
58 // to secondary by boosting into recoil nucleus "rest" frame.
59 // Don't need temporaries to create dummy "partners" for list.
60 // 20100622 M. Kelsey -- Restore use of "bindingEnergy()" function name, which
61 // is now a wrapper for G4NucleiProperties::GetBindingEnergy().
62 // 20100623 M. Kelsey -- Eliminate some temporaries terminating partner-list,
63 // discard recoil boost for now. Initialize all data
64 // members in ctors. Allow generateModel() to be called
65 // mutliple times, by clearing vectors each time through;
66 // avoid extra work by returning if A and Z are same as
67 // before.
68 // 20100628 M. Kelsey -- Two momentum-recoil bugs; don't subtract energies!
69 // 20100715 M. Kelsey -- Make G4InuclNuclei in generateModel(), use for
70 // balance checking (only verbose>2) in generateParticleFate().
71 // 20100721 M. Kelsey -- Use new G4CASCADE_CHECK_ECONS for conservation checks
72 // 20100723 M. Kelsey -- Move G4CollisionOutput buffer to .hh for reuse
73 // 20100726 M. Kelsey -- Preallocate arrays with number_of_zones dimension.
74 // 20100902 M. Kelsey -- Remove resize(3) directives from qdeutron/acsecs
75 // 20100907 M. Kelsey -- Limit interaction targets based on current nucleon
76 // counts (e.g., no pp if protonNumberCurrent < 2!)
77 // 20100914 M. Kelsey -- Migrate to integer A and Z
78 // 20100928 M. Kelsey -- worthToPropagate() use nuclear potential for hadrons.
79 // 20101005 M. Kelsey -- Annotate hardwired numbers for upcoming rationalizing,
80 // move hardwired numbers out to static data members, factor out
81 // all the model construction pieces to separate functions, fix
82 // wrong value for 4/3 pi.
83 // 20101019 M. Kelsey -- CoVerity reports: unitialized constructor, dtor leak
84 // 20101020 M. Kelsey -- Bug fixes to refactoring changes (5 Oct). Back out
85 // worthToPropagate() changes for better regression testing.
86 // 20101020 M. Kelsey -- Re-activate worthToPropagate() changes.
87 // 20101119 M. Kelsey -- Hide "negative path" and "no partners" messages in
88 // verbosity.
89 // 20110218 M. Kelsey -- Add crossSectionUnits and radiusUnits scale factors,
90 // use "theoretical" numbers for radii etc., multipled by scale
91 // factor; set scale factors using environment variables
92 // 20110303 M. Kelsey -- Add comments why using fabs() with B.E. differences?
93 // 20110321 M. Kelsey -- Replace strtof() with strtod() for envvar conversion
94 // 20110321 M. Kelsey -- Use fm and fm^2 as default units, Per D. Wright
95 // (NOTE: Restored from original 20110318 commit)
96 // 20110324 D. Wright -- Implement trailing effect
97 // 20110324 M. Kelsey -- Move ::reset() here, as it has more code.
98 // 20110519 M. Kelsey -- Used "rho" after assignment, instead of recomputing
99 // 20110525 M. Kelsey -- Revert scale factor changes (undo 20110321 changes)
100 // 20110617 M. Kelsey -- Apply scale factor to trailing-effect radius, make
101 // latter runtime adjustable (G4NUCMODEL_RAD_TRAILING)
102 // 20110720 M. Kelsey -- Follow interface change for cross-section tables,
103 // eliminating switch blocks.
104 // 20110806 M. Kelsey -- Reduce memory churn by pre-allocating buffers
105 // 20110823 M. Kelsey -- Remove local cross-section tables entirely
106 // 20110825 M. Kelsey -- Add comments regarding Fermi momentum scale, set of
107 // "best guess" parameter values
108 // 20110831 M. Kelsey -- Make "best guess" parameters the defaults
109 // 20110922 M. Kelsey -- Follow migrations G4InuclParticle::print(ostream&)
110 // and G4CascadParticle::print(ostream&)
111 // 20111018 M. Kelsey -- Correct kaon potential to be positive, not negative
112 // 20111107 M. Kelsey -- *** REVERT TO OLD NON-PHYSICAL PARAMETERS FOR 9.5 ***
113 // 20120306 M. Kelsey -- For incident projectile (CParticle incoming to
114 // nucleus from outside) don't use pw cut; force interaction.
115 // 20120307 M. Kelsey -- Compute zone volumes during setup, not during each
116 // interaction. Include photons in absorption by dibaryons.
117 // Discard unused code in totalCrossSection(). Encapsulate
118 // interaction path calculations in generateInteractionLength().
119 // 20120308 M. Kelsey -- Add binned photon absorption cross-section.
120 // 20120320 M. Kelsey -- Bug fix: fill zone_volumes for single-zone case
121 // 20120321 M. Kelsey -- Add check in isProjectile() for single-zone case.
122 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
123 // 20120822 M. Kelsey -- Move envvars to G4CascadeParameters.
124 // 20121205 M. Kelsey -- Bug fix in generateParticleFate(): daughters should
125 // have generation count incremented from parent. This allows
126 // isProjectile() to simply check generation number == 0.
127 // 20130221 M. Kelsey -- For incident photons, move to random point along
128 // trajectory through nucleus for first (forced) interaction.
129 // 20130223 M. Kelsey -- Fix bugs in weighting and selection of traj points
130 // 20130302 M. Kelsey -- Replace "isProjectile()" with "forceFirst()" in
131 // path-length calculation.
132 // 20130314 M. Kelsey -- Restore null initializer and if-block for _TLS_.
133 // 20130508 D. Wright -- Support muon-nuclear interactions.
134 // 20130510 M. Kelsey -- Check for zero-momentum in choosePointAlongTraj().
135 // 20130511 M. Kelsey -- Move choosePointAlongTraj() to initializeCascad(),
136 // instead of after generateIP(). Change "spath<path" to <= to
137 // handle at-rest particles. Skip mu-/neutron interactions.
138 // Hide "zone 0 transition" error message behind verbosity.
139 // 20130523 M. Kelsey -- Bug fix: Initialize dummy_convertor in inverseMFP.
140 // For capture-at-rest, replace exterior nucleus with outer zone.
141 // 20130524 M. Kelsey -- Move "large" and "small" cutoffs to static const.
142 // Check zone argument to inverseMFP() to ensure within range.
143 // 20130611 M. Kelsey -- Undo "spath<=path" change (20130511), replace with
144 // explicit check on spath==0.
145 // 20130619 A. Ribon -- Fixed reproducibility problem in the method
146 // generateParticleFate
147 // 20130627 M. Kelsey -- Check "path==0.", rather than spath.
148 // 20130628 M. Kelsey -- Print deuteron type code, rather than array index,
149 // Extend useQuasiDeuteron() to check good absorption
150 // 20130701 M. Kelsey -- Don't average 1/MFP for total interaction; just sum;
151 // inverseMFP() returns "large" value instead of input path.
152 // 20130806 M. Kelsey -- Move static G4InuclEP's to file scope to avoid
153 // thread collisions
154 // 20130924 M. Kelsey -- Use G4Log, G4Exp, G4Pow for CPU speedup
155 // 20130925 M. Kelsey -- Eliminate unnecessary use of pow() in integrals
156 // 20131001 M. Kelsey -- Move QDinterp object to data member, thread local
157 // 20140116 M. Kelsey -- Move all envvar parameters to const data members
158 // 20141001 M. Kelsey -- Change sign of "dv" in boundaryTransition
159 // 20150608 M. Kelsey -- Label all while loops as terminating.
160 // 20150622 M. Kelsey -- Use G4AutoDelete for _TLS_ buffers.
161 
162 #include "G4NucleiModel.hh"
163 #include "G4AutoDelete.hh"
164 #include "G4CascadeChannel.hh"
165 #include "G4CascadeChannelTables.hh"
166 #include "G4CascadeCheckBalance.hh"
167 #include "G4CascadeInterpolator.hh"
168 #include "G4CascadeParameters.hh"
169 #include "G4CollisionOutput.hh"
171 #include "G4Exp.hh"
172 #include "Randomize.hh"
173 #include "G4HadTmpUtil.hh"
174 #include "G4InuclNuclei.hh"
175 #include "G4InuclParticleNames.hh"
177 #include "G4Log.hh"
178 #include "G4LorentzConvertor.hh"
179 #include "G4Neutron.hh"
180 #include "G4ParticleDefinition.hh"
181 #include "G4ParticleLargerBeta.hh"
182 #include "G4PhysicalConstants.hh"
183 #include "G4Proton.hh"
184 #include "G4SystemOfUnits.hh"
185 #include <vector>
186 #include <algorithm>
187 #include <functional>
188 #include <numeric>
189 
190 using namespace G4InuclParticleNames;
191 using namespace G4InuclSpecialFunctions;
192 
193 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
194 
195 // Cut-off limits for extreme values in calculations
196 const G4double G4NucleiModel::small = 1.0e-9;
197 const G4double G4NucleiModel::large = 1000.;
198 
199 // For convenience in computing densities
201 
202 // Zone boundaries as fraction of nuclear radius (from outside in)
203 const G4double G4NucleiModel::alfa3[3] = { 0.7, 0.3, 0.01 };
204 const G4double G4NucleiModel::alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
205 
206 // Flat nuclear potentials for mesons and hyperons (GeV)
207 const G4double G4NucleiModel::pion_vp = 0.007;
209 const G4double G4NucleiModel::kaon_vp = 0.015;
210 const G4double G4NucleiModel::hyperon_vp = 0.030;
211 
212 namespace {
213  const G4double kebins[] =
214  { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
215  0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
216  2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
217 
218  const G4double gammaQDxsec[30] =
219  { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098,
220  0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005,
221  0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
222  0.0001, 0.0001, 0.0001 };
223 }
224 
225 
226 // Constructors and destructor
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()),
240  R_nucleon(G4CascadeParameters::radiusTrailing()),
241  gammaQDscale(G4CascadeParameters::gammaQDScale()),
242  neutronEP(neutron), protonEP(proton) {}
243 
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()),
257  R_nucleon(G4CascadeParameters::radiusTrailing()),
258  gammaQDscale(G4CascadeParameters::gammaQDScale()),
259  neutronEP(neutron), protonEP(proton) {
260  generateModel(a,z);
261 }
262 
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()),
276  R_nucleon(G4CascadeParameters::radiusTrailing()),
277  gammaQDscale(G4CascadeParameters::gammaQDScale()),
278  neutronEP(neutron), protonEP(proton) {
279  generateModel(nuclei);
280 }
281 
283  delete theNucleus;
284  theNucleus = 0;
285 }
286 
287 
288 // Initialize model state for new cascade
289 
290 void G4NucleiModel::reset(G4int nHitNeutrons, G4int nHitProtons,
291  const std::vector<G4ThreeVector>* hitPoints) {
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 }
299 
300 
301 // Generate nuclear model parameters for given nucleus
302 
304  generateModel(nuclei->getA(), nuclei->getZ());
305 }
306 
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 }
378 
379 
380 // Load binding energy array for current nucleus
381 
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 }
393 
394 // Load zone boundary radius array for current nucleus
395 
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 }
433 
434 // Compute zone-by-zone density integrals
435 
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 }
469 
470 // Load potentials for nucleons (using G4InuclParticle codes)
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 }
501 
502 // Zone integral of Woods-Saxon density function
504  G4double nuclearRadius) const {
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 }
553 
554 
555 // Zone integral of Gaussian density function
557  G4double nucRad) const {
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 }
602 
603 
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 }
624 
625 
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 }
637 
638 
641  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
643 
644  return generateWithRandomAngles(pmod, mass);
645 }
646 
647 
650  if (verboseLevel > 1) {
651  G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
652  }
653 
654  G4LorentzVector mom = generateNucleonMomentum(type, zone);
655  return G4InuclElementaryParticle(mom, type);
656 }
657 
658 
661  G4int zone) const {
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 }
683 
684 
685 void
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 }
854 
855 
856 void G4NucleiModel::
858  G4ElementaryParticleCollider* theEPCollider,
859  std::vector<G4CascadParticle>& outgoing_cparticles) {
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 
912  G4InuclElementaryParticle& target = thePartners[i].first;
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) {
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 =
940 
941  if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
942 
943  // Trailing effect: reject interaction at previously hit nucleon
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;
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;
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 }
1046 
1047 // Test if particle is suitable for absorption with dibaryons
1048 
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 }
1059 
1060 
1061 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles,
1062  G4int zone) {
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 }
1085 
1086 
1087 // Test here for trailing effect: loop over all previous collision
1088 // locations and test for d > R_nucleon
1089 
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 }
1105 
1106 
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 }
1170 
1171 
1172 // Select random point along full trajectory through nucleus
1173 // NOTE: Intended for projectile photons for initial interaction
1174 
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 }
1298 
1299 
1300 // Returns true if particle should interact immediately
1302  return (isProjectile(cparticle) &&
1303  (cparticle.getParticle().isPhoton() ||
1304  cparticle.getParticle().isMuon())
1305  );
1306 }
1307 
1309  return (cparticle.getGeneration() == 0); // Only initial state particles
1310 }
1311 
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 }
1341 
1342 
1344  if (verboseLevel > 4) {
1345  G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
1346  }
1347 
1348  switch (ip) {
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 }
1359 
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 }
1385 
1386 
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());
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 }
1412 
1414  G4InuclNuclei* target,
1415  modelLists& output) {
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());
1467  G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
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 }
1810 
1811 
1812 // Compute mean free path for inclusive interaction of projectile and target
1813 G4double
1815  const G4InuclElementaryParticle& target,
1816  G4int zone) {
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
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 }
1849 
1850 // Throw random distance for interaction of particle using path and MFP
1851 
1852 G4double
1854  G4double path, G4double invmfp) const {
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 }
1881 
1882 
1883 // Parametrized cross section for pion and photon absorption
1884 
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 }
1916 
1917 
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 }
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2, G4double nuclearRadius) const
void updateParticleMomentum(const G4LorentzVector &mom)
const G4double R_nucleon
G4InuclElementaryParticle generateQuasiDeuteron(G4int type1, G4int type2, G4int zone) const
std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
G4bool reflectedNow() const
static const G4CascadeChannel * GetTable(G4int initialState)
G4int getZ() const
G4LorentzVector getMomentum() const
void incrementReflectionCounter()
void incrementCurrentPath(G4double npath)
const G4double gammaQDscale
G4bool movingInsideNuclei() const
std::vector< G4double > vz
CLHEP::Hep3Vector G4ThreeVector
void generateInteractionPartners(G4CascadParticle &cparticle)
std::vector< G4double > zone_radii
G4int getGeneration() const
std::vector< G4double > zone_volumes
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4double z
Definition: TRTMaterials.hh:39
G4CascadeInterpolator< 30 > gammaQDinterp
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
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 v1[6]
G4double getPathToTheNextZone(G4double rz_in, G4double rz_out)
std::vector< partner > thePartners
G4double getEnergy() const
G4double nuclei_volume
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
G4double a
Definition: TRTMaterials.hh:39
void updatePosition(const G4ThreeVector &pos)
void boundaryTransition(G4CascadParticle &cparticle)
const G4InuclElementaryParticle protonEP
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
G4double ur[7]
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
#define G4ThreadLocal
Definition: tls.hh:89
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4CascadParticle initializeCascad(G4InuclElementaryParticle *particle)
G4InuclNuclei * theNucleus
G4bool young(G4double young_path_cut, G4double cpath) const
G4double getVolume(G4int izone) const
std::vector< G4double > pf
int G4int
Definition: G4Types.hh:78
const G4double radiusScale
G4double absorptionCrossSection(G4double e, G4int type) const
const G4InuclElementaryParticle & getParticle() const
static const G4double large
virtual void setVerboseLevel(G4int verbose=0)
virtual G4double getCrossSection(double ke) const =0
const G4double radiusScale2
G4bool isProjectile(const G4CascadParticle &cparticle) const
static const G4double piTimes4thirds
G4double getKineticEnergy() const
void fillBindingEnergies()
G4CollisionOutput EPCoutput
const G4double radScaleAlpha
void updateZone(G4int izone)
void Register(T *inst)
Definition: G4AutoDelete.hh:65
G4double getCurrentDensity(G4int ip, G4int izone) const
void fillZoneRadii(G4double nuclearRadius)
G4LorentzConvertor dummy_convertor
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4int getA() const
static const G4double kaon_vp
G4int protonNumberCurrent
void fillPotentials(G4int type, G4double tot_vol)
static const double deg
Definition: G4SIunits.hh:151
G4bool nucleon(G4int ityp)
static const G4double alfa6[6]
bool G4bool
Definition: G4Types.hh:79
std::vector< G4LorentzVector > momentums
G4double iz
Definition: TRTMaterials.hh:39
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
std::vector< G4ThreeVector > coordinates
G4int numberOfOutgoingParticles() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
std::vector< G4InuclElementaryParticle >::iterator particleIterator
G4double generateInteractionLength(const G4CascadParticle &cparticle, G4double path, G4double invmfp) const
G4double nuclei_radius
static const double second
Definition: G4SIunits.hh:156
std::vector< std::vector< G4double > > nucleon_densities
G4double fillZoneVolumes(G4double nuclearRadius)
static const double GeV
Definition: G4SIunits.hh:214
void generateModel(G4InuclNuclei *nuclei)
G4double v[6]
std::pair< G4InuclElementaryParticle, G4double > partner
static const G4double pion_vp_small
const G4double crossSectionUnits
G4double getKinEnergyInTheTRS() const
std::vector< G4double > binding_energies
const G4double radiusUnits
std::vector< G4InuclElementaryParticle > qdeutrons
static const G4double alfa3[3]
void generateParticleFate(G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
std::vector< G4double > rod
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double getRatio(G4int ip) const
G4double totalCrossSection(G4double ke, G4int rtype) const
static const G4double ab
static const double pi
Definition: G4SIunits.hh:74
static const G4double small
G4double getFermiMomentum(G4int ip, G4int izone) const
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
std::vector< G4InuclElementaryParticle > raw_particles
G4bool worthToPropagate(const G4CascadParticle &cparticle) const
static G4bool sortPartners(const partner &p1, const partner &p2)
static const G4double hyperon_vp
const G4double x[NPOINTSGL]
void choosePointAlongTraj(G4CascadParticle &cparticle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
G4int getCurrentZone() const
static const G4double b1
G4double getCharge() const
const G4InuclElementaryParticle neutronEP
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
std::vector< G4ThreeVector > collisionPts
void printModel() const
#define G4endl
Definition: G4ios.hh:61
std::vector< G4double > acsecs
G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const
const G4ThreeVector & getPosition() const
const G4double fermiMomentum
std::vector< std::vector< G4double > > fermi_momenta
double G4double
Definition: G4Types.hh:76
std::vector< std::vector< G4double > > zone_potentials
G4int neutronNumberCurrent
virtual ~G4NucleiModel()
static const G4double d2
G4double getDensity(G4int ip, G4int izone) const
G4double bindingEnergy(G4int A, G4int Z)
G4bool passTrailing(const G4ThreeVector &hit_position)
G4int getZone(G4double r) const
static const G4double pion_vp
const G4double radiusForSmall
const G4double skinDepth
double epsilon(double density, double temperature)
static const G4double pos
G4bool forceFirst(const G4CascadParticle &cparticle) const
G4GLOB_DLL std::ostream G4cerr
G4double getMass() const
void setTarget(const G4InuclParticle *target)
CLHEP::HepLorentzVector G4LorentzVector
G4double getPotential(G4int ip, G4int izone) const