Geant4  10.03
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< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
std::vector< G4double > zone_volumes
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4CascadeInterpolator< 30 > gammaQDinterp
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
void printCollisionOutput(std::ostream &os=G4cout) const
static const G4double d2
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
void updatePosition(const G4ThreeVector &pos)
static constexpr double second
Definition: G4SIunits.hh:157
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)
G4bool nucleon(G4int ityp)
static const G4double alfa6[6]
bool G4bool
Definition: G4Types.hh:79
std::vector< G4LorentzVector > momentums
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
std::vector< std::vector< G4double > > nucleon_densities
G4double fillZoneVolumes(G4double nuclearRadius)
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 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
void choosePointAlongTraj(G4CascadParticle &cparticle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
G4int getCurrentZone() const
static constexpr double GeV
Definition: G4SIunits.hh:217
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
static constexpr double pi
Definition: G4SIunits.hh:75
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
static constexpr double deg
Definition: G4SIunits.hh:152
std::vector< std::vector< G4double > > zone_potentials
G4int neutronNumberCurrent
virtual ~G4NucleiModel()
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