G4NucleiModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
00029 // 20100114  M. Kelsey -- Use G4ThreeVector for position
00030 // 20100309  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
00031 //              eliminate some unnecessary std::pow(), move ctor here
00032 // 20100407  M. Kelsey -- Replace std::vector<>::resize(0) with ::clear().
00033 //              Move output vectors from generateParticleFate() and 
00034 //              ::generateInteractionPartners() to be data members; return via
00035 //              const-ref instead of by value.
00036 // 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
00037 // 20100418  M. Kelsey -- Reference output particle lists via const-ref
00038 // 20100421  M. Kelsey -- Replace hardwired p/n masses with G4PartDef's
00039 // 20100517  M. Kelsey -- Use G4CascadeINterpolator for cross-section
00040 //              calculations.  use G4Cascade*Channel for total xsec in pi-N
00041 //              N-N channels.  Move absorptionCrossSection() from SpecialFunc.
00042 // 20100610  M. Kelsey -- Replace another random-angle code block; add some
00043 //              diagnostic output for partner-list production.
00044 // 20100617  M. Kelsey -- Replace preprocessor flag CHC_CHECK with
00045 //              G4CASCADE_DEBUG_CHARGE
00046 // 20100620  M. Kelsey -- Improve error message on empty partners list, add
00047 //              four-momentum checking after EPCollider
00048 // 20100621  M. Kelsey -- In boundaryTransition() account for momentum transfer
00049 //              to secondary by boosting into recoil nucleus "rest" frame.
00050 //              Don't need temporaries to create dummy "partners" for list.
00051 // 20100622  M. Kelsey -- Restore use of "bindingEnergy()" function name, which
00052 //              is now a wrapper for G4NucleiProperties::GetBindingEnergy().
00053 // 20100623  M. Kelsey -- Eliminate some temporaries terminating partner-list,
00054 //              discard recoil boost for now. Initialize all data
00055 //              members in ctors.  Allow generateModel() to be called
00056 //              mutliple times, by clearing vectors each time through;
00057 //              avoid extra work by returning if A and Z are same as
00058 //              before.
00059 // 20100628  M. Kelsey -- Two momentum-recoil bugs; don't subtract energies!
00060 // 20100715  M. Kelsey -- Make G4InuclNuclei in generateModel(), use for
00061 //              balance checking (only verbose>2) in generateParticleFate().
00062 // 20100721  M. Kelsey -- Use new G4CASCADE_CHECK_ECONS for conservation checks
00063 // 20100723  M. Kelsey -- Move G4CollisionOutput buffer to .hh for reuse
00064 // 20100726  M. Kelsey -- Preallocate arrays with number_of_zones dimension.
00065 // 20100902  M. Kelsey -- Remove resize(3) directives from qdeutron/acsecs
00066 // 20100907  M. Kelsey -- Limit interaction targets based on current nucleon
00067 //              counts (e.g., no pp if protonNumberCurrent < 2!)
00068 // 20100914  M. Kelsey -- Migrate to integer A and Z
00069 // 20100928  M. Kelsey -- worthToPropagate() use nuclear potential for hadrons.
00070 // 20101005  M. Kelsey -- Annotate hardwired numbers for upcoming rationalizing,
00071 //              move hardwired numbers out to static data members, factor out
00072 //              all the model construction pieces to separate functions, fix
00073 //              wrong value for 4/3 pi.
00074 // 20101019  M. Kelsey -- CoVerity reports: unitialized constructor, dtor leak
00075 // 20101020  M. Kelsey -- Bug fixes to refactoring changes (5 Oct).  Back out
00076 //              worthToPropagate() changes for better regression testing.
00077 // 20101020  M. Kelsey -- Re-activate worthToPropagate() changes.
00078 // 20101119  M. Kelsey -- Hide "negative path" and "no partners" messages in
00079 //              verbosity.
00080 // 20110218  M. Kelsey -- Add crossSectionUnits and radiusUnits scale factors,
00081 //              use "theoretical" numbers for radii etc., multipled by scale
00082 //              factor; set scale factors using environment variables
00083 // 20110303  M. Kelsey -- Add comments why using fabs() with B.E. differences?
00084 // 20110321  M. Kelsey -- Replace strtof() with strtod() for envvar conversion
00085 // 20110321  M. Kelsey -- Use fm and fm^2 as default units, Per D. Wright
00086 //              (NOTE: Restored from original 20110318 commit)
00087 // 20110324  D. Wright -- Implement trailing effect
00088 // 20110324  M. Kelsey -- Move ::reset() here, as it has more code.
00089 // 20110519  M. Kelsey -- Used "rho" after assignment, instead of recomputing
00090 // 20110525  M. Kelsey -- Revert scale factor changes (undo 20110321 changes)
00091 // 20110617  M. Kelsey -- Apply scale factor to trailing-effect radius, make
00092 //              latter runtime adjustable (G4NUCMODEL_RAD_TRAILING)
00093 // 20110720  M. Kelsey -- Follow interface change for cross-section tables,
00094 //              eliminating switch blocks.
00095 // 20110806  M. Kelsey -- Reduce memory churn by pre-allocating buffers
00096 // 20110823  M. Kelsey -- Remove local cross-section tables entirely
00097 // 20110825  M. Kelsey -- Add comments regarding Fermi momentum scale, set of
00098 //              "best guess" parameter values
00099 // 20110831  M. Kelsey -- Make "best guess" parameters the defaults
00100 // 20110922  M. Kelsey -- Follow migrations G4InuclParticle::print(ostream&)
00101 //              and G4CascadParticle::print(ostream&)
00102 // 20111018  M. Kelsey -- Correct kaon potential to be positive, not negative
00103 // 20111107  M. Kelsey -- *** REVERT TO OLD NON-PHYSICAL PARAMETERS FOR 9.5 ***
00104 // 20120306  M. Kelsey -- For incident projectile (CParticle incoming to
00105 //              nucleus from outside) don't use pw cut; force interaction.
00106 // 20120307  M. Kelsey -- Compute zone volumes during setup, not during each
00107 //              interaction.  Include photons in absorption by dibaryons.
00108 //              Discard unused code in totalCrossSection().  Encapsulate
00109 //              interaction path calculations in generateInteractionLength().
00110 // 20120308  M. Kelsey -- Add binned photon absorption cross-section.
00111 // 20120320  M. Kelsey -- Bug fix: fill zone_volumes for single-zone case
00112 // 20120321  M. Kelsey -- Add check in isProjectile() for single-zone case.
00113 // 20120608  M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
00114 // 20120822  M. Kelsey -- Move envvars to G4CascadeParameters.
00115 
00116 #include <numeric>
00117 
00118 #include "G4NucleiModel.hh"
00119 #include "G4PhysicalConstants.hh"
00120 #include "G4SystemOfUnits.hh"
00121 #include "G4CascadeChannel.hh"
00122 #include "G4CascadeChannelTables.hh"
00123 #include "G4CascadeCheckBalance.hh"
00124 #include "G4CascadeInterpolator.hh"
00125 #include "G4CascadeParameters.hh"
00126 #include "G4CollisionOutput.hh"
00127 #include "G4ElementaryParticleCollider.hh"
00128 #include "G4HadTmpUtil.hh"
00129 #include "G4InuclNuclei.hh"
00130 #include "G4InuclParticleNames.hh"
00131 #include "G4InuclSpecialFunctions.hh"
00132 #include "G4LorentzConvertor.hh"
00133 #include "G4Neutron.hh"
00134 #include "G4ParticleLargerBeta.hh"
00135 #include "G4ParticleDefinition.hh"
00136 #include "G4Proton.hh"
00137 
00138 using namespace G4InuclParticleNames;
00139 using namespace G4InuclSpecialFunctions;
00140 
00141 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
00142 
00143 // For the best approximation to a physical-units model, set the following:
00144 //      setenv G4NUCMODEL_XSEC_SCALE   0.1
00145 //      setenv G4NUCMODEL_RAD_SCALE    1.0
00146 //      setenv G4NUCMODEL_RAD_2PAR     1
00147 //      setenv G4NUCMODEL_RAD_SMALL    1.992
00148 //      setenv G4NUCMODEL_RAD_ALPHA    0.84
00149 //      setenv G4NUCMODEL_FERMI_SCALE  0.685
00150 //      setenv G4NUCMODEL_RAD_TRAILING 0.70
00151 
00152 // Scaling factors for radii and cross-sections, currently different!
00153 const G4double G4NucleiModel::crossSectionUnits = G4CascadeParameters::xsecScale();
00154 const G4double G4NucleiModel::radiusUnits =  G4CascadeParameters::radiusScale();
00155 const G4double G4NucleiModel::skinDepth = 0.611207*radiusUnits;
00156 
00157 // One- vs. two-parameter nuclear radius based on envvar
00158 // ==> radius = radiusScale*cbrt(A) + radiusScale2/cbrt(A)
00159 
00160 const G4double G4NucleiModel::radiusScale  = 
00161   (G4CascadeParameters::useTwoParam() ? 1.16 : 1.2) * radiusUnits;
00162 const G4double G4NucleiModel::radiusScale2 =
00163   (G4CascadeParameters::useTwoParam() ? -1.3456 : 0.) * radiusUnits;
00164 
00165 // NOTE:  Old code used R_small = 8.0 (~2.83*units), and R_alpha = 0.7*R_small
00166 // Published data suggests R_small ~ 1.992 fm, R_alpha = 0.84*R_small
00167 
00168 const G4double G4NucleiModel::radiusForSmall = G4CascadeParameters::radiusSmall();
00169 
00170 const G4double G4NucleiModel::radScaleAlpha = G4CascadeParameters::radiusAlpha();
00171 
00172 // Scale factor relating Fermi momentum to density of states, units GeV.fm
00173 // NOTE:  Old code has 0.685*units GeV.fm, literature suggests 0.470 GeV.fm,
00174 //        but this gives too small momentum; old value gives <P_F> ~ 270 MeV
00175 const G4double G4NucleiModel::fermiMomentum = G4CascadeParameters::fermiScale();
00176 
00177 // Effective radius (0.87 to 1.2 fm) of nucleon, for trailing effect
00178 const G4double G4NucleiModel::R_nucleon = G4CascadeParameters::radiusTrailing();
00179 
00180 // Zone boundaries as fraction of nuclear radius (from outside in)
00181 const G4double G4NucleiModel::alfa3[3] = { 0.7, 0.3, 0.01 };
00182 const G4double G4NucleiModel::alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
00183 
00184 // Flat nuclear potentials for mesons and hyperons (GeV)
00185 const G4double G4NucleiModel::pion_vp       = 0.007;
00186 const G4double G4NucleiModel::pion_vp_small = 0.007;
00187 const G4double G4NucleiModel::kaon_vp       = 0.015;
00188 const G4double G4NucleiModel::hyperon_vp    = 0.030;
00189 
00190 const G4double G4NucleiModel::piTimes4thirds = pi*4./3.;
00191 
00192 
00193 // Constructors
00194 G4NucleiModel::G4NucleiModel()
00195   : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00196     A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00197     neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00198     current_nucl2(0) {}
00199 
00200 G4NucleiModel::G4NucleiModel(G4int a, G4int z)
00201   : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00202     A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00203     neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00204     current_nucl2(0) {
00205   generateModel(a,z);
00206 }
00207 
00208 G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei)
00209   : verboseLevel(0), nuclei_radius(0.), nuclei_volume(0.), number_of_zones(0),
00210     A(0), Z(0), theNucleus(0), neutronNumber(0), protonNumber(0),
00211     neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
00212     current_nucl2(0) {
00213   generateModel(nuclei);
00214 }
00215 
00216 G4NucleiModel::~G4NucleiModel() {
00217   delete theNucleus;
00218   theNucleus = 0;
00219 }
00220 
00221 
00222 // Initialize model state for new cascade
00223 
00224 void G4NucleiModel::reset(G4int nHitNeutrons, G4int nHitProtons,
00225                           const std::vector<G4ThreeVector>* hitPoints) {
00226   neutronNumberCurrent = neutronNumber - nHitNeutrons;
00227   protonNumberCurrent  = protonNumber - nHitProtons;
00228   
00229   // zero or copy collision point array for trailing effect
00230   if (!hitPoints || !hitPoints->empty()) collisionPts.clear();
00231   else collisionPts = *hitPoints;
00232 }
00233 
00234 
00235 // Generate nuclear model parameters for given nucleus
00236 
00237 void G4NucleiModel::generateModel(G4InuclNuclei* nuclei) {
00238   generateModel(nuclei->getA(), nuclei->getZ());
00239 }
00240 
00241 void G4NucleiModel::generateModel(G4int a, G4int z) {
00242   if (verboseLevel) {
00243     G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
00244            << G4endl;
00245   }
00246 
00247   // If model already built, just return; otherwise intialize everything
00248   if (a == A && z == Z) {
00249     if (verboseLevel > 1) G4cout << " model already generated" << z << G4endl;
00250     reset();
00251     return;
00252   }
00253 
00254   A = a;
00255   Z = z;
00256   delete theNucleus;
00257   theNucleus = new G4InuclNuclei(A,Z);          // For conservation checking
00258 
00259   neutronNumber = A - Z;
00260   protonNumber = Z;
00261   reset();
00262 
00263   if (verboseLevel > 3) {
00264     G4cout << "  crossSectionUnits = " << crossSectionUnits << G4endl
00265            << "  radiusUnits = " << radiusUnits << G4endl
00266            << "  skinDepth = " << skinDepth << G4endl
00267            << "  radiusScale = " << radiusScale << G4endl
00268            << "  radiusScale2 = " << radiusScale2 << G4endl
00269            << "  radiusForSmall = " << radiusForSmall << G4endl
00270            << "  radScaleAlpha  = " << radScaleAlpha << G4endl
00271            << "  fermiMomentum = " << fermiMomentum << G4endl
00272            << "  piTimes4thirds = " << piTimes4thirds << G4endl;
00273   }
00274 
00275   G4double nuclearRadius;               // Nuclear radius computed from A
00276   if (A>4) nuclearRadius = radiusScale*G4cbrt(A) + radiusScale2/G4cbrt(A);
00277   else     nuclearRadius = radiusForSmall * (A==4 ? radScaleAlpha : 1.);
00278 
00279   // This will be used to pre-allocate lots of arrays below
00280   number_of_zones = (A < 5) ? 1 : (A < 100) ? 3 : 6;
00281 
00282   // Clear all parameters arrays for reloading
00283   binding_energies.clear();
00284   nucleon_densities.clear();
00285   zone_potentials.clear();
00286   fermi_momenta.clear();
00287   zone_radii.clear();
00288   zone_volumes.clear();
00289 
00290   fillBindingEnergies();
00291   fillZoneRadii(nuclearRadius);
00292 
00293   G4double tot_vol = fillZoneVolumes(nuclearRadius);    // Woods-Saxon integral
00294 
00295   fillPotentials(proton, tot_vol);              // Protons
00296   fillPotentials(neutron, tot_vol);             // Neutrons
00297 
00298   // Additional flat zone potentials for other hadrons
00299   const std::vector<G4double> vp(number_of_zones, (A>4)?pion_vp:pion_vp_small);
00300   const std::vector<G4double> kp(number_of_zones, kaon_vp);
00301   const std::vector<G4double> hp(number_of_zones, hyperon_vp);
00302 
00303   zone_potentials.push_back(vp);
00304   zone_potentials.push_back(kp);
00305   zone_potentials.push_back(hp);
00306 
00307   nuclei_radius = zone_radii.back();
00308   nuclei_volume = std::accumulate(zone_volumes.begin(),zone_volumes.end(),0.);
00309 
00310   if (verboseLevel > 3) printModel();
00311 }
00312 
00313 
00314 // Load binding energy array for current nucleus
00315 
00316 void G4NucleiModel::fillBindingEnergies() {
00317   if (verboseLevel > 1)
00318     G4cout << " >>> G4NucleiModel::fillBindingEnergies" << G4endl;
00319 
00320   G4double dm = bindingEnergy(A,Z);
00321 
00322   // Binding energy differences for proton and neutron loss, respectively
00323   // FIXME:  Why is fabs() used here instead of the signed difference?
00324   binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z-1)-dm)/GeV);
00325   binding_energies.push_back(std::fabs(bindingEnergy(A-1,Z)-dm)/GeV);
00326 }
00327 
00328 // Load zone boundary radius array for current nucleus
00329 
00330 void G4NucleiModel::fillZoneRadii(G4double nuclearRadius) {
00331   if (verboseLevel > 1)
00332     G4cout << " >>> G4NucleiModel::fillZoneRadii" << G4endl;
00333 
00334   G4double skinRatio = nuclearRadius/skinDepth;
00335   G4double skinDecay = std::exp(-skinRatio);    
00336 
00337   if (A < 5) {                  // Light ions treated as simple balls
00338     zone_radii.push_back(nuclearRadius);
00339     ur[0] = 0.;
00340     ur[1] = 1.;
00341   } else if (A < 12) {          // Small nuclei have Gaussian potential
00342     G4double rSq = nuclearRadius * nuclearRadius;
00343     G4double gaussRadius = std::sqrt(rSq * (1.0 - 1.0/A) + 6.4);
00344     
00345     ur[0] = 0.0;
00346     for (G4int i = 0; i < number_of_zones; i++) {
00347       G4double y = std::sqrt(-std::log(alfa3[i]));
00348       zone_radii.push_back(gaussRadius * y);
00349       ur[i+1] = y;
00350     }
00351   } else if (A < 100) {         // Intermediate nuclei have three-zone W-S
00352     ur[0] = -skinRatio;
00353     for (G4int i = 0; i < number_of_zones; i++) {
00354       G4double y = std::log((1.0 + skinDecay)/alfa3[i] - 1.0);
00355       zone_radii.push_back(nuclearRadius + skinDepth * y);
00356       ur[i+1] = y;
00357     }
00358   } else {                      // Heavy nuclei have six-zone W-S
00359     ur[0] = -skinRatio;
00360     for (G4int i = 0; i < number_of_zones; i++) {
00361       G4double y = std::log((1.0 + skinDecay)/alfa6[i] - 1.0);
00362       zone_radii.push_back(nuclearRadius + skinDepth * y);
00363       ur[i+1] = y;
00364     }
00365   }
00366 }
00367 
00368 // Compute zone-by-zone density integrals
00369 
00370 G4double G4NucleiModel::fillZoneVolumes(G4double nuclearRadius) {
00371   if (verboseLevel > 1)
00372     G4cout << " >>> G4NucleiModel::fillZoneVolumes" << G4endl;
00373 
00374   G4double tot_vol = 0.;        // Return value omitting 4pi/3 for sphere!
00375 
00376   if (A < 5) {                  // Light ions are treated as simple balls
00377     v[0] = v1[0] = 1.;
00378     tot_vol = zone_radii[0]*zone_radii[0]*zone_radii[0];
00379     zone_volumes.push_back(tot_vol*piTimes4thirds);     // True volume of sphere
00380 
00381     return tot_vol;
00382   }
00383 
00384   PotentialType usePotential = (A < 12) ? Gaussian : WoodsSaxon;
00385 
00386   for (G4int i = 0; i < number_of_zones; i++) {
00387     if (usePotential == WoodsSaxon) {
00388       v[i] = zoneIntegralWoodsSaxon(ur[i], ur[i+1], nuclearRadius);
00389     } else {
00390       v[i] = zoneIntegralGaussian(ur[i], ur[i+1], nuclearRadius);
00391     }
00392     
00393     tot_vol += v[i];
00394     
00395     v1[i] = zone_radii[i]*zone_radii[i]*zone_radii[i];
00396     if (i > 0) v1[i] -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
00397 
00398     zone_volumes.push_back(v1[i]*piTimes4thirds);       // True volume of shell
00399   }
00400 
00401   return tot_vol;       // Sum of zone integrals, not geometric volume
00402 }
00403 
00404 // Load potentials for nucleons (using G4InuclParticle codes)
00405 void G4NucleiModel::fillPotentials(G4int type, G4double tot_vol) {
00406   if (verboseLevel > 1) 
00407     G4cout << " >>> G4NucleiModel::fillZoneVolumes(" << type << ")" << G4endl;
00408 
00409   if (type != proton && type != neutron) return;
00410 
00411   const G4double mass = G4InuclElementaryParticle::getParticleMass(type);
00412 
00413   // FIXME:  This is the fabs() binding energy difference, not signed
00414   const G4double dm = binding_energies[type-1];
00415 
00416   rod.clear(); rod.reserve(number_of_zones);
00417   pf.clear();  pf.reserve(number_of_zones);
00418   vz.clear();  vz.reserve(number_of_zones);
00419 
00420   G4int nNucleons = (type==proton) ? protonNumber : neutronNumber;
00421   G4double dd0 = nNucleons / tot_vol / piTimes4thirds;
00422 
00423   for (G4int i = 0; i < number_of_zones; i++) {
00424     G4double rd = dd0 * v[i] / v1[i];
00425     rod.push_back(rd);
00426     G4double pff = fermiMomentum * G4cbrt(rd);
00427     pf.push_back(pff);
00428     vz.push_back(0.5 * pff * pff / mass + dm);
00429   }
00430   
00431   nucleon_densities.push_back(rod);
00432   fermi_momenta.push_back(pf);
00433   zone_potentials.push_back(vz);
00434 }
00435 
00436 // Zone integral of Woods-Saxon density function
00437 G4double G4NucleiModel::zoneIntegralWoodsSaxon(G4double r1, G4double r2, 
00438                                                G4double nuclearRadius) const {
00439   if (verboseLevel > 1) {
00440     G4cout << " >>> G4NucleiModel::zoneIntegralWoodsSaxon" << G4endl;
00441   }
00442 
00443   const G4double epsilon = 1.0e-3;
00444   const G4int itry_max = 1000;
00445 
00446   G4double skinRatio = nuclearRadius / skinDepth;
00447 
00448   G4double d2 = 2.0 * skinRatio;
00449   G4double dr = r2 - r1;
00450   G4double fr1 = r1 * (r1 + d2) / (1.0 + std::exp(r1));
00451   G4double fr2 = r2 * (r2 + d2) / (1.0 + std::exp(r2));
00452   G4double fi = (fr1 + fr2) / 2.;
00453   G4double fun1 = fi * dr;
00454   G4double fun;
00455   G4double jc = 1;
00456   G4double dr1 = dr;
00457   G4int itry = 0;
00458 
00459   while (itry < itry_max) {
00460     dr /= 2.;
00461     itry++;
00462 
00463     G4double r = r1 - dr;
00464     fi = 0.0;
00465     G4int jc1 = G4int(std::pow(2.0, jc - 1) + 0.1);
00466 
00467     for (G4int i = 0; i < jc1; i++) { 
00468       r += dr1; 
00469       fi += r * (r + d2) / (1.0 + std::exp(r));
00470     }
00471 
00472     fun = 0.5 * fun1 + fi * dr;
00473 
00474     if (std::fabs((fun - fun1) / fun) <= epsilon) break;
00475 
00476     jc++;
00477     dr1 = dr;
00478     fun1 = fun;
00479   }     // while (itry < itry_max)
00480 
00481   if (verboseLevel > 2 && itry == itry_max)
00482     G4cout << " zoneIntegralWoodsSaxon-> n iter " << itry_max << G4endl;
00483 
00484   G4double skinDepth3 = skinDepth*skinDepth*skinDepth;
00485 
00486   return skinDepth3 * (fun + skinRatio*skinRatio*std::log((1.0 + std::exp(-r1)) / (1.0 + std::exp(-r2))));
00487 }
00488 
00489 
00490 // Zone integral of Gaussian density function
00491 G4double G4NucleiModel::zoneIntegralGaussian(G4double r1, G4double r2, 
00492                                              G4double nucRad) const {
00493   if (verboseLevel > 1) {
00494     G4cout << " >>> G4NucleiModel::zoneIntegralGaussian" << G4endl;
00495   }
00496 
00497   G4double gaussRadius = std::sqrt(nucRad*nucRad * (1.0 - 1.0/A) + 6.4);
00498 
00499   const G4double epsilon = 1.0e-3;
00500   const G4int itry_max = 1000;
00501 
00502   G4double dr = r2 - r1;
00503   G4double fr1 = r1 * r1 * std::exp(-r1 * r1);
00504   G4double fr2 = r2 * r2 * std::exp(-r2 * r2);
00505   G4double fi = (fr1 + fr2) / 2.;
00506   G4double fun1 = fi * dr;
00507   G4double fun;
00508   G4double jc = 1;
00509   G4double dr1 = dr;
00510   G4int itry = 0;
00511 
00512   while (itry < itry_max) {
00513     dr /= 2.;
00514     itry++;
00515     G4double r = r1 - dr;
00516     fi = 0.0;
00517     G4int jc1 = int(std::pow(2.0, jc - 1) + 0.1);
00518 
00519     for (G4int i = 0; i < jc1; i++) { 
00520       r += dr1; 
00521       fi += r * r * std::exp(-r * r);
00522     }
00523 
00524     fun = 0.5 * fun1 + fi * dr;  
00525 
00526     if (std::fabs((fun - fun1) / fun) <= epsilon) break;
00527 
00528     jc++;
00529     dr1 = dr;
00530     fun1 = fun;
00531   }     // while (itry < itry_max)
00532 
00533   if (verboseLevel > 2 && itry == itry_max)
00534     G4cerr << " zoneIntegralGaussian-> n iter " << itry_max << G4endl;
00535 
00536   return gaussRadius*gaussRadius*gaussRadius * fun;
00537 }
00538 
00539 
00540 void G4NucleiModel::printModel() const {
00541   if (verboseLevel > 1) {
00542     G4cout << " >>> G4NucleiModel::printModel" << G4endl;
00543   }
00544 
00545   G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
00546          << " proton binding energy " << binding_energies[0]
00547          << " neutron binding energy " << binding_energies[1] << G4endl
00548          << " Nuclei radius " << nuclei_radius << " volume " << nuclei_volume
00549          << " number of zones " << number_of_zones << G4endl;
00550 
00551   for (G4int i = 0; i < number_of_zones; i++)
00552     G4cout << " zone " << i+1 << " radius " << zone_radii[i] 
00553            << " volume " << zone_volumes[i] << G4endl
00554            << " protons: density " << getDensity(1,i) << " PF " << 
00555       getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
00556            << " neutrons: density " << getDensity(2,i) << " PF " << 
00557       getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
00558            << " pions: VP " << getPotential(3,i) << G4endl;
00559 }
00560 
00561 
00562 G4double G4NucleiModel::getFermiKinetic(G4int ip, G4int izone) const {
00563   G4double ekin = 0.0;
00564   
00565   if (ip < 3 && izone < number_of_zones) {      // ip for proton/neutron only
00566     G4double pfermi = fermi_momenta[ip - 1][izone]; 
00567     G4double mass = G4InuclElementaryParticle::getParticleMass(ip);
00568     
00569     ekin = std::sqrt(pfermi*pfermi + mass*mass) - mass;
00570   }  
00571   return ekin;
00572 }
00573 
00574 
00575 G4LorentzVector 
00576 G4NucleiModel::generateNucleonMomentum(G4int type, G4int zone) const {
00577   G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
00578   G4double mass = G4InuclElementaryParticle::getParticleMass(type);
00579 
00580   return generateWithRandomAngles(pmod, mass);
00581 }
00582 
00583 
00584 G4InuclElementaryParticle 
00585 G4NucleiModel::generateNucleon(G4int type, G4int zone) const {
00586   if (verboseLevel > 1) {
00587     G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
00588   }
00589 
00590   G4LorentzVector mom = generateNucleonMomentum(type, zone);
00591   return G4InuclElementaryParticle(mom, type);
00592 }
00593 
00594 
00595 G4InuclElementaryParticle
00596 G4NucleiModel::generateQuasiDeutron(G4int type1, G4int type2,
00597                                     G4int zone) const {
00598   if (verboseLevel > 1) {
00599     G4cout << " >>> G4NucleiModel::generateQuasiDeutron" << G4endl;
00600   }
00601 
00602   // Quasideuteron consists of an unbound but associated nucleon pair
00603   
00604   // FIXME:  Why generate two separate nucleon momenta (randomly!) and
00605   //         add them, instead of just throwing a net momentum for the
00606   //         dinulceon state?  And why do I have to capture the two
00607   //         return values into local variables?
00608   G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
00609   G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
00610   G4LorentzVector dmom = mom1+mom2;
00611 
00612   G4int dtype = 0;
00613        if (type1*type2 == pro*pro) dtype = 111;
00614   else if (type1*type2 == pro*neu) dtype = 112;
00615   else if (type1*type2 == neu*neu) dtype = 122;
00616 
00617   return G4InuclElementaryParticle(dmom, dtype);
00618 }
00619 
00620 
00621 void
00622 G4NucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) {
00623   if (verboseLevel > 1) {
00624     G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
00625   }
00626 
00627   const G4double small = 1.0e-10;       // Cutoff for identifying 0.
00628 
00629   thePartners.clear();          // Reset buffer for next cycle
00630 
00631   G4int ptype = cparticle.getParticle().type();
00632   G4int zone = cparticle.getCurrentZone();
00633 
00634   G4double r_in;                // Destination radius if inbound
00635   G4double r_out;               // Destination radius if outbound
00636 
00637   if (zone == number_of_zones) {
00638     r_in = nuclei_radius;
00639     r_out = 0.0;
00640   } else if (zone == 0) {                       // particle is outside core
00641     r_in = 0.0;
00642     r_out = zone_radii[0];
00643   } else {
00644     r_in = zone_radii[zone - 1];
00645     r_out = zone_radii[zone];
00646   }  
00647 
00648   G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
00649 
00650   if (verboseLevel > 2) {
00651     if (isProjectile(cparticle)) G4cout << " incident particle" << G4endl;
00652     G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path
00653            << G4endl;
00654   }
00655 
00656   if (path < -small) {                  // something wrong
00657     if (verboseLevel)
00658       G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
00659     return;
00660   }
00661 
00662   if (std::fabs(path) < small) {        // just on the boundary
00663     if (verboseLevel > 3)
00664       G4cout << " generateInteractionPartners-> zero path" << G4endl;
00665 
00666     thePartners.push_back(partner());   // Dummy list terminator with zero path
00667     return;
00668   }
00669 
00670   // Initialize frame-transformations for current particle
00671   dummy_convertor.setBullet(cparticle.getParticle());
00672 
00673   G4double invmfp = 0.;                 // Buffers for interaction probability
00674   G4double spath  = 0.;
00675   for (G4int ip = 1; ip < 3; ip++) {
00676     // Only process nucleons which remain active in target
00677     if (ip==proton  && protonNumberCurrent < 1) continue;
00678     if (ip==neutron && neutronNumberCurrent < 1) continue;
00679 
00680     // All nucleons are assumed to be at rest when colliding
00681     G4InuclElementaryParticle particle = generateNucleon(ip, zone);
00682     invmfp = inverseMeanFreePath(cparticle, particle);
00683     spath  = generateInteractionLength(cparticle, path, invmfp);
00684 
00685     if (spath < path) {
00686       if (verboseLevel > 3) {
00687         G4cout << " adding partner[" << thePartners.size() << "]: "
00688                << particle << G4endl;
00689       }
00690       thePartners.push_back(partner(particle, spath));
00691     }
00692   }     // for (G4int ip...
00693   
00694   if (verboseLevel > 2) {
00695     G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
00696   }
00697 
00698   // Absorption possible for pions or photons interacting with dibaryons
00699   if (cparticle.getParticle().pion() || cparticle.getParticle().isPhoton()) {
00700     if (verboseLevel > 2) {
00701       G4cout << " trying quasi-deuterons with bullet: "
00702              << cparticle.getParticle() << G4endl;
00703     }
00704   
00705     // Initialize buffers for quasi-deuteron results
00706     qdeutrons.clear();
00707     acsecs.clear();
00708   
00709     G4double tot_invmfp = 0.0;          // Total inv. mean-free-path for all QDs
00710 
00711     // Proton-proton state interacts with pi- or neutrals
00712     if (protonNumberCurrent >= 2 && ptype != pip) {
00713       G4InuclElementaryParticle ppd = generateQuasiDeutron(pro, pro, zone);
00714       if (verboseLevel > 2)
00715         G4cout << " ptype=" << ptype << " using pp target\n" << ppd << G4endl;
00716       
00717       invmfp = inverseMeanFreePath(cparticle, ppd);      
00718       tot_invmfp += invmfp;
00719       acsecs.push_back(invmfp);
00720       qdeutrons.push_back(ppd);
00721     }
00722     
00723     // Proton-neutron state interacts with any pion type or photon
00724     if (protonNumberCurrent >= 1 && neutronNumberCurrent >= 1) {
00725       G4InuclElementaryParticle npd = generateQuasiDeutron(pro, neu, zone);
00726       if (verboseLevel > 2) 
00727         G4cout << " ptype=" << ptype << " using np target\n" << npd << G4endl;
00728       
00729       invmfp = inverseMeanFreePath(cparticle, npd);
00730       tot_invmfp += invmfp;
00731       acsecs.push_back(invmfp);
00732       qdeutrons.push_back(npd);
00733     }
00734     
00735     // Neutron-neutron state interacts with pi+ or neutrals
00736     if (neutronNumberCurrent >= 2 && ptype != pim) {
00737       G4InuclElementaryParticle nnd = generateQuasiDeutron(neu, neu, zone);
00738       if (verboseLevel > 2)
00739         G4cout << " ptype=" << ptype << " using nn target\n" << nnd << G4endl;
00740       
00741       invmfp = inverseMeanFreePath(cparticle, nnd);
00742       tot_invmfp += invmfp;
00743       acsecs.push_back(invmfp);
00744       qdeutrons.push_back(nnd);
00745     } 
00746     
00747     // Select quasideuteron interaction from non-zero cross-section choices
00748     if (verboseLevel > 2) {
00749       for (size_t i=0; i<qdeutrons.size(); i++) {
00750         G4cout << " acsecs[" << qdeutrons[i].getDefinition()->GetParticleName()
00751                << "] " << acsecs[i];
00752       }
00753       G4cout << G4endl;
00754     }
00755   
00756     if (tot_invmfp > small) {           // Must have absorption cross-section
00757       G4double apath = generateInteractionLength(cparticle, path, tot_invmfp);
00758       
00759       if (apath < path) {               // choose the qdeutron
00760         G4double sl = inuclRndm() * tot_invmfp;
00761         G4double as = 0.0;
00762         
00763         for (size_t i = 0; i < qdeutrons.size(); i++) {
00764           as += acsecs[i];
00765           if (sl < as) { 
00766             if (verboseLevel > 2) G4cout << " deut type " << i << G4endl; 
00767             thePartners.push_back(partner(qdeutrons[i], apath));
00768             break;
00769           }
00770         }       // for (qdeutrons...
00771       }         // if (apath < path)
00772     }           // if (tot_invmfp > small)
00773   }             // pion or photon
00774   
00775   if (verboseLevel > 2) {
00776     G4cout << " after deuterons " << thePartners.size() << " partners"
00777            << G4endl;
00778   }
00779   
00780   if (thePartners.size() > 1) {         // Sort list by path length
00781     std::sort(thePartners.begin(), thePartners.end(), sortPartners);
00782   }
00783   
00784   G4InuclElementaryParticle particle;           // Total path at end of list
00785   thePartners.push_back(partner(particle, path));
00786 }
00787 
00788 
00789 void G4NucleiModel::
00790 generateParticleFate(G4CascadParticle& cparticle,
00791                      G4ElementaryParticleCollider* theEPCollider,
00792                      std::vector<G4CascadParticle>& outgoing_cparticles) {
00793   if (verboseLevel > 1)
00794     G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
00795 
00796   if (verboseLevel > 2) G4cout << " cparticle: " << cparticle << G4endl;
00797 
00798   // Create four-vector checking
00799 #ifdef G4CASCADE_CHECK_ECONS
00800   G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel");  // Second arg is in GeV
00801   balance.setVerboseLevel(verboseLevel);
00802 #endif
00803 
00804   outgoing_cparticles.clear();          // Clear return buffer for this event
00805   generateInteractionPartners(cparticle);       // Fills "thePartners" data
00806 
00807   if (thePartners.empty()) { // smth. is wrong -> needs special treatment
00808     if (verboseLevel)
00809       G4cerr << " generateParticleFate-> got empty interaction-partners list "
00810              << G4endl;
00811     return;
00812   }
00813 
00814   G4int npart = thePartners.size();     // Last item is a total-path placeholder
00815 
00816   if (npart == 1) {             // cparticle is on the next zone entry
00817     cparticle.propagateAlongThePath(thePartners[0].second);
00818     cparticle.incrementCurrentPath(thePartners[0].second);
00819     boundaryTransition(cparticle);
00820     outgoing_cparticles.push_back(cparticle);
00821     
00822     if (verboseLevel > 2) G4cout << " next zone \n" << cparticle << G4endl;
00823   } else {                      // there are possible interactions
00824     if (verboseLevel > 1)
00825       G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
00826 
00827     G4ThreeVector old_position = cparticle.getPosition();
00828     G4InuclElementaryParticle& bullet = cparticle.getParticle();
00829     G4bool no_interaction = true;
00830     G4int zone = cparticle.getCurrentZone();
00831     
00832     for (G4int i=0; i<npart-1; i++) {   // Last item is a total-path placeholder
00833       if (i > 0) cparticle.updatePosition(old_position); 
00834       
00835       G4InuclElementaryParticle& target = thePartners[i].first; 
00836 
00837       if (verboseLevel > 3) {
00838         if (target.quasi_deutron()) G4cout << " try absorption: ";
00839         G4cout << " target " << target.type() << " bullet " << bullet.type()
00840                << G4endl;
00841       }
00842 
00843       EPCoutput.reset();
00844       theEPCollider->collide(&bullet, &target, EPCoutput);
00845 
00846       // If collision failed, exit loop over partners
00847       if (EPCoutput.numberOfOutgoingParticles() == 0) break;
00848 
00849       if (verboseLevel > 2) {
00850         EPCoutput.printCollisionOutput();
00851 #ifdef G4CASCADE_CHECK_ECONS
00852         balance.collide(&bullet, &target, EPCoutput);
00853         balance.okay();         // Do checks, but ignore result
00854 #endif
00855       }
00856 
00857       // Get list of outgoing particles for evaluation
00858       std::vector<G4InuclElementaryParticle>& outgoing_particles = 
00859         EPCoutput.getOutgoingParticles();
00860       
00861       if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
00862 
00863       // Trailing effect: reject interaction at previously hit nucleon
00864       cparticle.propagateAlongThePath(thePartners[i].second);
00865       G4ThreeVector new_position = cparticle.getPosition();
00866 
00867       if (!passTrailing(new_position)) continue;
00868       collisionPts.push_back(new_position);
00869 
00870       // Sort particles according to beta (fastest first)
00871       std::sort(outgoing_particles.begin(), outgoing_particles.end(),
00872                 G4ParticleLargerBeta() );
00873 
00874       if (verboseLevel > 2)
00875         G4cout << " adding " << outgoing_particles.size()
00876                << " output particles" << G4endl;
00877 
00878       // NOTE:  Embedded temporary is optimized away (no copying gets done)
00879       for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 
00880         outgoing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
00881                                                        new_position, zone,
00882                                                        0.0, 0));
00883       }
00884       
00885       no_interaction = false;
00886       current_nucl1 = 0;
00887       current_nucl2 = 0;
00888       
00889 #ifdef G4CASCADE_DEBUG_CHARGE
00890       {
00891         G4double out_charge = 0.0;
00892         
00893         for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) 
00894           out_charge += outgoing_particles[ip].getCharge();
00895         
00896         G4cout << " multiplicity " << outgoing_particles.size()
00897                << " bul type " << bullet.type()
00898                << " targ type " << target.type()
00899                << "\n initial charge "
00900                << bullet.getCharge() + target.getCharge() 
00901                << " out charge " << out_charge << G4endl;  
00902       }
00903 #endif
00904       
00905       if (verboseLevel > 2)
00906         G4cout << " partner type " << target.type() << G4endl;
00907       
00908       if (target.nucleon()) {
00909         current_nucl1 = target.type();
00910       } else {
00911         if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
00912         current_nucl1 = (target.type() - 100) / 10;
00913         current_nucl2 = target.type() - 100 - 10 * current_nucl1;
00914       }   
00915       
00916       if (current_nucl1 == 1) {
00917         if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
00918         protonNumberCurrent--;
00919       } else {
00920         if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
00921         neutronNumberCurrent--;
00922       } 
00923       
00924       if (current_nucl2 == 1) {
00925         if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
00926         protonNumberCurrent--;
00927       } else if (current_nucl2 == 2) {
00928         if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
00929         neutronNumberCurrent--;
00930       }
00931       
00932       break;
00933     }           // loop over partners
00934     
00935     if (no_interaction) {               // still no interactions
00936       if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
00937 
00938       // For conservation checking (below), get particle before updating
00939       static G4InuclElementaryParticle prescatCP;       // Avoid memory churn
00940       prescatCP = cparticle.getParticle();
00941 
00942       // Last "partner" is just a total-path placeholder
00943       cparticle.updatePosition(old_position); 
00944       cparticle.propagateAlongThePath(thePartners[npart-1].second);
00945       cparticle.incrementCurrentPath(thePartners[npart-1].second);
00946       boundaryTransition(cparticle);
00947       outgoing_cparticles.push_back(cparticle);
00948 
00949       // Check conservation for simple scattering (ignore target nucleus!)
00950 #ifdef G4CASCADE_CHECK_ECONS
00951       if (verboseLevel > 2) {
00952         balance.collide(&prescatCP, 0, outgoing_cparticles);
00953         balance.okay();         // Report violations, but don't act on them
00954       }
00955 #endif
00956     }   // if (no_interaction)
00957   }     // if (npart == 1) [else]
00958 
00959   return;
00960 }
00961 
00962 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles, 
00963                                 G4int zone) {
00964   if (verboseLevel > 1) {
00965     G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
00966   }
00967 
00968   // Only check Fermi momenta for nucleons
00969   for (G4int i = 0; i < G4int(particles.size()); i++) {
00970     if (!particles[i].nucleon()) continue;
00971 
00972     G4int type   = particles[i].type();
00973     G4double mom = particles[i].getMomModule();
00974     G4double pfermi  = fermi_momenta[type-1][zone];
00975 
00976     if (verboseLevel > 2)
00977       G4cout << " type " << type << " p " << mom << " pf " << pfermi << G4endl;
00978     
00979     if (mom < pfermi) {
00980         if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
00981         return false;
00982     }
00983   }
00984   return true; 
00985 }
00986 
00987 
00988 // Test here for trailing effect: loop over all previous collision
00989 // locations and test for d > R_nucleon
00990 
00991 G4bool G4NucleiModel::passTrailing(const G4ThreeVector& hit_position) {
00992   if (verboseLevel > 1)
00993     G4cout << " >>> G4NucleiModel::passTrailing " << hit_position << G4endl;
00994 
00995   G4double dist;
00996   for (G4int i = 0; i < G4int(collisionPts.size() ); i++) {
00997     dist = (collisionPts[i] - hit_position).mag();
00998     if (verboseLevel > 2) G4cout << " dist " << dist << G4endl;
00999     if (dist < R_nucleon) {
01000       if (verboseLevel > 2) G4cout << " rejected by Trailing" << G4endl;
01001       return false;
01002     }
01003   }
01004   return true;          // New point far enough away to be used
01005 }
01006 
01007 
01008 void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
01009   if (verboseLevel > 1) {
01010     G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
01011   }
01012 
01013   G4int zone = cparticle.getCurrentZone();
01014 
01015   if (cparticle.movingInsideNuclei() && zone == 0) {
01016     G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
01017     return;
01018   }
01019 
01020   G4LorentzVector mom = cparticle.getMomentum();
01021   G4ThreeVector pos = cparticle.getPosition();
01022   
01023   G4int type = cparticle.getParticle().type();
01024   
01025   G4double r = pos.mag();
01026   G4double pr = pos.dot(mom.vect()) / r;
01027   
01028   G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
01029   
01030   G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
01031   if (verboseLevel > 3) {
01032     G4cout << "Potentials for type " << type << " = " 
01033            << getPotential(type,zone) << " , "
01034            << getPotential(type,next_zone) << G4endl;
01035   }
01036 
01037   G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
01038   
01039   G4double p1r;
01040   
01041   if (verboseLevel > 3) {
01042     G4cout << " type " << type << " zone " << zone << " next " << next_zone
01043            << " qv " << qv << " dv " << dv << G4endl;
01044   }
01045 
01046   if (qv <= 0.0) {      // reflection
01047     if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
01048     p1r = -pr;
01049     cparticle.incrementReflectionCounter();
01050   } else {              // transition
01051     if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
01052     p1r = std::sqrt(qv);
01053     if(pr < 0.0) p1r = -p1r;
01054     cparticle.updateZone(next_zone);
01055     cparticle.resetReflection();
01056   }
01057   
01058   G4double prr = (p1r - pr) / r;  
01059   
01060   if (verboseLevel > 3) {
01061     G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
01062            << prr*pos.y()  << " pz " << prr*pos.z() << " mag "
01063            << std::fabs(prr*r) << G4endl;
01064   }
01065 
01066   mom.setVect(mom.vect() + pos*prr);
01067   cparticle.updateParticleMomentum(mom);
01068 }
01069 
01070 
01071 G4bool G4NucleiModel::isProjectile(const G4CascadParticle& cparticle) const {
01072   if (verboseLevel > 2) {
01073     G4cout << " isProjectile(cparticle): zone "
01074            << cparticle.getCurrentZone() << " of " << number_of_zones
01075            << " path " << cparticle.getCurrentPath()
01076            << " movingInside " << cparticle.movingInsideNuclei()
01077            << " nReflections " << cparticle.getNumberOfReflections()
01078            << G4endl;
01079   }
01080 
01081   return (cparticle.getCurrentZone() == number_of_zones-1 &&    // At surface
01082           cparticle.getCurrentPath() == 1000. &&                // Input primary
01083           cparticle.getNumberOfReflections() <= 0 &&            // first pass
01084           (cparticle.movingInsideNuclei()||number_of_zones==1)  // inbound
01085           );
01086 }
01087 
01088 G4bool G4NucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const {
01089   if (verboseLevel > 1) {
01090     G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
01091   }
01092 
01093   const G4double ekin_scale = 2.0;
01094 
01095   G4bool worth = true;
01096 
01097   if (cparticle.reflectedNow()) {       // Just reflected -- keep going?
01098     G4int zone = cparticle.getCurrentZone();
01099     G4int ip = cparticle.getParticle().type();
01100 
01101     // NOTE:  Temporarily backing out use of potential for non-nucleons
01102     G4double ekin_cut = (cparticle.getParticle().nucleon()) ?
01103       getFermiKinetic(ip, zone) : 0.; //*** getPotential(ip, zone);
01104 
01105     worth = cparticle.getParticle().getKineticEnergy()/ekin_scale > ekin_cut;
01106 
01107     if (verboseLevel > 3) {
01108       G4cout << " type=" << ip
01109              << " ekin=" << cparticle.getParticle().getKineticEnergy()
01110              << " potential=" << ekin_cut
01111              << " : worth? " << worth << G4endl;
01112     }
01113   }
01114 
01115   return worth;
01116 }
01117 
01118 
01119 G4double G4NucleiModel::getRatio(G4int ip) const {
01120   if (verboseLevel > 3) {
01121     G4cout << " >>> G4NucleiModel::getRatio " << ip << G4endl;
01122   }
01123 
01124   switch (ip) {
01125   case proton:    return G4double(protonNumberCurrent)/G4double(protonNumber);
01126   case neutron:   return G4double(neutronNumberCurrent)/G4double(neutronNumber);
01127   case diproton:  return getRatio(proton)*getRatio(proton);
01128   case unboundPN: return getRatio(proton)*getRatio(neutron);
01129   case dineutron: return getRatio(neutron)*getRatio(neutron);
01130   default:        return 0.;
01131   }
01132 
01133   return 0.;
01134 }
01135 
01136 G4double G4NucleiModel::getCurrentDensity(G4int ip, G4int izone) const {
01137   const G4double pn_spec = 1.0;         // Scale factor for pn vs. pp/nn
01138   //const G4double pn_spec = 0.5;
01139 
01140   G4double dens = 0.;
01141 
01142   if (ip < 100) dens = getDensity(ip,izone);
01143   else {        // For dibaryons, remove extra 1/volume term in density product
01144     switch (ip) {
01145     case diproton:  
01146       dens = getDensity(proton,izone) * getDensity(proton,izone);
01147       break;
01148     case unboundPN: 
01149       dens = getDensity(proton,izone) * getDensity(neutron,izone) * pn_spec;
01150       break;
01151     case dineutron:
01152       dens = getDensity(neutron,izone) * getDensity(neutron,izone);
01153       break;
01154     default: dens = 0.;
01155     }
01156     dens *= getVolume(izone);
01157   }
01158 
01159   return getRatio(ip) * dens;
01160 }
01161 
01162 
01163 G4CascadParticle 
01164 G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
01165   if (verboseLevel > 1) {
01166     G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
01167   }
01168 
01169   const G4double large = 1000.0;
01170 
01171   // FIXME:  Previous version generated random sin(theta), then used -cos(theta)
01172   //         Using generateWithRandomAngles changes result!
01173   // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
01174   G4double costh = std::sqrt(1.0 - inuclRndm());
01175   G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
01176 
01177   G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0);
01178 
01179   if (verboseLevel > 2) G4cout << cpart << G4endl;
01180 
01181   return cpart;
01182 }
01183 
01184 void G4NucleiModel::initializeCascad(G4InuclNuclei* bullet, 
01185                                      G4InuclNuclei* target,
01186                                      modelLists& output) {
01187   if (verboseLevel) {
01188     G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
01189            << G4endl;
01190   }
01191 
01192   const G4double large = 1000.0;
01193   const G4double max_a_for_cascad = 5.0;
01194   const G4double ekin_cut = 2.0;
01195   const G4double small_ekin = 1.0e-6;
01196   const G4double r_large2for3 = 62.0;
01197   const G4double r0forAeq3 = 3.92;
01198   const G4double s3max = 6.5;
01199   const G4double r_large2for4 = 69.14;
01200   const G4double r0forAeq4 = 4.16;
01201   const G4double s4max = 7.0;
01202   const G4int itry_max = 100;
01203 
01204   // Convenient references to input buffer contents
01205   std::vector<G4CascadParticle>& casparticles = output.first;
01206   std::vector<G4InuclElementaryParticle>& particles = output.second;
01207 
01208   casparticles.clear();         // Reset input buffer to fill this event
01209   particles.clear();
01210 
01211   // first decide whether it will be cascad or compound final nuclei
01212   G4int ab = bullet->getA();
01213   G4int zb = bullet->getZ();
01214   G4int at = target->getA();
01215   G4int zt = target->getZ();
01216 
01217   G4double massb = bullet->getMass();   // For creating LorentzVectors below
01218 
01219   if (ab < max_a_for_cascad) {
01220 
01221     G4double benb = bindingEnergy(ab,zb)/GeV / G4double(ab);
01222     G4double bent = bindingEnergy(at,zt)/GeV / G4double(at);
01223     G4double ben = benb < bent ? bent : benb;
01224 
01225     if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
01226       G4int itryg = 0;
01227 
01228       while (casparticles.size() == 0 && itryg < itry_max) {      
01229         itryg++;
01230         particles.clear();
01231       
01232         //    nucleons coordinates and momenta in nuclei rest frame
01233         coordinates.clear();
01234         momentums.clear();
01235      
01236         if (ab < 3) { // deuteron, simplest case
01237           G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
01238           G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
01239           coordinates.push_back(coord1);
01240           coordinates.push_back(-coord1);
01241 
01242           G4double p = 0.0;
01243           G4bool bad = true;
01244           G4int itry = 0;
01245 
01246           while (bad && itry < itry_max) {
01247             itry++;
01248             p = 456.0 * inuclRndm();
01249 
01250             if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
01251                p * r > 312.0) bad = false;
01252           }
01253 
01254           if (itry == itry_max)
01255             if (verboseLevel > 2){ 
01256               G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;   
01257             }
01258 
01259           p = 0.0005 * p;
01260 
01261           if (verboseLevel > 2){ 
01262             G4cout << " p nuc " << p << G4endl;
01263           }
01264 
01265           G4LorentzVector mom = generateWithRandomAngles(p, massb);
01266 
01267           momentums.push_back(mom);
01268           mom.setVect(-mom.vect());
01269           momentums.push_back(-mom);
01270         } else {
01271           G4ThreeVector coord1;
01272 
01273           G4bool badco = true;
01274 
01275           G4int itry = 0;
01276         
01277           if (ab == 3) {
01278             while (badco && itry < itry_max) {
01279               if (itry > 0) coordinates.clear();
01280               itry++;   
01281               G4int i(0);    
01282 
01283               for (i = 0; i < 2; i++) {
01284                 G4int itry1 = 0;
01285                 G4double ss, u, rho; 
01286                 G4double fmax = std::exp(-0.5) / std::sqrt(0.5);
01287 
01288                 while (itry1 < itry_max) {
01289                   itry1++;
01290                   ss = -std::log(inuclRndm());
01291                   u = fmax * inuclRndm();
01292                   rho = std::sqrt(ss) * std::exp(-ss);
01293 
01294                   if (rho > u && ss < s3max) {
01295                     ss = r0forAeq3 * std::sqrt(ss);
01296                     coord1 = generateWithRandomAngles(ss).vect();
01297                     coordinates.push_back(coord1);
01298 
01299                     if (verboseLevel > 2){
01300                       G4cout << " i " << i << " r " << coord1.mag() << G4endl;
01301                     }
01302                     break;
01303                   }
01304                 }
01305 
01306                 if (itry1 == itry_max) { // bad case
01307                   coord1.set(10000.,10000.,10000.);
01308                   coordinates.push_back(coord1);
01309                   break;
01310                 }
01311               }
01312 
01313               coord1 = -coordinates[0] - coordinates[1]; 
01314               if (verboseLevel > 2) {
01315                 G4cout << " 3  r " << coord1.mag() << G4endl;
01316               }
01317 
01318               coordinates.push_back(coord1);        
01319             
01320               G4bool large_dist = false;
01321 
01322               for (i = 0; i < 2; i++) {
01323                 for (G4int j = i+1; j < 3; j++) {
01324                   G4double r2 = (coordinates[i]-coordinates[j]).mag2();
01325 
01326                   if (verboseLevel > 2) {
01327                     G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
01328                   }
01329 
01330                   if (r2 > r_large2for3) {
01331                     large_dist = true;
01332 
01333                     break; 
01334                   }      
01335                 }
01336 
01337                 if (large_dist) break;
01338               } 
01339 
01340               if(!large_dist) badco = false;
01341 
01342             }
01343 
01344           } else { // a >= 4
01345             G4double b = 3./(ab - 2.0);
01346             G4double b1 = 1.0 - b / 2.0;
01347             G4double u = b1 + std::sqrt(b1 * b1 + b);
01348             G4double fmax = (1.0 + u/b) * u * std::exp(-u);
01349           
01350             while (badco && itry < itry_max) {
01351 
01352               if (itry > 0) coordinates.clear();
01353               itry++;
01354               G4int i(0);
01355             
01356               for (i = 0; i < ab-1; i++) {
01357                 G4int itry1 = 0;
01358                 G4double ss; 
01359 
01360                 while (itry1 < itry_max) {
01361                   itry1++;
01362                   ss = -std::log(inuclRndm());
01363                   u = fmax * inuclRndm();
01364 
01365                   if (std::sqrt(ss) * std::exp(-ss) * (1.0 + ss/b) > u
01366                       && ss < s4max) {
01367                     ss = r0forAeq4 * std::sqrt(ss);
01368                     coord1 = generateWithRandomAngles(ss).vect();
01369                     coordinates.push_back(coord1);
01370 
01371                     if (verboseLevel > 2) {
01372                       G4cout << " i " << i << " r " << coord1.mag() << G4endl;
01373                     }
01374 
01375                     break;
01376                   }
01377                 }
01378 
01379                 if (itry1 == itry_max) { // bad case
01380                   coord1.set(10000.,10000.,10000.);
01381                   coordinates.push_back(coord1);
01382                   break;
01383                 }
01384               }
01385 
01386               coord1 *= 0.0;    // Cheap way to reset
01387               for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
01388 
01389               coordinates.push_back(coord1);   
01390 
01391               if (verboseLevel > 2){
01392                 G4cout << " last r " << coord1.mag() << G4endl;
01393               }
01394             
01395               G4bool large_dist = false;
01396 
01397               for (i = 0; i < ab-1; i++) {
01398                 for (G4int j = i+1; j < ab; j++) {
01399              
01400                   G4double r2 = (coordinates[i]-coordinates[j]).mag2();
01401 
01402                   if (verboseLevel > 2){
01403                     G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
01404                   }
01405 
01406                   if (r2 > r_large2for4) {
01407                     large_dist = true;
01408 
01409                     break; 
01410                   }      
01411                 }
01412 
01413                 if (large_dist) break;
01414               } 
01415 
01416               if (!large_dist) badco = false;
01417             }
01418           } 
01419 
01420           if(badco) {
01421             G4cout << " can not generate the nucleons coordinates for a "
01422                    << ab << G4endl;     
01423 
01424             casparticles.clear();       // Return empty buffer on error
01425             particles.clear();
01426             return;
01427 
01428           } else { // momentums
01429             G4double p, u, x;
01430             G4LorentzVector mom;
01431             //G4bool badp = True;
01432 
01433             for (G4int i = 0; i < ab - 1; i++) {
01434               G4int itry2 = 0;
01435 
01436               while(itry2 < itry_max) {
01437                 itry2++;
01438                 u = -std::log(0.879853 - 0.8798502 * inuclRndm());
01439                 x = u * std::exp(-u);
01440 
01441                 if(x > inuclRndm()) {
01442                   p = std::sqrt(0.01953 * u);
01443                   mom = generateWithRandomAngles(p, massb);
01444                   momentums.push_back(mom);
01445 
01446                   break;
01447                 }
01448               } // while (itry2 < itry_max)
01449 
01450               if(itry2 == itry_max) {
01451                 G4cout << " can not generate proper momentum for a "
01452                        << ab << G4endl;
01453 
01454                 casparticles.clear();   // Return empty buffer on error
01455                 particles.clear();
01456                 return;
01457               } 
01458             }   // for (i=0 ...
01459             // last momentum
01460 
01461             mom *= 0.;          // Cheap way to reset
01462             mom.setE(bullet->getEnergy()+target->getEnergy());
01463 
01464             for(G4int j=0; j< ab-1; j++) mom -= momentums[j]; 
01465 
01466             momentums.push_back(mom);
01467           } 
01468         }
01469  
01470         // Coordinates and momenta at rest are generated, now back to the lab
01471         G4double rb = 0.0;
01472         G4int i(0);
01473 
01474         for(i = 0; i < G4int(coordinates.size()); i++) {      
01475           G4double rp = coordinates[i].mag2();
01476 
01477           if(rp > rb) rb = rp;
01478         }
01479 
01480         // nuclei i.p. as a whole
01481         G4double s1 = std::sqrt(inuclRndm()); 
01482         G4double phi = randomPHI();
01483         G4double rz = (nuclei_radius + rb) * s1;
01484         G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
01485                                  -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
01486 
01487         for (i = 0; i < G4int(coordinates.size()); i++) {
01488           coordinates[i] += global_pos;
01489         }  
01490 
01491         // all nucleons at rest
01492         raw_particles.clear();
01493 
01494         for (G4int ipa = 0; ipa < ab; ipa++) {
01495           G4int knd = ipa < zb ? 1 : 2;
01496           raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
01497         } 
01498       
01499         G4InuclElementaryParticle dummy(small_ekin, 1);
01500         G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
01501         toTheBulletRestFrame.toTheTargetRestFrame();
01502 
01503         particleIterator ipart;
01504 
01505         for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
01506           ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 
01507         }
01508 
01509         // fill cascad particles and outgoing particles
01510 
01511         for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
01512           G4LorentzVector mom = raw_particles[ip].getMomentum();
01513           G4double pmod = mom.rho();
01514           G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
01515           G4double det = t0 * t0 + nuclei_radius * nuclei_radius
01516                        - coordinates[ip].mag2();
01517           G4double tr = -1.0;
01518 
01519           if(det > 0.0) {
01520             G4double t1 = t0 + std::sqrt(det);
01521             G4double t2 = t0 - std::sqrt(det);
01522 
01523             if(std::fabs(t1) <= std::fabs(t2)) {         
01524               if(t1 > 0.0) {
01525                 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
01526               }
01527 
01528               if(tr < 0.0 && t2 > 0.0) {
01529 
01530                 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
01531               }
01532 
01533             } else {
01534               if(t2 > 0.0) {
01535 
01536                 if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
01537               }
01538 
01539               if(tr < 0.0 && t1 > 0.0) {
01540                 if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
01541               }
01542             } 
01543 
01544           }
01545 
01546           if(tr >= 0.0) { // cascad particle
01547             coordinates[ip] += mom.vect()*tr / pmod;
01548             casparticles.push_back(G4CascadParticle(raw_particles[ip], 
01549                                                     coordinates[ip], 
01550                                                     number_of_zones, large, 0));
01551 
01552           } else {
01553             particles.push_back(raw_particles[ip]); 
01554           } 
01555         }
01556       }    
01557 
01558       if(casparticles.size() == 0) {
01559         particles.clear();
01560 
01561         G4cout << " can not generate proper distribution for " << itry_max
01562                << " steps " << G4endl;
01563       }    
01564     }
01565   }
01566 
01567   if(verboseLevel > 2){
01568     G4int ip(0);
01569 
01570     G4cout << " cascad particles: " << casparticles.size() << G4endl;
01571     for(ip = 0; ip < G4int(casparticles.size()); ip++)
01572       G4cout << casparticles[ip] << G4endl;
01573 
01574     G4cout << " outgoing particles: " << particles.size() << G4endl;
01575     for(ip = 0; ip < G4int(particles.size()); ip++)
01576       G4cout << particles[ip] << G4endl;
01577   }
01578 
01579   return;       // Buffer has been filled
01580 }
01581 
01582 
01583 // Compute mean free path for inclusive interaction of projectile and target
01584 G4double 
01585 G4NucleiModel::inverseMeanFreePath(const G4CascadParticle& cparticle,
01586                                    const G4InuclElementaryParticle& target) {
01587   G4int ptype = cparticle.getParticle().type();
01588   G4int ip = target.type();
01589   G4int zone = cparticle.getCurrentZone();
01590 
01591   // Configure frame transformation to get kinetic energy for lookups
01592   dummy_convertor.setBullet(cparticle.getParticle());
01593   dummy_convertor.setTarget(&target);
01594   G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
01595 
01596   // Get cross section for interacting with target (dibaryons are absorptive)
01597   G4double csec = (ip < 100) ? totalCrossSection(ekin, ptype*ip)
01598                              : absorptionCrossSection(ekin, ptype);
01599 
01600   if(verboseLevel > 2) {
01601     G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
01602   }
01603 
01604   if (csec <= 0.) return 0.;            // No interaction, avoid unnecessary work
01605 
01606   // Get nuclear density and compute mean free path
01607   return csec * getCurrentDensity(ip, zone);
01608 }
01609 
01610 // Throw random distance for interaction of particle using path and MFP
01611 
01612 G4double 
01613 G4NucleiModel::generateInteractionLength(const G4CascadParticle& cparticle,
01614                                          G4double path, G4double invmfp) const {
01615   // Delay interactions of newly formed secondaries (minimum int. length)
01616   const G4double young_cut = std::sqrt(10.0) * 0.25;
01617   const G4double huge_num = 50.0;       // Argument to exponential
01618   const G4double small = 1.0e-10;
01619 
01620   if (invmfp < small) return path;      // No interaction, avoid unnecessary work
01621 
01622   G4double pw = -path * invmfp;
01623   if (pw < -huge_num) pw = -huge_num;
01624   pw = 1.0 - std::exp(pw);
01625     
01626   if (verboseLevel > 2) 
01627     G4cout << " mfp " << 1./invmfp << " pw " << pw << G4endl;
01628     
01629   G4double spath = path;                // Buffer for return value
01630 
01631   // Primary particle(s) should always interact at least once
01632   if (isProjectile(cparticle) || (inuclRndm() < pw)) {
01633     spath = -std::log(1.0 - pw * inuclRndm()) / invmfp;
01634     if (cparticle.young(young_cut, spath)) spath = path;
01635     
01636     if (verboseLevel > 2) 
01637       G4cout << " spath " << spath << " path " << path << G4endl;
01638   }
01639 
01640   return spath;
01641 }
01642 
01643 
01644 // Parametrized cross section for pion and photon absorption
01645 
01646 G4double G4NucleiModel::absorptionCrossSection(G4double ke, G4int type) const {
01647   if (type != pionPlus && type != pionMinus && type != pionZero &&
01648       type != photon) {
01649     G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
01650     return 0.;
01651   }
01652 
01653   G4double csec = 0.;
01654 
01655   // Pion absorption is parametrized for low vs. medium energy
01656   if (type == pionPlus || type == pionMinus || type == pionZero) {
01657     if (ke < 0.3) csec = (0.1106 / std::sqrt(ke) - 0.8
01658                           + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056) );
01659     else if (ke < 1.0) csec = 3.6735 * (1.0-ke)*(1.0-ke);     
01660   }
01661 
01662   // Photon cross-section is binned from parametrization by Mi. Kossov
01663   if (type == photon) {
01664     static const G4double gammaQDscale = G4CascadeParameters::gammaQDScale();
01665     static const G4double kebins[] =
01666       {  0.0,  0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
01667          0.13, 0.18, 0.24,  0.32,  0.42,  0.56,  0.75,  1.0,   1.3,   1.8,
01668          2.4,  3.2,  4.2,   5.6,   7.5,   10.0,  13.0,  18.0,  24.0, 32.0 };
01669     static const G4double gammaD[] =
01670       { 2.8,    1.3,    0.89,   0.56,   0.38,   0.27,   0.19,   0.14,   0.098,
01671         0.071, 0.054,  0.0003, 0.0007, 0.0027, 0.0014, 0.001,  0.0012, 0.0005,
01672         0.0003, 0.0002,0.0002, 0.0002, 0.0002, 0.0002, 0.0001, 0.0001, 0.0001,
01673         0.0001, 0.0001, 0.0001 };
01674     static G4CascadeInterpolator<30> interp(kebins);
01675     csec = interp.interpolate(ke, gammaD) * gammaQDscale;
01676   }
01677 
01678   if (csec < 0.0) csec = 0.0;
01679 
01680   if (verboseLevel > 2) {
01681     G4cout << " ekin " << ke << " abs. csec " << csec << " mb" << G4endl;   
01682   }
01683 
01684   return crossSectionUnits * csec;
01685 }
01686 
01687 
01688 G4double G4NucleiModel::totalCrossSection(G4double ke, G4int rtype) const
01689 {
01690   // All scattering cross-sections are available from tables
01691   const G4CascadeChannel* xsecTable = G4CascadeChannelTables::GetTable(rtype);
01692   if (!xsecTable) {
01693     G4cerr << " unknown collison type = " << rtype << G4endl;
01694     return 0.;
01695   }
01696 
01697   return (crossSectionUnits * xsecTable->getCrossSection(ke));
01698 }

Generated on Mon May 27 17:49:07 2013 for Geant4 by  doxygen 1.4.7