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