Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4NucleiModel.cc
이 파일의 문서화 페이지로 가기
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 // 20180227 A. Ribon -- Replaced obsolete std::bind2nd with std::bind
162 
163 #include "G4NucleiModel.hh"
164 #include "G4AutoDelete.hh"
165 #include "G4CascadeChannel.hh"
166 #include "G4CascadeChannelTables.hh"
167 #include "G4CascadeCheckBalance.hh"
168 #include "G4CascadeInterpolator.hh"
169 #include "G4CascadeParameters.hh"
170 #include "G4CollisionOutput.hh"
172 #include "G4Exp.hh"
173 #include "Randomize.hh"
174 #include "G4HadTmpUtil.hh"
175 #include "G4InuclNuclei.hh"
176 #include "G4InuclParticleNames.hh"
178 #include "G4Log.hh"
179 #include "G4LorentzConvertor.hh"
180 #include "G4Neutron.hh"
181 #include "G4ParticleDefinition.hh"
182 #include "G4ParticleLargerBeta.hh"
183 #include "G4PhysicalConstants.hh"
184 #include "G4Proton.hh"
185 #include "G4SystemOfUnits.hh"
186 #include <vector>
187 #include <algorithm>
188 #include <functional>
189 #include <numeric>
190 
191 using namespace G4InuclParticleNames;
192 using namespace G4InuclSpecialFunctions;
193 
194 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
195 
196 // Cut-off limits for extreme values in calculations
197 const G4double G4NucleiModel::small = 1.0e-9;
198 const G4double G4NucleiModel::large = 1000.;
199 
200 // For convenience in computing densities
202 
203 // Zone boundaries as fraction of nuclear radius (from outside in)
204 const G4double G4NucleiModel::alfa3[3] = { 0.7, 0.3, 0.01 };
205 const G4double G4NucleiModel::alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
206 
207 // Flat nuclear potentials for mesons and hyperons (GeV)
208 const G4double G4NucleiModel::pion_vp = 0.007;
210 const G4double G4NucleiModel::kaon_vp = 0.015;
211 const G4double G4NucleiModel::hyperon_vp = 0.030;
212 
213 namespace {
214  const G4double kebins[] =
215  { 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
216  0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
217  2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
218 
219  const G4double gammaQDxsec[30] =
220  { 2.8, 1.3, 0.89, 0.56, 0.38, 0.27, 0.19, 0.14, 0.098,
221  0.071, 0.054, 0.0003, 0.0007, 0.0027, 0.0014, 0.001, 0.0012, 0.0005,
222  0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
223  0.0001, 0.0001, 0.0001 };
224 }
225 
226 
227 // Constructors and destructor
229  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
230  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
231  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
232  current_nucl2(0), gammaQDinterp(kebins),
233  crossSectionUnits(G4CascadeParameters::xsecScale()),
234  radiusUnits(G4CascadeParameters::radiusScale()),
235  skinDepth(0.611207*radiusUnits),
236  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
237  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
238  radiusForSmall(G4CascadeParameters::radiusSmall()),
239  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
240  fermiMomentum(G4CascadeParameters::fermiScale()),
241  R_nucleon(G4CascadeParameters::radiusTrailing()),
242  gammaQDscale(G4CascadeParameters::gammaQDScale()),
243  neutronEP(neutron), protonEP(proton) {}
244 
246  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
247  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
248  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
249  current_nucl2(0), gammaQDinterp(kebins),
250  crossSectionUnits(G4CascadeParameters::xsecScale()),
251  radiusUnits(G4CascadeParameters::radiusScale()),
252  skinDepth(0.611207*radiusUnits),
253  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
254  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
255  radiusForSmall(G4CascadeParameters::radiusSmall()),
256  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
257  fermiMomentum(G4CascadeParameters::fermiScale()),
258  R_nucleon(G4CascadeParameters::radiusTrailing()),
259  gammaQDscale(G4CascadeParameters::gammaQDScale()),
260  neutronEP(neutron), protonEP(proton) {
261  generateModel(a,z);
262 }
263 
265  : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
266  A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
267  neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
268  current_nucl2(0), gammaQDinterp(kebins),
269  crossSectionUnits(G4CascadeParameters::xsecScale()),
270  radiusUnits(G4CascadeParameters::radiusScale()),
271  skinDepth(0.611207*radiusUnits),
272  radiusScale((G4CascadeParameters::useTwoParam()?1.16:1.2)*radiusUnits),
273  radiusScale2((G4CascadeParameters::useTwoParam()?-1.3456:0.)*radiusUnits),
274  radiusForSmall(G4CascadeParameters::radiusSmall()),
275  radScaleAlpha(G4CascadeParameters::radiusAlpha()),
276  fermiMomentum(G4CascadeParameters::fermiScale()),
277  R_nucleon(G4CascadeParameters::radiusTrailing()),
278  gammaQDscale(G4CascadeParameters::gammaQDScale()),
279  neutronEP(neutron), protonEP(proton) {
280  generateModel(nuclei);
281 }
282 
284  delete theNucleus;
285  theNucleus = 0;
286 }
287 
288 
289 // Initialize model state for new cascade
290 
291 void G4NucleiModel::reset(G4int nHitNeutrons, G4int nHitProtons,
292  const std::vector<G4ThreeVector>* hitPoints) {
293  neutronNumberCurrent = neutronNumber - nHitNeutrons;
294  protonNumberCurrent = protonNumber - nHitProtons;
295 
296  // zero or copy collision point array for trailing effect
297  if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
298  else collisionPts = *hitPoints;
299 }
300 
301 
302 // Generate nuclear model parameters for given nucleus
303 
305  generateModel(nuclei->getA(), nuclei->getZ());
306 }
307 
309  if (verboseLevel) {
310  G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
311  << G4endl;
312  }
313 
314  // If model already built, just return; otherwise intialize everything
315  if (a == A && z == Z) {
316  if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
317  reset();
318  return;
319  }
320 
321  A = a;
322  Z = z;
323  delete theNucleus;
324  theNucleus = new G4InuclNuclei(A,Z); // For conservation checking
325 
326  neutronNumber = A - Z;
327  protonNumber = Z;
328  reset();
329 
330  if (verboseLevel > 3) {
331  G4cout << " crossSectionUnits = " << crossSectionUnits << G4endl
332  << " radiusUnits = " << radiusUnits << G4endl
333  << " skinDepth = " << skinDepth << G4endl
334  << " radiusScale = " << radiusScale << G4endl
335  << " radiusScale2 = " << radiusScale2 << G4endl
336  << " radiusForSmall = " << radiusForSmall << G4endl
337  << " radScaleAlpha = " << radScaleAlpha << G4endl
338  << " fermiMomentum = " << fermiMomentum << G4endl
339  << " piTimes4thirds = " << piTimes4thirds << G4endl;
340  }
341 
342  G4double nuclearRadius; // Nuclear radius computed from A
343  if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
344  else nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
345 
346  // This will be used to pre-allocate lots of arrays below
347  number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
348 
349  // Clear all parameters arrays for reloading
350  binding_energies.clear();
351  nucleon_densities.clear();
352  zone_potentials.clear();
353  fermi_momenta.clear();
354  zone_radii.clear();
355  zone_volumes.clear();
356 
358  fillZoneRadii(nuclearRadius);
359 
360  G4double tot_vol = fillZoneVolumes(nuclearRadius); // Woods-Saxon integral
361 
362  fillPotentials(proton, tot_vol); // Protons
363  fillPotentials(neutron, tot_vol); // Neutrons
364 
365  // Additional flat zone potentials for other hadrons
366  const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
367  const std::vector<G4double> kp(number_of_zones, kaon_vp);
368  const std::vector<G4double> hp(number_of_zones, hyperon_vp);
369 
370  zone_potentials.push_back(vp);
371  zone_potentials.push_back(kp);
372  zone_potentials.push_back(hp);
373 
374  nuclei_radius = zone_radii.back();
375  nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
376 
377  if (verboseLevel > 3) printModel();
378 }
379 
380 
381 // Load binding energy array for current nucleus
382 
384  if (verboseLevel > 1)
385  G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
386 
387  G4double dm = bindingEnergy(A,Z);
388 
389  // Binding energy differences for proton and neutron loss, respectively
390  // FIXME: Why is fabs() used here instead of the signed difference?
391  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
392  binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
393 }
394 
395 // Load zone boundary radius array for current nucleus
396 
398  if (verboseLevel > 1)
399  G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
400 
401  G4double skinRatio = nuclearRadius/skinDepth;
402  G4double skinDecay = G4Exp(-skinRatio);
403 
404  if (A < 5) { // Light ions treated as simple balls
405  zone_radii.push_back(nuclearRadius);
406  ur[0] = 0.;
407  ur[1] = 1.;
408  } else if (A < 12) { // Small nuclei have Gaussian potential
409  G4double rSq = nuclearRadius * nuclearRadius;
410  G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
411 
412  ur[0] = 0.0;
413  for (G4int i = 0; i < number_of_zones; i++) {
414  G4double y = std::sqrt(-G4Log(alfa3[i]));
415  zone_radii.push_back(gaussRadius * y);
416  ur[i+1] = y;
417  }
418  } else if (A < 100) { // Intermediate nuclei have three-zone W-S
419  ur[0] = -skinRatio;
420  for (G4int i = 0; i < number_of_zones; i++) {
421  G4double y = G4Log((1.0 + skinDecay)/alfa3[i] - 1.0);
422  zone_radii.push_back(nuclearRadius + skinDepth * y);
423  ur[i+1] = y;
424  }
425  } else { // Heavy nuclei have six-zone W-S
426  ur[0] = -skinRatio;
427  for (G4int i = 0; i < number_of_zones; i++) {
428  G4double y = G4Log((1.0 + skinDecay)/alfa6[i] - 1.0);
429  zone_radii.push_back(nuclearRadius + skinDepth * y);
430  ur[i+1] = y;
431  }
432  }
433 }
434 
435 // Compute zone-by-zone density integrals
436 
438  if (verboseLevel > 1)
439  G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
440 
441  G4double tot_vol = 0.; // Return value omitting 4pi/3 for sphere!
442 
443  if (A < 5) { // Light ions are treated as simple balls
444  v[0] = v1[0] = 1.;
445  tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
446  zone_volumes.push_back(tot_vol*piTimes4thirds); // True volume of sphere
447 
448  return tot_vol;
449  }
450 
451  PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
452 
453  for (G4int i = 0; i < number_of_zones; i++) {
454  if (usePotential == WoodsSaxon) {
455  v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
456  } else {
457  v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
458  }
459 
460  tot_vol += v[i];
461 
462  v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
463  if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
464 
465  zone_volumes.push_back(v1[i]*piTimes4thirds); // True volume of shell
466  }
467 
468  return tot_vol; // Sum of zone integrals, not geometric volume
469 }
470 
471 // Load potentials for nucleons (using G4InuclParticle codes)
473  if (verboseLevel > 1)
474  G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
475 
476  if (type != proton && type != neutron) return;
477 
479 
480  // FIXME: This is the fabs() binding energy difference, not signed
481  const G4double dm = binding_energies[type-1];
482 
483  rod.clear(); rod.reserve(number_of_zones);
484  pf.clear(); pf.reserve(number_of_zones);
485  vz.clear(); vz.reserve(number_of_zones);
486 
487  G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
488  G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
489 
490  for (G4int i = 0; i < number_of_zones; i++) {
491  G4double rd = dd0 * v[i] / v1[i];
492  rod.push_back(rd);
493  G4double pff = fermiMomentum * G4cbrt(rd);
494  pf.push_back(pff);
495  vz.push_back(0.5 * pff * pff / mass + dm);
496  }
497 
498  nucleon_densities.push_back(rod);
499  fermi_momenta.push_back(pf);
500  zone_potentials.push_back(vz);
501 }
502 
503 // Zone integral of Woods-Saxon density function
505  G4double nuclearRadius) const {
506  if (verboseLevel > 1) {
507  G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
508  }
509 
510  const G4double epsilon = 1.0e-3;
511  const G4int itry_max = 1000;
512 
513  G4double skinRatio = nuclearRadius / skinDepth;
514 
515  G4double d2 = 2.0 * skinRatio;
516  G4double dr = r2 - r1;
517  G4double fr1 = r1 * (r1 + d2) / (1.0 + G4Exp(r1));
518  G4double fr2 = r2 * (r2 + d2) / (1.0 + G4Exp(r2));
519  G4double fi = (fr1 + fr2) / 2.;
520  G4double fun1 = fi * dr;
521  G4double fun;
522  G4int jc = 1;
523  G4double dr1 = dr;
524  G4int itry = 0;
525 
526  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
527  dr /= 2.;
528  itry++;
529 
530  G4double r = r1 - dr;
531  fi = 0.0;
532 
533  for (G4int i = 0; i < jc; i++) {
534  r += dr1;
535  fi += r * (r + d2) / (1.0 + G4Exp(r));
536  }
537 
538  fun = 0.5 * fun1 + fi * dr;
539 
540  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
541 
542  jc *= 2;
543  dr1 = dr;
544  fun1 = fun;
545  } // while (itry < itry_max)
546 
547  if (verboseLevel > 2 && itry == itry_max)
548  G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
549 
550  G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
551 
552  return skinDepth3 * (fun + skinRatio*skinRatio*G4Log((1.0 + G4Exp(-r1)) / (1.0 + G4Exp(-r2))));
553 }
554 
555 
556 // Zone integral of Gaussian density function
558  G4double nucRad) const {
559  if (verboseLevel > 1) {
560  G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
561  }
562 
563  G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
564 
565  const G4double epsilon = 1.0e-3;
566  const G4int itry_max = 1000;
567 
568  G4double dr = r2 - r1;
569  G4double fr1 = r1 * r1 * G4Exp(-r1 * r1);
570  G4double fr2 = r2 * r2 * G4Exp(-r2 * r2);
571  G4double fi = (fr1 + fr2) / 2.;
572  G4double fun1 = fi * dr;
573  G4double fun;
574  G4int jc = 1;
575  G4double dr1 = dr;
576  G4int itry = 0;
577 
578  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
579  dr /= 2.;
580  itry++;
581  G4double r = r1 - dr;
582  fi = 0.0;
583 
584  for (G4int i = 0; i < jc; i++) {
585  r += dr1;
586  fi += r * r * G4Exp(-r * r);
587  }
588 
589  fun = 0.5 * fun1 + fi * dr;
590 
591  if (std::fabs((fun - fun1) / fun) <= epsilon) break;
592 
593  jc *= 2;
594  dr1 = dr;
595  fun1 = fun;
596  } // while (itry < itry_max)
597 
598  if (verboseLevel > 2 && itry == itry_max)
599  G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
600 
601  return gaussRadius*gaussRadius*gaussRadius * fun;
602 }
603 
604 
606  if (verboseLevel > 1) {
607  G4cout << " >>> G4NucleiModel::printModel" << G4endl;
608  }
609 
610  G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
611  << " proton binding energy " << binding_energies[0]
612  << " neutron binding energy " << binding_energies[1] << G4endl
613  << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
614  << " number of zones " << number_of_zones << G4endl;
615 
616  for (G4int i = 0; i < number_of_zones; i++)
617  G4cout << " zone " << i+1 << " radius " << zone_radii[i]
618  << " volume " << zone_volumes[i] << G4endl
619  << " protons: density " << getDensity(1,i) << " PF " <<
620  getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
621  << " neutrons: density " << getDensity(2,i) << " PF " <<
622  getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
623  << " pions: VP " << getPotential(3,i) << G4endl;
624 }
625 
626 
628  G4double ekin = 0.0;
629 
630  if (ip < 3 && izone < number_of_zones) { // ip for proton/neutron only
631  G4double pfermi = fermi_momenta[ip - 1][izone];
633 
634  ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
635  }
636  return ekin;
637 }
638 
639 
642  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
644 
645  return generateWithRandomAngles(pmod, mass);
646 }
647 
648 
651  if (verboseLevel > 1) {
652  G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
653  }
654 
655  G4LorentzVector mom = generateNucleonMomentum(type, zone);
656  return G4InuclElementaryParticle(mom, type);
657 }
658 
659 
662  G4int zone) const {
663  if (verboseLevel > 1) {
664  G4cout << " >>> G4NucleiModel::generateQuasiDeuteron" << G4endl;
665  }
666 
667  // Quasideuteron consists of an unbound but associated nucleon pair
668 
669  // FIXME: Why generate two separate nucleon momenta (randomly!) and
670  // add them, instead of just throwing a net momentum for the
671  // dinulceon state? And why do I have to capture the two
672  // return values into local variables?
673  G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
674  G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
675  G4LorentzVector dmom = mom1+mom2;
676 
677  G4int dtype = 0;
678  if (type1*type2 == pro*pro) dtype = 111;
679  else if (type1*type2 == pro*neu) dtype = 112;
680  else if (type1*type2 == neu*neu) dtype = 122;
681 
682  return G4InuclElementaryParticle(dmom, dtype);
683 }
684 
685 
686 void
688  if (verboseLevel > 1) {
689  G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
690  }
691 
692  thePartners.clear(); // Reset buffer for next cycle
693 
694  G4int ptype = cparticle.getParticle().type();
695  G4int zone = cparticle.getCurrentZone();
696 
697  G4double r_in; // Destination radius if inbound
698  G4double r_out; // Destination radius if outbound
699 
700  if (zone == number_of_zones) {
701  r_in = nuclei_radius;
702  r_out = 0.0;
703  } else if (zone == 0) { // particle is outside core
704  r_in = 0.0;
705  r_out = zone_radii[0];
706  } else {
707  r_in = zone_radii[zone - 1];
708  r_out = zone_radii[zone];
709  }
710 
711  G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
712 
713  if (verboseLevel > 2) {
714  if (isProjectile(cparticle)) G4cout << " incident particle: ";
715  G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
716  << G4endl;
717  }
718 
719  if (path < -small) { // something wrong
720  if (verboseLevel)
721  G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
722  return;
723  }
724 
725  if (std::fabs(path) < small) { // Not moving, or just at boundary
726  if (cparticle.getMomentum().vect().mag() > small) {
727  if (verboseLevel > 3)
728  G4cout << " generateInteractionPartners-> zero path" << G4endl;
729 
730  thePartners.push_back(partner()); // Dummy list terminator with zero path
731  return;
732  }
733 
734  if (zone >= number_of_zones) // Place captured-at-rest in nucleus
735  zone = number_of_zones-1;
736  }
737 
738  G4double invmfp = 0.; // Buffers for interaction probability
739  G4double spath = 0.;
740  for (G4int ip = 1; ip < 3; ip++) {
741  // Only process nucleons which remain active in target
742  if (ip==proton && protonNumberCurrent < 1) continue;
743  if (ip==neutron && neutronNumberCurrent < 1) continue;
744  if (ip==neutron && ptype==muonMinus) continue; // mu-/n forbidden
745 
746  // All nucleons are assumed to be at rest when colliding
747  G4InuclElementaryParticle particle = generateNucleon(ip, zone);
748  invmfp = inverseMeanFreePath(cparticle, particle);
749  spath = generateInteractionLength(cparticle, path, invmfp);
750 
751  if (path<small || spath < path) {
752  if (verboseLevel > 3) {
753  G4cout << " adding partner[" << thePartners.size() << "]: "
754  << particle << G4endl;
755  }
756  thePartners.push_back(partner(particle, spath));
757  }
758  } // for (G4int ip...
759 
760  if (verboseLevel > 2) {
761  G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
762  }
763 
764  // Absorption possible for pions or photons interacting with dibaryons
765  if (useQuasiDeuteron(cparticle.getParticle().type())) {
766  if (verboseLevel > 2) {
767  G4cout << " trying quasi-deuterons with bullet: "
768  << cparticle.getParticle() << G4endl;
769  }
770 
771  // Initialize buffers for quasi-deuteron results
772  qdeutrons.clear();
773  acsecs.clear();
774 
775  G4double tot_invmfp = 0.0; // Total inv. mean-free-path for all QDs
776 
777  // Proton-proton state interacts with pi-, mu- or neutrals
778  if (protonNumberCurrent >= 2 && ptype != pip) {
780  if (verboseLevel > 2)
781  G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
782 
783  invmfp = inverseMeanFreePath(cparticle, ppd);
784  tot_invmfp += invmfp;
785  acsecs.push_back(invmfp);
786  qdeutrons.push_back(ppd);
787  }
788 
789  // Proton-neutron state interacts with any pion type or photon
790  if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
792  if (verboseLevel > 2)
793  G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
794 
795  invmfp = inverseMeanFreePath(cparticle, npd);
796  tot_invmfp += invmfp;
797  acsecs.push_back(invmfp);
798  qdeutrons.push_back(npd);
799  }
800 
801  // Neutron-neutron state interacts with pi+ or neutrals
802  if (neutronNumberCurrent >= 2 && ptype != pim && ptype != mum) {
804  if (verboseLevel > 2)
805  G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
806 
807  invmfp = inverseMeanFreePath(cparticle, nnd);
808  tot_invmfp += invmfp;
809  acsecs.push_back(invmfp);
810  qdeutrons.push_back(nnd);
811  }
812 
813  // Select quasideuteron interaction from non-zero cross-section choices
814  if (verboseLevel > 2) {
815  for (size_t i=0; i<qdeutrons.size(); i++) {
816  G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
817  << "] " << acsecs[i];
818  }
819  G4cout << G4endl;
820  }
821 
822  if (tot_invmfp > small) { // Must have absorption cross-section
823  G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
824 
825  if (path<small || apath < path) { // choose the qdeutron
826  G4double sl = inuclRndm() * tot_invmfp;
827  G4double as = 0.0;
828 
829  for (size_t i = 0; i < qdeutrons.size(); i++) {
830  as += acsecs[i];
831  if (sl < as) {
832  if (verboseLevel > 2)
833  G4cout << " deut type " << qdeutrons[i] << G4endl;
834 
835  thePartners.push_back(partner(qdeutrons[i], apath));
836  break;
837  }
838  } // for (qdeutrons...
839  } // if (apath < path)
840  } // if (tot_invmfp > small)
841  } // if (useQuasiDeuteron) [pion, muon or photon]
842 
843  if (verboseLevel > 2) {
844  G4cout << " after deuterons " << thePartners.size() << " partners"
845  << G4endl;
846  }
847 
848  if (thePartners.size() > 1) { // Sort list by path length
849  std::sort(thePartners.begin(), thePartners.end(), sortPartners);
850  }
851 
852  G4InuclElementaryParticle particle; // Total path at end of list
853  thePartners.push_back(partner(particle, path));
854 }
855 
856 
857 void G4NucleiModel::
859  G4ElementaryParticleCollider* theEPCollider,
860  std::vector<G4CascadParticle>& outgoing_cparticles) {
861  if (verboseLevel > 1)
862  G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
863 
864  if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
865 
866  // Create four-vector checking
867 #ifdef G4CASCADE_CHECK_ECONS
868  G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel"); // Second arg is in GeV
869  balance.setVerboseLevel(verboseLevel);
870 #endif
871 
872  outgoing_cparticles.clear(); // Clear return buffer for this event
873  generateInteractionPartners(cparticle); // Fills "thePartners" data
874 
875  if (thePartners.empty()) { // smth. is wrong -> needs special treatment
876  if (verboseLevel)
877  G4cerr << " generateParticleFate-> got empty interaction-partners list "
878  << G4endl;
879  return;
880  }
881 
882  G4int npart = thePartners.size(); // Last item is a total-path placeholder
883 
884  if (npart == 1) { // cparticle is on the next zone entry
885  if (verboseLevel > 1)
886  G4cout << " no interactions; moving to next zone" << G4endl;
887 
889  cparticle.incrementCurrentPath(thePartners[0].second);
890  boundaryTransition(cparticle);
891  outgoing_cparticles.push_back(cparticle);
892 
893  if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
894 
895  //A.R. 19-Jun-2013: Fixed rare cases of non-reproducibility.
896  current_nucl1 = 0;
897  current_nucl2 = 0;
898 
899  return;
900  } // if (npart == 1)
901 
902  if (verboseLevel > 1)
903  G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
904 
905  G4ThreeVector old_position = cparticle.getPosition();
906  G4InuclElementaryParticle& bullet = cparticle.getParticle();
907  G4bool no_interaction = true;
908  G4int zone = cparticle.getCurrentZone();
909 
910  for (G4int i=0; i<npart-1; i++) { // Last item is a total-path placeholder
911  if (i > 0) cparticle.updatePosition(old_position);
912 
914 
915  if (verboseLevel > 3) {
916  if (target.quasi_deutron()) G4cout << " try absorption: ";
917  G4cout << " target " << target.type() << " bullet " << bullet.type()
918  << G4endl;
919  }
920 
921  EPCoutput.reset();
922  // Pass current (A,Z) configuration for possible recoils
923  G4int massNumberCurrent = protonNumberCurrent+neutronNumberCurrent;
924  theEPCollider->setNucleusState(massNumberCurrent, protonNumberCurrent);
925  theEPCollider->collide(&bullet, &target, EPCoutput);
926 
927  // If collision failed, exit loop over partners
928  if (EPCoutput.numberOfOutgoingParticles() == 0) break;
929 
930  if (verboseLevel > 2) {
932 #ifdef G4CASCADE_CHECK_ECONS
933  balance.collide(&bullet, &target, EPCoutput);
934  balance.okay(); // Do checks, but ignore result
935 #endif
936  }
937 
938  // Get list of outgoing particles for evaluation
939  std::vector<G4InuclElementaryParticle>& outgoing_particles =
941 
942  if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
943 
944  // Trailing effect: reject interaction at previously hit nucleon
946  const G4ThreeVector& new_position = cparticle.getPosition();
947 
948  if (!passTrailing(new_position)) continue;
949  collisionPts.push_back(new_position);
950 
951  // Sort particles according to beta (fastest first)
952  std::sort(outgoing_particles.begin(), outgoing_particles.end(),
954 
955  if (verboseLevel > 2)
956  G4cout << " adding " << outgoing_particles.size()
957  << " output particles" << G4endl;
958 
959  // NOTE: Embedded temporary is optimized away (no copying gets done)
960  G4int nextGen = cparticle.getGeneration()+1;
961  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) {
962  outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
963  new_position, zone,
964  0.0, nextGen));
965  }
966 
967  no_interaction = false;
968  current_nucl1 = 0;
969  current_nucl2 = 0;
970 
971 #ifdef G4CASCADE_DEBUG_CHARGE
972  {
973  G4double out_charge = 0.0;
974 
975  for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
976  out_charge += outgoing_particles[ip].getCharge();
977 
978  G4cout << " multiplicity " << outgoing_particles.size()
979  << " bul type " << bullet.type()
980  << " targ type " << target.type()
981  << "\n initial charge "
982  << bullet.getCharge() + target.getCharge()
983  << " out charge " << out_charge << G4endl;
984  }
985 #endif
986 
987  if (verboseLevel > 2)
988  G4cout << " partner type " << target.type() << G4endl;
989 
990  if (target.nucleon()) {
991  current_nucl1 = target.type();
992  } else {
993  if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
994  current_nucl1 = (target.type() - 100) / 10;
995  current_nucl2 = target.type() - 100 - 10 * current_nucl1;
996  }
997 
998  if (current_nucl1 == 1) {
999  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
1001  } else {
1002  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
1003  neutronNumberCurrent--;
1004  }
1005 
1006  if (current_nucl2 == 1) {
1007  if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
1009  } else if (current_nucl2 == 2) {
1010  if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
1011  neutronNumberCurrent--;
1012  }
1013 
1014  break;
1015  } // loop over partners
1016 
1017  if (no_interaction) { // still no interactions
1018  if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
1019 
1020  // For conservation checking (below), get particle before updating
1021  static G4ThreadLocal G4InuclElementaryParticle *prescatCP_G4MT_TLS_ = 0;
1022  if (!prescatCP_G4MT_TLS_) {
1023  prescatCP_G4MT_TLS_ = new G4InuclElementaryParticle;
1024  G4AutoDelete::Register(prescatCP_G4MT_TLS_);
1025  }
1026  G4InuclElementaryParticle &prescatCP = *prescatCP_G4MT_TLS_; // Avoid memory churn
1027  prescatCP = cparticle.getParticle();
1028 
1029  // Last "partner" is just a total-path placeholder
1030  cparticle.updatePosition(old_position);
1031  cparticle.propagateAlongThePath(thePartners[npart-1].second);
1032  cparticle.incrementCurrentPath(thePartners[npart-1].second);
1033  boundaryTransition(cparticle);
1034  outgoing_cparticles.push_back(cparticle);
1035 
1036  // Check conservation for simple scattering (ignore target nucleus!)
1037 #ifdef G4CASCADE_CHECK_ECONS
1038  if (verboseLevel > 2) {
1039  balance.collide(&prescatCP, 0, outgoing_cparticles);
1040  balance.okay(); // Report violations, but don't act on them
1041  }
1042 #endif
1043  } // if (no_interaction)
1044 
1045  return;
1046 }
1047 
1048 // Test if particle is suitable for absorption with dibaryons
1049 
1051  if (qdtype==pn || qdtype==0) // All absorptive particles
1052  return (ptype==pi0 || ptype==pip || ptype==pim || ptype==gam || ptype==mum);
1053  else if (qdtype==pp) // Negative or neutral only
1054  return (ptype==pi0 || ptype==pim || ptype==gam || ptype==mum);
1055  else if (qdtype==nn) // Positive or neutral only
1056  return (ptype==pi0 || ptype==pip || ptype==gam);
1057 
1058  return false; // Not a quasideuteron target
1059 }
1060 
1061 
1062 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles,
1063  G4int zone) {
1064  if (verboseLevel > 1) {
1065  G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
1066  }
1067 
1068  // Only check Fermi momenta for nucleons
1069  for (G4int i = 0; i < G4int(particles.size()); i++) {
1070  if (!particles[i].nucleon()) continue;
1071 
1072  G4int type = particles[i].type();
1073  G4double mom = particles[i].getMomModule();
1074  G4double pfermi = fermi_momenta[type-1][zone];
1075 
1076  if (verboseLevel > 2)
1077  G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
1078 
1079  if (mom < pfermi) {
1080  if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
1081  return false;
1082  }
1083  }
1084  return true;
1085 }
1086 
1087 
1088 // Test here for trailing effect: loop over all previous collision
1089 // locations and test for d > R_nucleon
1090 
1092  if (verboseLevel > 1)
1093  G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
1094 
1095  G4double dist;
1096  for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
1097  dist = (collisionPts[i] - hit_position).mag();
1098  if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
1099  if (dist < R_nucleon) {
1100  if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
1101  return false;
1102  }
1103  }
1104  return true; // New point far enough away to be used
1105 }
1106 
1107 
1109  if (verboseLevel > 1) {
1110  G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
1111  }
1112 
1113  G4int zone = cparticle.getCurrentZone();
1114 
1115  if (cparticle.movingInsideNuclei() && zone == 0) {
1116  if (verboseLevel) G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
1117  return;
1118  }
1119 
1120  G4LorentzVector mom = cparticle.getMomentum();
1121  G4ThreeVector pos = cparticle.getPosition();
1122 
1123  G4int type = cparticle.getParticle().type();
1124 
1125  G4double r = pos.mag();
1126  G4double pr = pos.dot(mom.vect()) / r;
1127 
1128  G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
1129 
1130  // NOTE: dv is the height of the wall seen by the particle
1131  G4double dv = getPotential(type,next_zone) - getPotential(type,zone);
1132  if (verboseLevel > 3) {
1133  G4cout << "Potentials for type " << type << " = "
1134  << getPotential(type,zone) << " , "
1135  << getPotential(type,next_zone) << G4endl;
1136  }
1137 
1138  G4double qv = dv * dv + 2.0 * dv * mom.e() + pr * pr;
1139 
1140  G4double p1r = 0.;
1141 
1142  if (verboseLevel > 3) {
1143  G4cout << " type " << type << " zone " << zone << " next " << next_zone
1144  << " qv " << qv << " dv " << dv << G4endl;
1145  }
1146 
1147  if (qv <= 0.0) { // reflection
1148  if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
1149  p1r = -pr;
1150  cparticle.incrementReflectionCounter();
1151  } else { // transition
1152  if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
1153  p1r = std::sqrt(qv);
1154  if (pr < 0.0) p1r = -p1r;
1155 
1156  cparticle.updateZone(next_zone);
1157  cparticle.resetReflection();
1158  }
1159 
1160  G4double prr = (p1r - pr) / r;
1161 
1162  if (verboseLevel > 3) {
1163  G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
1164  << prr*pos.y() << " pz " << prr*pos.z() << " mag "
1165  << std::fabs(prr*r) << G4endl;
1166  }
1167 
1168  mom.setVect(mom.vect() + pos*prr);
1169  cparticle.updateParticleMomentum(mom);
1170 }
1171 
1172 
1173 // Select random point along full trajectory through nucleus
1174 // NOTE: Intended for projectile photons for initial interaction
1175 
1177  if (verboseLevel > 1)
1178  G4cout << " >>> G4NucleiModel::choosePointAlongTraj" << G4endl;
1179 
1180  // Get trajectory through nucleus by computing exit point of line,
1181  // assuming that current position is on surface
1182 
1183  // FIXME: These need to be reusable (static) buffers
1184  G4ThreeVector pos = cparticle.getPosition();
1185  G4ThreeVector rhat = pos.unit();
1186 
1187  G4ThreeVector phat = cparticle.getMomentum().vect().unit();
1188  if (cparticle.getMomentum().vect().mag() < small) phat.set(0.,0.,1.);
1189 
1190  if (verboseLevel > 3)
1191  G4cout << " pos " << pos << " phat " << phat << " rhat " << rhat << G4endl;
1192 
1193  G4ThreeVector posout = pos;
1194  G4double prang = rhat.angle(-phat);
1195 
1196  if (prang < 1e-6) posout = -pos; // Check for radial incidence
1197  else {
1198  G4double posrot = 2.*prang - pi;
1199  posout.rotate(posrot, phat.cross(rhat));
1200  if (verboseLevel > 3) G4cout << " posrot " << posrot/deg << " deg";
1201  }
1202 
1203  if (verboseLevel > 3) G4cout << " posout " << posout << G4endl;
1204 
1205  // Get list of zone crossings along trajectory
1206  G4ThreeVector posmid = (pos+posout)/2.; // Midpoint of trajectory
1207  G4double r2mid = posmid.mag2();
1208  G4double lenmid = (posout-pos).mag()/2.; // Half-length of trajectory
1209 
1210  G4int zoneout = number_of_zones-1;
1211  G4int zonemid = getZone(posmid.mag()); // Middle zone traversed
1212 
1213  // Every zone is entered then exited, so preallocate vector
1214  G4int ncross = (number_of_zones-zonemid)*2;
1215 
1216  if (verboseLevel > 3) {
1217  G4cout << " posmid " << posmid << " lenmid " << lenmid
1218  << " zoneout " << zoneout << " zonemid " << zonemid
1219  << " ncross " << ncross << G4endl;
1220  }
1221 
1222  // FIXME: These need to be reusable (static) buffers
1223  std::vector<G4double> wtlen(ncross,0.); // CDF from entry point
1224  std::vector<G4double> len(ncross,0.); // Distance from entry point
1225 
1226  // Work from outside in, to accumulate CDF steps properly
1227  G4int i; // Loop variable, used multiple times
1228  for (i=0; i<ncross/2; i++) {
1229  G4int iz = zoneout-i;
1230  G4double ds = std::sqrt(zone_radii[iz]*zone_radii[iz]-r2mid);
1231 
1232  len[i] = lenmid - ds; // Distance from entry to crossing
1233  len[ncross-1-i] = lenmid + ds; // Distance to outbound crossing
1234 
1235  if (verboseLevel > 3) {
1236  G4cout << " i " << i << " iz " << iz << " ds " << ds
1237  << " len " << len[i] << G4endl;
1238  }
1239  }
1240 
1241  // Compute weights for each zone along trajectory and integrate
1242  for (i=1; i<ncross; i++) {
1243  G4int iz = (i<ncross/2) ? zoneout-i+1 : zoneout-ncross+i+1;
1244 
1245  G4double dlen = len[i]-len[i-1]; // Distance from previous crossing
1246 
1247  // Uniform probability across entire zone
1248  //*** G4double wt = dlen*(getDensity(1,iz)+getDensity(2,iz));
1249 
1250  // Probability based on interaction length through zone
1251  G4double invmfp = (inverseMeanFreePath(cparticle, neutronEP, iz)
1252  + inverseMeanFreePath(cparticle, protonEP, iz));
1253 
1254  // Integral of exp(-len/mfp) from start of zone to end
1255  G4double wt = (G4Exp(-len[i-1]*invmfp)-G4Exp(-len[i]*invmfp)) / invmfp;
1256 
1257  wtlen[i] = wtlen[i-1] + wt;
1258 
1259  if (verboseLevel > 3) {
1260  G4cout << " i " << i << " iz " << iz << " avg.mfp " << 1./invmfp
1261  << " dlen " << dlen << " wt " << wt << " wtlen " << wtlen[i]
1262  << G4endl;
1263  }
1264  }
1265 
1266  // Normalize CDF to unit integral
1267  std::transform(wtlen.begin(), wtlen.end(), wtlen.begin(),
1268  std::bind(std::divides<G4double>(), std::placeholders::_1, wtlen.back()));
1269 
1270  if (verboseLevel > 3) {
1271  G4cout << " weights";
1272  for (i=0; i<ncross; i++) G4cout << " " << wtlen[i];
1273  G4cout << G4endl;
1274  }
1275 
1276  // Choose random point along trajectory, weighted by density
1277  G4double rand = G4UniformRand();
1278  G4int ir = std::upper_bound(wtlen.begin(),wtlen.end(),rand) - wtlen.begin();
1279 
1280  G4double frac = (rand-wtlen[ir-1]) / (wtlen[ir]-wtlen[ir-1]);
1281  G4double drand = (1.-frac)*len[ir-1] + frac*len[ir];
1282 
1283  if (verboseLevel > 3) {
1284  G4cout << " rand " << rand << " ir " << ir << " frac " << frac
1285  << " drand " << drand << G4endl;
1286  }
1287 
1288  pos += drand * phat; // Distance from entry point along trajectory
1289 
1290  // Update input particle with new location
1291  cparticle.updatePosition(pos);
1292  cparticle.updateZone(getZone(pos.mag()));
1293 
1294  if (verboseLevel > 2) {
1295  G4cout << " moved particle to zone " << cparticle.getCurrentZone()
1296  << " @ " << pos << G4endl;
1297  }
1298 }
1299 
1300 
1301 // Returns true if particle should interact immediately
1303  return (isProjectile(cparticle) &&
1304  (cparticle.getParticle().isPhoton() ||
1305  cparticle.getParticle().isMuon())
1306  );
1307 }
1308 
1310  return (cparticle.getGeneration() == 0); // Only initial state particles
1311 }
1312 
1314  if (verboseLevel > 1) {
1315  G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
1316  }
1317 
1318  const G4double ekin_scale = 2.0;
1319 
1320  G4bool worth = true;
1321 
1322  if (cparticle.reflectedNow()) { // Just reflected -- keep going?
1323  G4int zone = cparticle.getCurrentZone();
1324  G4int ip = cparticle.getParticle().type();
1325 
1326  // NOTE: Temporarily backing out use of potential for non-nucleons
1327  G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
1328  getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
1329 
1330  worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
1331 
1332  if (verboseLevel > 3) {
1333  G4cout << " type=" << ip
1334  << " ekin=" << cparticle.getParticle().getKineticEnergy()
1335  << " potential=" << ekin_cut
1336  << " : worth? " << worth << G4endl;
1337  }
1338  }
1339 
1340  return worth;
1341 }
1342 
1343 
1345  if (verboseLevel > 4) {
1346  G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
1347  }
1348 
1349  switch (ip) {
1352  case diproton: return getRatio(proton)*getRatio(proton);
1353  case unboundPN: return getRatio(proton)*getRatio(neutron);
1354  case dineutron: return getRatio(neutron)*getRatio(neutron);
1355  default: return 0.;
1356  }
1357 
1358  return 0.;
1359 }
1360 
1362  const G4double pn_spec = 1.0; // Scale factor for pn vs. pp/nn
1363  //const G4double pn_spec = 0.5;
1364 
1365  G4double dens = 0.;
1366 
1367  if (ip < 100) dens = getDensity(ip,izone);
1368  else { // For dibaryons, remove extra 1/volume term in density product
1369  switch (ip) {
1370  case diproton:
1371  dens = getDensity(proton,izone) * getDensity(proton,izone);
1372  break;
1373  case unboundPN:
1374  dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
1375  break;
1376  case dineutron:
1377  dens = getDensity(neutron,izone) * getDensity(neutron,izone);
1378  break;
1379  default: dens = 0.;
1380  }
1381  dens *= getVolume(izone);
1382  }
1383 
1384  return getRatio(ip) * dens;
1385 }
1386 
1387 
1390  if (verboseLevel > 1) {
1391  G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
1392  }
1393 
1394  // FIXME: Previous version generated random sin(theta), then used -cos(theta)
1395  // Using generateWithRandomAngles changes result!
1396  // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
1397  G4double costh = std::sqrt(1.0 - inuclRndm());
1399 
1400  // Start particle outside nucleus, unless capture-at-rest
1401  G4int zone = number_of_zones;
1402  if (particle->getKineticEnergy() < small) zone--;
1403 
1404  G4CascadParticle cpart(*particle, pos, zone, large, 0);
1405 
1406  // SPECIAL CASE: Inbound photons are emplanted along through-path
1407  if (forceFirst(cpart)) choosePointAlongTraj(cpart);
1408 
1409  if (verboseLevel > 2) G4cout << cpart << G4endl;
1410 
1411  return cpart;
1412 }
1413 
1416  modelLists& output) {
1417  if (verboseLevel) {
1418  G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
1419  << G4endl;
1420  }
1421 
1422  const G4double max_a_for_cascad = 5.0;
1423  const G4double ekin_cut = 2.0;
1424  const G4double small_ekin = 1.0e-6;
1425  const G4double r_large2for3 = 62.0;
1426  const G4double r0forAeq3 = 3.92;
1427  const G4double s3max = 6.5;
1428  const G4double r_large2for4 = 69.14;
1429  const G4double r0forAeq4 = 4.16;
1430  const G4double s4max = 7.0;
1431  const G4int itry_max = 100;
1432 
1433  // Convenient references to input buffer contents
1434  std::vector<G4CascadParticle>& casparticles = output.first;
1435  std::vector<G4InuclElementaryParticle>& particles = output.second;
1436 
1437  casparticles.clear(); // Reset input buffer to fill this event
1438  particles.clear();
1439 
1440  // first decide whether it will be cascad or compound final nuclei
1441  G4int ab = bullet->getA();
1442  G4int zb = bullet->getZ();
1443  G4int at = target->getA();
1444  G4int zt = target->getZ();
1445 
1446  G4double massb = bullet->getMass(); // For creating LorentzVectors below
1447 
1448  if (ab < max_a_for_cascad) {
1449 
1450  G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
1451  G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
1452  G4double ben = benb < bent ? bent : benb;
1453 
1454  if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
1455  G4int itryg = 0;
1456 
1457  /* Loop checking 08.06.2015 MHK */
1458  while (casparticles.size() == 0 && itryg < itry_max) {
1459  itryg++;
1460  particles.clear();
1461 
1462  // nucleons coordinates and momenta in nuclei rest frame
1463  coordinates.clear();
1464  momentums.clear();
1465 
1466  if (ab < 3) { // deuteron, simplest case
1467  G4double r = 2.214 - 3.4208 * G4Log(1.0 - 0.981 * inuclRndm());
1469  coordinates.push_back(coord1);
1470  coordinates.push_back(-coord1);
1471 
1472  G4double p = 0.0;
1473  G4bool bad = true;
1474  G4int itry = 0;
1475 
1476  while (bad && itry < itry_max) { /* Loop checking 08.06.2015 MHK */
1477  itry++;
1478  p = 456.0 * inuclRndm();
1479 
1480  if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1481  p * r > 312.0) bad = false;
1482  }
1483 
1484  if (itry == itry_max)
1485  if (verboseLevel > 2){
1486  G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;
1487  }
1488 
1489  p = 0.0005 * p;
1490 
1491  if (verboseLevel > 2){
1492  G4cout << " p nuc " << p << G4endl;
1493  }
1494 
1495  G4LorentzVector mom = generateWithRandomAngles(p, massb);
1496 
1497  momentums.push_back(mom);
1498  mom.setVect(-mom.vect());
1499  momentums.push_back(-mom);
1500  } else {
1501  G4ThreeVector coord1;
1502 
1503  G4bool badco = true;
1504 
1505  G4int itry = 0;
1506 
1507  if (ab == 3) {
1508  while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
1509  if (itry > 0) coordinates.clear();
1510  itry++;
1511  G4int i(0);
1512 
1513  for (i = 0; i < 2; i++) {
1514  G4int itry1 = 0;
1515  G4double ss, u, rho;
1516  G4double fmax = G4Exp(-0.5) / std::sqrt(0.5);
1517 
1518  while (itry1 < itry_max) { /* Loop checking 08.06.2015 MHK */
1519  itry1++;
1520  ss = -G4Log(inuclRndm());
1521  u = fmax * inuclRndm();
1522  rho = std::sqrt(ss) * G4Exp(-ss);
1523 
1524  if (rho > u && ss < s3max) {
1525  ss = r0forAeq3 * std::sqrt(ss);
1526  coord1 = generateWithRandomAngles(ss).vect();
1527  coordinates.push_back(coord1);
1528 
1529  if (verboseLevel > 2){
1530  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1531  }
1532  break;
1533  }
1534  }
1535 
1536  if (itry1 == itry_max) { // bad case
1537  coord1.set(10000.,10000.,10000.);
1538  coordinates.push_back(coord1);
1539  break;
1540  }
1541  }
1542 
1543  coord1 = -coordinates[0] - coordinates[1];
1544  if (verboseLevel > 2) {
1545  G4cout << " 3 r " << coord1.mag() << G4endl;
1546  }
1547 
1548  coordinates.push_back(coord1);
1549 
1550  G4bool large_dist = false;
1551 
1552  for (i = 0; i < 2; i++) {
1553  for (G4int j = i+1; j < 3; j++) {
1554  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1555 
1556  if (verboseLevel > 2) {
1557  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1558  }
1559 
1560  if (r2 > r_large2for3) {
1561  large_dist = true;
1562 
1563  break;
1564  }
1565  }
1566 
1567  if (large_dist) break;
1568  }
1569 
1570  if(!large_dist) badco = false;
1571 
1572  }
1573 
1574  } else { // a >= 4
1575  G4double b = 3./(ab - 2.0);
1576  G4double b1 = 1.0 - b / 2.0;
1577  G4double u = b1 + std::sqrt(b1 * b1 + b);
1578  G4double fmax = (1.0 + u/b) * u * G4Exp(-u);
1579 
1580  while (badco && itry < itry_max) {/* Loop checking 08.06.2015 MHK */
1581 
1582  if (itry > 0) coordinates.clear();
1583  itry++;
1584  G4int i(0);
1585 
1586  for (i = 0; i < ab-1; i++) {
1587  G4int itry1 = 0;
1588  G4double ss;
1589 
1590  while (itry1 < itry_max) { /* Loop checking 08.06.2015 MHK */
1591  itry1++;
1592  ss = -G4Log(inuclRndm());
1593  u = fmax * inuclRndm();
1594 
1595  if (std::sqrt(ss) * G4Exp(-ss) * (1.0 + ss/b) > u
1596  && ss < s4max) {
1597  ss = r0forAeq4 * std::sqrt(ss);
1598  coord1 = generateWithRandomAngles(ss).vect();
1599  coordinates.push_back(coord1);
1600 
1601  if (verboseLevel > 2) {
1602  G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1603  }
1604 
1605  break;
1606  }
1607  }
1608 
1609  if (itry1 == itry_max) { // bad case
1610  coord1.set(10000.,10000.,10000.);
1611  coordinates.push_back(coord1);
1612  break;
1613  }
1614  }
1615 
1616  coord1 *= 0.0; // Cheap way to reset
1617  for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
1618 
1619  coordinates.push_back(coord1);
1620 
1621  if (verboseLevel > 2){
1622  G4cout << " last r " << coord1.mag() << G4endl;
1623  }
1624 
1625  G4bool large_dist = false;
1626 
1627  for (i = 0; i < ab-1; i++) {
1628  for (G4int j = i+1; j < ab; j++) {
1629 
1630  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1631 
1632  if (verboseLevel > 2){
1633  G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1634  }
1635 
1636  if (r2 > r_large2for4) {
1637  large_dist = true;
1638 
1639  break;
1640  }
1641  }
1642 
1643  if (large_dist) break;
1644  }
1645 
1646  if (!large_dist) badco = false;
1647  }
1648  }
1649 
1650  if(badco) {
1651  G4cout << " can not generate the nucleons coordinates for a "
1652  << ab << G4endl;
1653 
1654  casparticles.clear(); // Return empty buffer on error
1655  particles.clear();
1656  return;
1657 
1658  } else { // momentums
1659  G4double p, u, x;
1660  G4LorentzVector mom;
1661  //G4bool badp = True;
1662 
1663  for (G4int i = 0; i < ab - 1; i++) {
1664  G4int itry2 = 0;
1665 
1666  while(itry2 < itry_max) { /* Loop checking 08.06.2015 MHK */
1667  itry2++;
1668  u = -G4Log(0.879853 - 0.8798502 * inuclRndm());
1669  x = u * G4Exp(-u);
1670 
1671  if(x > inuclRndm()) {
1672  p = std::sqrt(0.01953 * u);
1673  mom = generateWithRandomAngles(p, massb);
1674  momentums.push_back(mom);
1675 
1676  break;
1677  }
1678  } // while (itry2 < itry_max)
1679 
1680  if(itry2 == itry_max) {
1681  G4cout << " can not generate proper momentum for a "
1682  << ab << G4endl;
1683 
1684  casparticles.clear(); // Return empty buffer on error
1685  particles.clear();
1686  return;
1687  }
1688  } // for (i=0 ...
1689  // last momentum
1690 
1691  mom *= 0.; // Cheap way to reset
1692  mom.setE(bullet->getEnergy()+target->getEnergy());
1693 
1694  for(G4int j=0; j< ab-1; j++) mom -= momentums[j];
1695 
1696  momentums.push_back(mom);
1697  }
1698  }
1699 
1700  // Coordinates and momenta at rest are generated, now back to the lab
1701  G4double rb = 0.0;
1702  G4int i(0);
1703 
1704  for(i = 0; i < G4int(coordinates.size()); i++) {
1705  G4double rp = coordinates[i].mag2();
1706 
1707  if(rp > rb) rb = rp;
1708  }
1709 
1710  // nuclei i.p. as a whole
1711  G4double s1 = std::sqrt(inuclRndm());
1712  G4double phi = randomPHI();
1713  G4double rz = (nuclei_radius + rb) * s1;
1714  G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
1715  -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
1716 
1717  for (i = 0; i < G4int(coordinates.size()); i++) {
1718  coordinates[i] += global_pos;
1719  }
1720 
1721  // all nucleons at rest
1722  raw_particles.clear();
1723 
1724  for (G4int ipa = 0; ipa < ab; ipa++) {
1725  G4int knd = ipa < zb ? 1 : 2;
1726  raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1727  }
1728 
1729  G4InuclElementaryParticle dummy(small_ekin, 1);
1730  G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
1731  toTheBulletRestFrame.toTheTargetRestFrame();
1732 
1734 
1735  for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1736  ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum()));
1737  }
1738 
1739  // fill cascad particles and outgoing particles
1740 
1741  for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1742  G4LorentzVector mom = raw_particles[ip].getMomentum();
1743  G4double pmod = mom.rho();
1744  G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
1745  G4double det = t0 * t0 + nuclei_radius * nuclei_radius
1746  - coordinates[ip].mag2();
1747  G4double tr = -1.0;
1748 
1749  if(det > 0.0) {
1750  G4double t1 = t0 + std::sqrt(det);
1751  G4double t2 = t0 - std::sqrt(det);
1752 
1753  if(std::fabs(t1) <= std::fabs(t2)) {
1754  if(t1 > 0.0) {
1755  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1756  }
1757 
1758  if(tr < 0.0 && t2 > 0.0) {
1759 
1760  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1761  }
1762 
1763  } else {
1764  if(t2 > 0.0) {
1765 
1766  if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1767  }
1768 
1769  if(tr < 0.0 && t1 > 0.0) {
1770  if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1771  }
1772  }
1773 
1774  }
1775 
1776  if(tr >= 0.0) { // cascad particle
1777  coordinates[ip] += mom.vect()*tr / pmod;
1778  casparticles.push_back(G4CascadParticle(raw_particles[ip],
1779  coordinates[ip],
1780  number_of_zones, large, 0));
1781 
1782  } else {
1783  particles.push_back(raw_particles[ip]);
1784  }
1785  }
1786  }
1787 
1788  if(casparticles.size() == 0) {
1789  particles.clear();
1790 
1791  G4cout << " can not generate proper distribution for " << itry_max
1792  << " steps " << G4endl;
1793  }
1794  }
1795  }
1796 
1797  if(verboseLevel > 2){
1798  G4int ip(0);
1799 
1800  G4cout << " cascad particles: " << casparticles.size() << G4endl;
1801  for(ip = 0; ip < G4int(casparticles.size()); ip++)
1802  G4cout << casparticles[ip] << G4endl;
1803 
1804  G4cout << " outgoing particles: " << particles.size() << G4endl;
1805  for(ip = 0; ip < G4int(particles.size()); ip++)
1806  G4cout << particles[ip] << G4endl;
1807  }
1808 
1809  return; // Buffer has been filled
1810 }
1811 
1812 
1813 // Compute mean free path for inclusive interaction of projectile and target
1814 G4double
1817  G4int zone) {
1818  G4int ptype = cparticle.getParticle().type();
1819  G4int ip = target.type();
1820 
1821  // Ensure that zone specified is within nucleus, for array lookups
1822  if (zone<0) zone = cparticle.getCurrentZone();
1823  if (zone>=number_of_zones) zone = number_of_zones-1;
1824 
1825  // Special cases: neutrinos, and muon-on-neutron, have infinite path
1826  if (cparticle.getParticle().isNeutrino()) return 0.;
1827  if (ptype == muonMinus && ip == neutron) return 0.;
1828 
1829  // Configure frame transformation to get kinetic energy for lookups
1830  dummy_convertor.setBullet(cparticle.getParticle());
1831  dummy_convertor.setTarget(&target);
1832  dummy_convertor.toTheCenterOfMass(); // Fill internal kinematics
1834 
1835  // Get cross section for interacting with target (dibaryons are absorptive)
1836  G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
1837  : absorptionCrossSection(ekin, ptype);
1838 
1839  if (verboseLevel > 2) {
1840  G4cout << " ip " << ip << " zone " << zone << " ekin " << ekin
1841  << " dens " << getCurrentDensity(ip, zone)
1842  << " csec " << csec << G4endl;
1843  }
1844 
1845  if (csec <= 0.) return 0.; // No interaction, avoid unnecessary work
1846 
1847  // Get nuclear density and compute mean free path
1848  return csec * getCurrentDensity(ip, zone);
1849 }
1850 
1851 // Throw random distance for interaction of particle using path and MFP
1852 
1853 G4double
1855  G4double path, G4double invmfp) const {
1856  // Delay interactions of newly formed secondaries (minimum int. length)
1857  const G4double young_cut = std::sqrt(10.0) * 0.25;
1858  const G4double huge_num = 50.0; // Argument to exponential
1859 
1860  G4double spath = large; // Buffer for return value
1861 
1862  if (invmfp < small) return spath; // No interaction, avoid unnecessary work
1863 
1864  G4double pw = -path * invmfp; // Ratio of path in zone to MFP
1865  if (pw < -huge_num) pw = -huge_num;
1866  pw = 1.0 - G4Exp(pw);
1867 
1868  if (verboseLevel > 2)
1869  G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
1870 
1871  // Primary particle(s) should always interact at least once
1872  if (forceFirst(cparticle) || (inuclRndm() < pw)) {
1873  spath = -G4Log(1.0 - pw * inuclRndm()) / invmfp;
1874  if (cparticle.young(young_cut, spath)) spath = large;
1875 
1876  if (verboseLevel > 2)
1877  G4cout << " spath " << spath << " path " << path << G4endl;
1878  }
1879 
1880  return spath;
1881 }
1882 
1883 
1884 // Parametrized cross section for pion and photon absorption
1885 
1887  if (!useQuasiDeuteron(type)) {
1888  G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
1889  return 0.;
1890  }
1891 
1892  G4double csec = 0.;
1893 
1894  // Pion absorption is parametrized for low vs. medium energy
1895  // ... use for muon capture as well
1896  if (type == pionPlus || type == pionMinus || type == pionZero ||
1897  type == muonMinus) {
1898  if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
1899  + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
1900  else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);
1901  }
1902 
1903  // Photon cross-section is binned from parametrization by Mi. Kossov
1904  // See below for initialization of gammaQDinterp, gammaQDxsec
1905  if (type == photon) {
1906  csec = gammaQDinterp.interpolate(ke, gammaQDxsec) * gammaQDscale;
1907  }
1908 
1909  if (csec < 0.0) csec = 0.0;
1910 
1911  if (verboseLevel > 2) {
1912  G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;
1913  }
1914 
1915  return crossSectionUnits * csec;
1916 }
1917 
1918 
1920 {
1921  // All scattering cross-sections are available from tables
1922  const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
1923  if (!xsecTable) {
1924  G4cerr << " unknown collison type = " << rtype << G4endl;
1925  return 0.;
1926  }
1927 
1928  return (crossSectionUnits * xsecTable->getCrossSection(ke));
1929 }
Float_t x
Definition: compare.C:6
const XML_Char int len
Definition: expat.h:262
std::vector< G4ThreeVector > collisionPts
G4double nuclei_radius
void set(double x, double y, double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
std::vector< G4double > acsecs
G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
const G4InuclElementaryParticle & getParticle() const
const G4InuclElementaryParticle protonEP
static const G4double pion_vp_small
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static G4double getParticleMass(G4int type)
static const G4double pos
const XML_Char * target
Definition: expat.h:268
TTree * t1
Definition: plottest35.C:26
void fillZoneRadii(G4double nuclearRadius)
const G4double fermiMomentum
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
Int_t ipart
Definition: Style.C:10
std::vector< G4double > vz
G4int getA() const
static const G4double large
std::vector< std::vector< G4double > > fermi_momenta
const G4double radiusForSmall
static const G4double hyperon_vp
void setBullet(const G4InuclParticle *bullet)
G4double getRatio(G4int ip) const
G4bool isProjectile(const G4CascadParticle &cparticle) const
G4CascadeInterpolator< 30 > gammaQDinterp
#define G4endl
Definition: G4ios.hh:61
void updateZone(G4int izone)
Float_t y
Definition: compare.C:6
void fillPotentials(G4int type, G4double tot_vol)
const char * p
Definition: xmltok.h:285
std::vector< std::vector< G4double > > nucleon_densities
Double_t z
void incrementReflectionCounter()
G4bool nucleon(G4int ityp)
G4bool worthToPropagate(const G4CascadParticle &cparticle) const
const G4double R_nucleon
G4bool reflectedNow() const
void printModel() const
double z() const
std::vector< G4double > binding_energies
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
double dot(const Hep3Vector &) const
void setVect(const Hep3Vector &)
G4double v1[6]
void updateParticleMomentum(const G4LorentzVector &mom)
G4bool passTrailing(const G4ThreeVector &hit_position)
static constexpr double second
Definition: G4SIunits.hh:157
G4int getZ() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
std::vector< G4double > zone_volumes
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
void updatePosition(const G4ThreeVector &pos)
#define G4ThreadLocal
Definition: tls.hh:69
std::vector< std::vector< G4double > > zone_potentials
void Register(T *inst)
Definition: G4AutoDelete.hh:65
const G4ThreeVector & getPosition() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double angle(const Hep3Vector &) const
std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
Float_t Z
G4double fillZoneVolumes(G4double nuclearRadius)
G4double totalCrossSection(G4double ke, G4int rtype) const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
static const G4double ab
G4double generateInteractionLength(const G4CascadParticle &cparticle, G4double path, G4double invmfp) const
static const G4CascadeChannel * GetTable(G4int initialState)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
virtual ~G4NucleiModel()
Int_t npart
Definition: Style.C:7
G4double getKineticEnergy() const
G4bool young(G4double young_path_cut, G4double cpath) const
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
static const G4double d2
void propagateAlongThePath(G4double path)
G4double getFermiMomentum(G4int ip, G4int izone) const
static constexpr double deg
Definition: G4SIunits.hh:152
std::vector< G4double > pf
void boundaryTransition(G4CascadParticle &cparticle)
G4double getPotential(G4int ip, G4int izone) const
static const G4double pion_vp
std::vector< G4double > zone_radii
const G4double skinDepth
G4CollisionOutput EPCoutput
#define G4UniformRand()
Definition: Randomize.hh:53
void incrementCurrentPath(G4double npath)
G4double getDensity(G4int ip, G4int izone) const
double A(double temperature)
G4double ur[7]
G4int getZone(G4double r) const
double mag2() const
const G4InuclElementaryParticle neutronEP
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
G4bool movingInsideNuclei() const
static G4bool sortPartners(const partner &p1, const partner &p2)
const G4double crossSectionUnits
G4GLOB_DLL std::ostream G4cerr
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
G4int numberOfOutgoingParticles() const
double rho() const
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector unit() const
TTree * t2
Definition: plottest35.C:36
G4bool passFermi(const std::vector< G4InuclElementaryParticle > &particles, G4int zone)
G4double getVolume(G4int izone) const
static const G4double alfa6[6]
double epsilon(double density, double temperature)
G4int getCurrentZone() const
std::vector< partner > thePartners
double mag() const
void setTarget(const G4InuclParticle *target)
static const G4double alfa3[3]
int G4int
Definition: G4Types.hh:78
const G4double radiusScale2
std::pair< G4InuclElementaryParticle, G4double > partner
void generateInteractionPartners(G4CascadParticle &cparticle)
G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2, G4double nuclearRadius) const
const G4double radiusScale
std::vector< G4InuclElementaryParticle > raw_particles
G4double getCurrentDensity(G4int ip, G4int izone) const
virtual G4double getCrossSection(double ke) const =0
G4InuclNuclei * theNucleus
G4bool forceFirst(const G4CascadParticle &cparticle) const
G4double getCharge() const
G4InuclElementaryParticle generateQuasiDeuteron(G4int type1, G4int type2, G4int zone) const
const G4double radiusUnits
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
G4double getKinEnergyInTheTRS() const
G4GLOB_DLL std::ostream G4cout
double x() const
G4CascadParticle initializeCascad(G4InuclElementaryParticle *particle)
virtual void setVerboseLevel(G4int verbose=0)
static const G4double piTimes4thirds
static const G4double kaon_vp
G4double v[6]
const G4double radScaleAlpha
G4double getFermiKinetic(G4int ip, G4int izone) const
Hep3Vector vect() const
std::vector< G4double > rod
std::vector< G4InuclElementaryParticle > qdeutrons
static constexpr double pi
Definition: G4SIunits.hh:75
G4double getPathToTheNextZone(G4double rz_in, G4double rz_out)
G4double nuclei_volume
double y() const
G4double zoneIntegralGaussian(G4double ur1, G4double ur2, G4double nuclearRadius) const
G4LorentzConvertor dummy_convertor
void printCollisionOutput(std::ostream &os=G4cout) const
G4LorentzVector getMomentum() const
G4double bindingEnergy(G4int A, G4int Z)
static const G4double small
void choosePointAlongTraj(G4CascadParticle &cparticle)
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
G4int neutronNumberCurrent
G4double getEnergy() const
std::vector< G4ThreeVector > coordinates
std::vector< G4LorentzVector > momentums
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4double gammaQDscale
G4int protonNumberCurrent
void generateModel(G4InuclNuclei *nuclei)
G4int getGeneration() const
void fillBindingEnergies()
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
G4double absorptionCrossSection(G4double e, G4int type) const
G4double getMass() const
void generateParticleFate(G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)