Geant4  10.00.p03
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 
76 #ifndef G4NUCLEI_MODEL_HH
77 #define G4NUCLEI_MODEL_HH
78 
79 #include <algorithm>
80 #include <vector>
81 
83 #include "G4CascadParticle.hh"
84 #include "G4CascadeInterpolator.hh"
85 #include "G4CollisionOutput.hh"
86 #include "G4LorentzConvertor.hh"
87 
88 class G4InuclNuclei;
90 
92 public:
93  G4NucleiModel();
96 
97  virtual ~G4NucleiModel();
98 
99  void setVerboseLevel(G4int verbose) { verboseLevel = verbose; }
100 
101  void generateModel(G4InuclNuclei* nuclei);
102  void generateModel(G4int a, G4int z);
103 
104  // Arguments here are used for rescattering (::Propagate)
105  void reset(G4int nHitNeutrons=0, G4int nHitProtons=0,
106  const std::vector<G4ThreeVector>* hitPoints=0);
107 
108  void printModel() const;
109 
110  G4double getDensity(G4int ip, G4int izone) const {
111  return nucleon_densities[ip - 1][izone];
112  }
113 
115  return fermi_momenta[ip - 1][izone];
116  }
117 
118  G4double getFermiKinetic(G4int ip, G4int izone) const;
119 
120  G4double getPotential(G4int ip, G4int izone) const {
121  if (ip == 9 || ip < 0) return 0.0; // Photons and leptons
122  G4int ip0 = ip < 3 ? ip - 1 : 2;
123  if (ip > 10 && ip < 18) ip0 = 3;
124  if (ip > 20) ip0 = 4;
125  return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
126  }
127 
128  // Factor to convert GEANT4 lengths to internal units
130 
131  G4double getRadius() const { return nuclei_radius; }
132  G4double getRadius(G4int izone) const {
133  return ( (izone<0) ? 0.
134  : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
135  }
136  G4double getVolume(G4int izone) const {
137  return ( (izone<0) ? 0.
138  : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
139  }
140 
142  G4int getZone(G4double r) const {
143  for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
144  return number_of_zones;
145  }
146 
149 
150  G4bool empty() const {
151  return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
152  }
153 
154  G4bool stillInside(const G4CascadParticle& cparticle) {
155  return cparticle.getCurrentZone() < number_of_zones;
156  }
157 
158 
160 
161  typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > modelLists;
162 
163  void initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target,
164  modelLists& output);
165 
166 
167  std::pair<G4int, G4int> getTypesOfNucleonsInvolved() const {
168  return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
169  }
170 
171  void generateParticleFate(G4CascadParticle& cparticle,
172  G4ElementaryParticleCollider* theEPCollider,
173  std::vector<G4CascadParticle>& cascade);
174 
175  G4bool forceFirst(const G4CascadParticle& cparticle) const;
176  G4bool isProjectile(const G4CascadParticle& cparticle) const;
177  G4bool worthToPropagate(const G4CascadParticle& cparticle) const;
178 
180 
182 
184  G4double totalCrossSection(G4double ke, G4int rtype) const;
185 
186  // Identify whether given particle can interact with dibaryons
187  static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0);
188 
189 protected:
190  G4bool passFermi(const std::vector<G4InuclElementaryParticle>& particles,
191  G4int zone);
192 
193  G4bool passTrailing(const G4ThreeVector& hit_position);
194 
195  void boundaryTransition(G4CascadParticle& cparticle);
196 
197  void choosePointAlongTraj(G4CascadParticle& cparticle);
198 
200  G4int type2,
201  G4int zone) const;
202 
203  typedef std::pair<G4InuclElementaryParticle, G4double> partner;
204 
205  std::vector<partner> thePartners; // Buffer for output below
207 
208  // Function for std::sort() to use in organizing partners by path length
209  static G4bool sortPartners(const partner& p1, const partner& p2) {
210  return (p2.second > p1.second);
211  }
212 
213  // Functions used to generate model nuclear structure
214  void fillBindingEnergies();
215 
216  void fillZoneRadii(G4double nuclearRadius);
217 
218  G4double fillZoneVolumes(G4double nuclearRadius);
219 
220  void fillPotentials(G4int type, G4double tot_vol);
221 
223  G4double nuclearRadius) const;
224 
226  G4double nuclearRadius) const;
227 
228  G4double getRatio(G4int ip) const; // Fraction of nucleons still available
229 
230  // Scale nuclear density with interactions so far
231  G4double getCurrentDensity(G4int ip, G4int izone) const;
232 
233  // Average path length for any interaction of projectile and target
234  // = 1. / (density * cross-section)
236  const G4InuclElementaryParticle& target,
237  G4int zone = -1); // Override zone value
238  // NOTE: Function is non-const in order to use dummy_converter
239 
240  // Use path-length and MFP (above) to throw random distance to collision
242  G4double path, G4double invmfp) const;
243 
244 private:
246 
247  // Buffers for processing interactions on each cycle
248  G4LorentzConvertor dummy_convertor; // For getting collision frame
249  G4CollisionOutput EPCoutput; // For N-body inelastic collisions
250 
251  std::vector<G4InuclElementaryParticle> qdeutrons; // For h+(NN) trials
252  std::vector<G4double> acsecs;
253 
254  std::vector<G4ThreeVector> coordinates; // for initializeCascad()
255  std::vector<G4LorentzVector> momentums;
256  std::vector<G4InuclElementaryParticle> raw_particles;
257 
258  std::vector<G4ThreeVector> collisionPts;
259 
260  // Temporary buffers for computing nuclear model
261  G4double ur[7]; // Number of skin depths for each zone
262  G4double v[6]; // Density integrals by zone
263  G4double v1[6]; // Pseudo-volume (delta r^3) by zone
264  std::vector<G4double> rod; // Nucleon density
265  std::vector<G4double> pf; // Fermi momentum
266  std::vector<G4double> vz; // Nucleon potential
267 
268  // Nuclear model configuration
269  std::vector<std::vector<G4double> > nucleon_densities;
270  std::vector<std::vector<G4double> > zone_potentials;
271  std::vector<std::vector<G4double> > fermi_momenta;
272  std::vector<G4double> zone_radii;
273  std::vector<G4double> zone_volumes;
274  std::vector<G4double> binding_energies;
278 
282 
285 
288 
291 
292  G4CascadeInterpolator<30> gammaQDinterp; // quasideuteron interpolator
293 
294  // Symbolic names for nuclear potentials
296 
297  // Cutoffs for extreme values
298  static const G4double small;
299  static const G4double large;
300 
301  // Parameters for nuclear structure
302  static const G4double skinDepth;
303  static const G4double radiusScale; // Coefficients for two-parameter fit
304  static const G4double radiusScale2; // R = 1.16*cbrt(A) - 1.3456/cbrt(A)
305  static const G4double radiusForSmall; // Average radius of light A<5 nuclei
306  static const G4double radScaleAlpha; // Scaling factor R_alpha/R_small
307  static const G4double fermiMomentum;
308  static const G4double R_nucleon;
309  static const G4double alfa3[3], alfa6[6];
310  static const G4double pion_vp;
311  static const G4double pion_vp_small;
312  static const G4double kaon_vp;
313  static const G4double hyperon_vp;
314 
315  // Scale parameter for gamma-quasideuteron scattering
316  static const G4double gammaQDscale;
317 
318  // FIXME: We should not be using this!
319  static const G4double piTimes4thirds;
321  static const G4double radiusUnits;
322 
323  // Neutrons and protons, for computing trajectory placements
326 };
327 
328 #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
static 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
static 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
static const G4double radiusScale
G4double absorptionCrossSection(G4double e, G4int type) const
static const G4double large
static const G4double radiusScale2
G4bool isProjectile(const G4CascadParticle &cparticle) const
static const G4double piTimes4thirds
void fillBindingEnergies()
G4CollisionOutput EPCoutput
static 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
static const G4double crossSectionUnits
std::vector< G4double > binding_energies
static 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
static 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
static const G4double radiusForSmall
static const G4double skinDepth
G4bool forceFirst(const G4CascadParticle &cparticle) const
static const double fermi
Definition: G4SIunits.hh:93
CLHEP::HepLorentzVector G4LorentzVector
G4double getPotential(G4int ip, G4int izone) const