Geant4  10.02.p01
G4NucleiModel.hh
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.hh 71940 2013-06-28 19:04:44Z mkelsey $
27 // GEANT4 tag: $Name: $
28 //
29 // 20100319 M. Kelsey -- Remove "using" directory and unnecessary #includes,
30 // move ctor to .cc file
31 // 20100407 M. Kelsey -- Create "partners thePartners" data member to act
32 // as buffer between ::generateInteractionPartners() and
33 // ::generateParticleFate(), and make "outgoing_cparticles" a
34 // data member returned from the latter by const-ref. Replace
35 // return-by-value of initializeCascad() with an input buffer.
36 // 20100409 M. Kelsey -- Add function to sort list of partnerts by pathlen,
37 // move non-inlinable code to .cc.
38 // 20100421 M. Kelsey -- Move getFermiKinetic() to .cc, no hardwired masses.
39 // 20100517 M. Kelsey -- Change cross-section tables to static arrays. Move
40 // absorptionCrossSection() from SpecialFunc.
41 // 20100520 M. Kelsey -- Add function to separate momentum from nucleon
42 // 20100617 M. Kelsey -- Add setVerboseLevel() function, add generateModel()
43 // with particle input, and ctor with A/Z input.
44 // 20100715 M. Kelsey -- Add G4InuclNuclei object for use with balance checks
45 // 20100723 M. Kelsey -- Move G4CollisionOutput buffer, along with all
46 // std::vector<> buffers, here for reuse
47 // 20100914 M. Kelsey -- Migrate to integer A and Z
48 // 20101004 M. Kelsey -- Rename and create functions used to generate model
49 // 20101019 M. Kelsey -- CoVerity report: dtor leak; move dtor to .cc file
50 // 20110223 M. Kelsey -- Add static parameters for radius and cross-section
51 // scaling factors.
52 // 20110303 M. Kelsey -- Add accessors for model parameters and units
53 // 20110304 M. Kelsey -- Extend reset() to fill neutron and proton counts
54 // 20110324 D. Wright -- Add list of nucleon interaction points, and nucleon
55 // effective radius, for trailing effect.
56 // 20110324 M. Kelsey -- Extend reset() to pass list of points; move
57 // implementation to .cc file.
58 // 20110405 M. Kelsey -- Add "passTrailing()" to encapsulate trailing effect
59 // 20110808 M. Kelsey -- Pass buffer into generateParticleFate instead of
60 // returning a vector to be copied.
61 // 20110823 M. Kelsey -- Remove local cross-section tables entirely
62 // 20120125 M. Kelsey -- Add special case for photons to have zero potential.
63 // 20120307 M. Kelsey -- Add zone volume array for use with quasideuteron
64 // density, functions to identify projectile, compute interaction
65 // distance.
66 // 20130129 M. Kelsey -- Add static objects for gamma-quasideuteron scattering
67 // 20130221 M. Kelsey -- Add function to emplant particle along trajectory
68 // 20130226 M. Kelsey -- Allow forcing zone selection in MFP calculation.
69 // 20130302 M. Kelsey -- Add forceFirst() wrapper function (allows configuring)
70 // 20130628 M. Kelsey -- Extend useQuasiDeuteron() to check good absorption,
71 // fix spelling of "Deutron" -> "Deuteron"
72 // 20130808 M. Kelsey -- To avoid thread collisions, move static neutronEP
73 // and protonEP objects to const data members.
74 // 20131001 M. Kelsey -- Move QDinterp object to data member, thread local
75 // 20140116 M. Kelsey -- Move statics to const data members to avoid weird
76 // interactions with MT.
77 
78 #ifndef G4NUCLEI_MODEL_HH
79 #define G4NUCLEI_MODEL_HH
80 
81 #include <algorithm>
82 #include <vector>
83 
85 #include "G4CascadParticle.hh"
86 #include "G4CascadeInterpolator.hh"
87 #include "G4CollisionOutput.hh"
88 #include "G4LorentzConvertor.hh"
89 
90 class G4InuclNuclei;
92 
94 public:
95  G4NucleiModel();
98 
99  virtual ~G4NucleiModel();
100 
101  void setVerboseLevel(G4int verbose) { verboseLevel = verbose; }
102 
103  void generateModel(G4InuclNuclei* nuclei);
104  void generateModel(G4int a, G4int z);
105 
106  // Arguments here are used for rescattering (::Propagate)
107  void reset(G4int nHitNeutrons=0, G4int nHitProtons=0,
108  const std::vector<G4ThreeVector>* hitPoints=0);
109 
110  void printModel() const;
111 
112  G4double getDensity(G4int ip, G4int izone) const {
113  return nucleon_densities[ip - 1][izone];
114  }
115 
117  return fermi_momenta[ip - 1][izone];
118  }
119 
120  G4double getFermiKinetic(G4int ip, G4int izone) const;
121 
122  G4double getPotential(G4int ip, G4int izone) const {
123  if (ip == 9 || ip < 0) return 0.0; // Photons and leptons
124  G4int ip0 = ip < 3 ? ip - 1 : 2;
125  if (ip > 10 && ip < 18) ip0 = 3;
126  if (ip > 20) ip0 = 4;
127  return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
128  }
129 
130  // Factor to convert GEANT4 lengths to internal units
132 
133  G4double getRadius() const { return nuclei_radius; }
134  G4double getRadius(G4int izone) const {
135  return ( (izone<0) ? 0.
136  : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
137  }
138  G4double getVolume(G4int izone) const {
139  return ( (izone<0) ? 0.
140  : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
141  }
142 
144  G4int getZone(G4double r) const {
145  for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
146  return number_of_zones;
147  }
148 
151 
152  G4bool empty() const {
153  return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
154  }
155 
156  G4bool stillInside(const G4CascadParticle& cparticle) {
157  return cparticle.getCurrentZone() < number_of_zones;
158  }
159 
160 
162 
163  typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > modelLists;
164 
165  void initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target,
166  modelLists& output);
167 
168 
169  std::pair<G4int, G4int> getTypesOfNucleonsInvolved() const {
170  return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
171  }
172 
173  void generateParticleFate(G4CascadParticle& cparticle,
174  G4ElementaryParticleCollider* theEPCollider,
175  std::vector<G4CascadParticle>& cascade);
176 
177  G4bool forceFirst(const G4CascadParticle& cparticle) const;
178  G4bool isProjectile(const G4CascadParticle& cparticle) const;
179  G4bool worthToPropagate(const G4CascadParticle& cparticle) const;
180 
182 
184 
186  G4double totalCrossSection(G4double ke, G4int rtype) const;
187 
188  // Identify whether given particle can interact with dibaryons
189  static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0);
190 
191 protected:
192  G4bool passFermi(const std::vector<G4InuclElementaryParticle>& particles,
193  G4int zone);
194 
195  G4bool passTrailing(const G4ThreeVector& hit_position);
196 
197  void boundaryTransition(G4CascadParticle& cparticle);
198 
199  void choosePointAlongTraj(G4CascadParticle& cparticle);
200 
202  G4int type2,
203  G4int zone) const;
204 
205  typedef std::pair<G4InuclElementaryParticle, G4double> partner;
206 
207  std::vector<partner> thePartners; // Buffer for output below
209 
210  // Function for std::sort() to use in organizing partners by path length
211  static G4bool sortPartners(const partner& p1, const partner& p2) {
212  return (p2.second > p1.second);
213  }
214 
215  // Functions used to generate model nuclear structure
216  void fillBindingEnergies();
217 
218  void fillZoneRadii(G4double nuclearRadius);
219 
220  G4double fillZoneVolumes(G4double nuclearRadius);
221 
222  void fillPotentials(G4int type, G4double tot_vol);
223 
225  G4double nuclearRadius) const;
226 
228  G4double nuclearRadius) const;
229 
230  G4double getRatio(G4int ip) const; // Fraction of nucleons still available
231 
232  // Scale nuclear density with interactions so far
233  G4double getCurrentDensity(G4int ip, G4int izone) const;
234 
235  // Average path length for any interaction of projectile and target
236  // = 1. / (density * cross-section)
238  const G4InuclElementaryParticle& target,
239  G4int zone = -1); // Override zone value
240  // NOTE: Function is non-const in order to use dummy_converter
241 
242  // Use path-length and MFP (above) to throw random distance to collision
244  G4double path, G4double invmfp) const;
245 
246 private:
248 
249  // Buffers for processing interactions on each cycle
250  G4LorentzConvertor dummy_convertor; // For getting collision frame
251  G4CollisionOutput EPCoutput; // For N-body inelastic collisions
252 
253  std::vector<G4InuclElementaryParticle> qdeutrons; // For h+(NN) trials
254  std::vector<G4double> acsecs;
255 
256  std::vector<G4ThreeVector> coordinates; // for initializeCascad()
257  std::vector<G4LorentzVector> momentums;
258  std::vector<G4InuclElementaryParticle> raw_particles;
259 
260  std::vector<G4ThreeVector> collisionPts;
261 
262  // Temporary buffers for computing nuclear model
263  G4double ur[7]; // Number of skin depths for each zone
264  G4double v[6]; // Density integrals by zone
265  G4double v1[6]; // Pseudo-volume (delta r^3) by zone
266  std::vector<G4double> rod; // Nucleon density
267  std::vector<G4double> pf; // Fermi momentum
268  std::vector<G4double> vz; // Nucleon potential
269 
270  // Nuclear model configuration
271  std::vector<std::vector<G4double> > nucleon_densities;
272  std::vector<std::vector<G4double> > zone_potentials;
273  std::vector<std::vector<G4double> > fermi_momenta;
274  std::vector<G4double> zone_radii;
275  std::vector<G4double> zone_volumes;
276  std::vector<G4double> binding_energies;
280 
284 
287 
290 
293 
294  G4CascadeInterpolator<30> gammaQDinterp; // quasideuteron interpolator
295 
296  // Symbolic names for nuclear potentials
298 
299  // Parameters for nuclear structure
300  const G4double crossSectionUnits; // Scale from internal to natural units
302  const G4double skinDepth; // Fraction of radius for outer skin
303  const G4double radiusScale; // Coefficients for two-parameter fit
304  const G4double radiusScale2; // R = 1.16*cbrt(A) - 1.3456/cbrt(A)
305  const G4double radiusForSmall; // Average radius of light A<5 nuclei
306  const G4double radScaleAlpha; // Scaling factor R_alpha/R_small
309  const G4double gammaQDscale; // Gamma/cluster scattering rescaling
310 
311  // Cutoffs for extreme values
312  static const G4double small;
313  static const G4double large;
314  static const G4double piTimes4thirds; // FIXME: We should not be using this!
315 
316  static const G4double alfa3[3], alfa6[6]; // Zone boundaries in radii
317  static const G4double pion_vp; // Flat potentials for pi, K, Y
318  static const G4double pion_vp_small;
319  static const G4double kaon_vp;
320  static const G4double hyperon_vp;
321 
322  // Neutrons and protons, for computing trajectory placements
325 };
326 
327 #endif // G4NUCLEI_MODEL_HH
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2, G4double nuclearRadius) const
G4double getRadius(G4int izone) const
const G4double R_nucleon
G4InuclElementaryParticle generateQuasiDeuteron(G4int type1, G4int type2, G4int zone) const
std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
std::pair< G4int, G4int > getTypesOfNucleonsInvolved() const
const G4double gammaQDscale
std::vector< G4double > vz
CLHEP::Hep3Vector G4ThreeVector
void generateInteractionPartners(G4CascadParticle &cparticle)
std::vector< G4double > zone_radii
std::vector< G4double > zone_volumes
G4double z
Definition: TRTMaterials.hh:39
G4CascadeInterpolator< 30 > gammaQDinterp
G4double inverseMeanFreePath(const G4CascadParticle &cparticle, const G4InuclElementaryParticle &target, G4int zone=-1)
G4double getFermiKinetic(G4int ip, G4int izone) const
G4bool passFermi(const std::vector< G4InuclElementaryParticle > &particles, G4int zone)
G4bool stillInside(const G4CascadParticle &cparticle)
G4double zoneIntegralGaussian(G4double ur1, G4double ur2, G4double nuclearRadius) const
G4double v1[6]
std::vector< partner > thePartners
G4double nuclei_volume
G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const
G4double a
Definition: TRTMaterials.hh:39
void boundaryTransition(G4CascadParticle &cparticle)
const G4InuclElementaryParticle protonEP
static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0)
G4double ur[7]
G4CascadParticle initializeCascad(G4InuclElementaryParticle *particle)
G4InuclNuclei * theNucleus
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
static const G4double large
const G4double radiusScale2
G4bool isProjectile(const G4CascadParticle &cparticle) const
static const G4double piTimes4thirds
void fillBindingEnergies()
G4CollisionOutput EPCoutput
const G4double radScaleAlpha
G4double getCurrentDensity(G4int ip, G4int izone) const
void fillZoneRadii(G4double nuclearRadius)
G4LorentzConvertor dummy_convertor
static const G4double kaon_vp
G4int protonNumberCurrent
void fillPotentials(G4int type, G4double tot_vol)
static const G4double alfa6[6]
bool G4bool
Definition: G4Types.hh:79
std::vector< G4LorentzVector > momentums
G4double iz
Definition: TRTMaterials.hh:39
std::vector< G4ThreeVector > coordinates
G4double generateInteractionLength(const G4CascadParticle &cparticle, G4double path, G4double invmfp) const
G4double nuclei_radius
std::vector< std::vector< G4double > > nucleon_densities
G4double fillZoneVolumes(G4double nuclearRadius)
void generateModel(G4InuclNuclei *nuclei)
G4double v[6]
std::pair< G4InuclElementaryParticle, G4double > partner
static const G4double pion_vp_small
const G4double crossSectionUnits
std::vector< G4double > binding_energies
const G4double radiusUnits
G4double getRadius() const
std::vector< G4InuclElementaryParticle > qdeutrons
G4int getNumberOfZones() const
static const G4double alfa3[3]
void generateParticleFate(G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)
std::vector< G4double > rod
G4double getRatio(G4int ip) const
G4int getNumberOfProtons() const
G4double totalCrossSection(G4double ke, G4int rtype) const
G4bool empty() const
static const G4double small
G4double getFermiMomentum(G4int ip, G4int izone) const
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)
G4int getCurrentZone() const
G4int getNumberOfNeutrons() const
void setVerboseLevel(G4int verbose)
const G4InuclElementaryParticle neutronEP
std::vector< G4ThreeVector > collisionPts
void printModel() const
std::vector< G4double > acsecs
G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) 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()
G4double getDensity(G4int ip, G4int izone) const
G4double getRadiusUnits() const
G4bool passTrailing(const G4ThreeVector &hit_position)
G4int getZone(G4double r) const
static const G4double pion_vp
const G4double radiusForSmall
const G4double skinDepth
G4bool forceFirst(const G4CascadParticle &cparticle) const
static const double fermi
Definition: G4SIunits.hh:102
CLHEP::HepLorentzVector G4LorentzVector
G4double getPotential(G4int ip, G4int izone) const