include/Geant4/G4NucleiModel.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // 20100319  M. Kelsey -- Remove "using" directory and unnecessary #includes,
0028 //      move ctor to .cc file
0029 // 20100407  M. Kelsey -- Create "partners thePartners" data member to act
0030 //      as buffer between ::generateInteractionPartners() and
0031 //      ::generateParticleFate(), and make "outgoing_cparticles" a
0032 //      data member returned from the latter by const-ref.  Replace
0033 //      return-by-value of initializeCascad() with an input buffer.
0034 // 20100409  M. Kelsey -- Add function to sort list of partnerts by pathlen,
0035 //      move non-inlinable code to .cc.
0036 // 20100421  M. Kelsey -- Move getFermiKinetic() to .cc, no hardwired masses.
0037 // 20100517  M. Kelsey -- Change cross-section tables to static arrays.  Move
0038 //      absorptionCrossSection() from SpecialFunc.
0039 // 20100520  M. Kelsey -- Add function to separate momentum from nucleon
0040 // 20100617  M. Kelsey -- Add setVerboseLevel() function, add generateModel()
0041 //      with particle input, and ctor with A/Z input.
0042 // 20100715  M. Kelsey -- Add G4InuclNuclei object for use with balance checks
0043 // 20100723  M. Kelsey -- Move G4CollisionOutput buffer, along with all
0044 //      std::vector<> buffers, here for reuse
0045 // 20100914  M. Kelsey -- Migrate to integer A and Z
0046 // 20101004  M. Kelsey -- Rename and create functions used to generate model
0047 // 20101019  M. Kelsey -- CoVerity report: dtor leak; move dtor to .cc file
0048 // 20110223  M. Kelsey -- Add static parameters for radius and cross-section
0049 //      scaling factors.
0050 // 20110303  M. Kelsey -- Add accessors for model parameters and units
0051 // 20110304  M. Kelsey -- Extend reset() to fill neutron and proton counts
0052 // 20110324  D. Wright -- Add list of nucleon interaction points, and nucleon
0053 //      effective radius, for trailing effect.
0054 // 20110324  M. Kelsey -- Extend reset() to pass list of points; move
0055 //      implementation to .cc file.
0056 // 20110405  M. Kelsey -- Add "passTrailing()" to encapsulate trailing effect
0057 // 20110808  M. Kelsey -- Pass buffer into generateParticleFate instead of
0058 //      returning a vector to be copied.
0059 // 20110823  M. Kelsey -- Remove local cross-section tables entirely
0060 // 20120125  M. Kelsey -- Add special case for photons to have zero potential.
0061 // 20120307  M. Kelsey -- Add zone volume array for use with quasideuteron
0062 //      density, functions to identify projectile, compute interaction
0063 //      distance.
0064 // 20130129  M. Kelsey -- Add static objects for gamma-quasideuteron scattering
0065 // 20130221  M. Kelsey -- Add function to emplant particle along trajectory
0066 // 20130226  M. Kelsey -- Allow forcing zone selection in MFP calculation.
0067 // 20130302  M. Kelsey -- Add forceFirst() wrapper function (allows configuring)
0068 // 20130628  M. Kelsey -- Extend useQuasiDeuteron() to check good absorption,
0069 //      fix spelling of "Deutron" -> "Deuteron"
0070 // 20130808  M. Kelsey -- To avoid thread collisions, move static neutronEP
0071 //      and protonEP objects to const data members.
0072 // 20131001  M. Kelsey -- Move QDinterp object to data member, thread local
0073 // 20140116  M. Kelsey -- Move statics to const data members to avoid weird
0074 //      interactions with MT.
0075
0076 #ifndef G4NUCLEI_MODEL_HH
0077 #define G4NUCLEI_MODEL_HH
0078
0079 #include <algorithm>
0080 #include <vector>
0081
0082 #include "G4InuclElementaryParticle.hh"
0083 #include "G4CascadParticle.hh"
0084 #include "G4CascadeInterpolator.hh"
0085 #include "G4CollisionOutput.hh"
0086 #include "G4LorentzConvertor.hh"
0087
0088 class G4InuclNuclei;
0089 class G4ElementaryParticleCollider;
0090
0091 class G4NucleiModel {
0092 public:
0093   G4NucleiModel();
0094   G4NucleiModel(G4int a, G4int z);
0095   explicit G4NucleiModel(G4InuclNuclei* nuclei);
0096
0097   virtual ~G4NucleiModel();
0098
0099   void setVerboseLevel(G4int verbose) { verboseLevel = verbose; }
0100
0101   void generateModel(G4InuclNuclei* nuclei);
0102   void generateModel(G4int a, G4int z);
0103
0104   // Arguments here are used for rescattering (::Propagate)
0105   void reset(G4int nHitNeutrons=0, G4int nHitProtons=0,
0106          const std::vector<G4ThreeVector>* hitPoints=0);
0107
0108   void printModel() const;
0109
0110   G4double getDensity(G4int ip, G4int izone) const {
0111     return nucleon_densities[ip - 1][izone];
0112   }
0113
0114   G4double getFermiMomentum(G4int ip, G4int izone) const {
0115     return fermi_momenta[ip - 1][izone];
0116   }
0117
0118   G4double getFermiKinetic(G4int ip, G4int izone) const;
0119
0120   G4double getPotential(G4int ip, G4int izone) const {
0121     if (ip == 9 || ip < 0) return 0.0;      // Photons and leptons
0122     G4int ip0 = ip < 3 ? ip - 1 : 2;
0123     if (ip > 10 && ip < 18) ip0 = 3;
0124     if (ip > 20) ip0 = 4;
0125     return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
0126   }
0127
0128   // Factor to convert GEANT4 lengths to internal units
0129   G4double getRadiusUnits() const { return radiusUnits*CLHEP::fermi; }
0130
0131   G4double getRadius() const { return nuclei_radius; }
0132   G4double getRadius(G4int izone) const {
0133     return ( (izone<0) ? 0.
0134          : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
0135   }
0136   G4double getVolume(G4int izone) const {
0137     return ( (izone<0) ? 0.
0138          : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
0139   }
0140
0141   G4int getNumberOfZones() const { return number_of_zones; }
0142   G4int getZone(G4double r) const {
0143     for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
0144     return number_of_zones;
0145   }
0146
0147   G4int getNumberOfNeutrons() const { return neutronNumberCurrent; }
0148   G4int getNumberOfProtons() const  { return protonNumberCurrent; }
0149
0150   G4bool empty() const {
0151     return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
0152   }
0153
0154   G4bool stillInside(const G4CascadParticle& cparticle) {
0155     return cparticle.getCurrentZone() < number_of_zones;
0156   }
0157
0158
0159   G4CascadParticle initializeCascad(G4InuclElementaryParticle* particle);
0160
0161   typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > modelLists;
0162
0163   void initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target,
0164             modelLists& output);
0165
0166
0167   std::pair<G4int, G4int> getTypesOfNucleonsInvolved() const {
0168     return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
0169   }
0170
0171   void generateParticleFate(G4CascadParticle& cparticle,
0172                 G4ElementaryParticleCollider* theEPCollider,
0173                 std::vector<G4CascadParticle>& cascade);
0174
0175   G4bool forceFirst(const G4CascadParticle& cparticle) const;
0176   G4bool isProjectile(const G4CascadParticle& cparticle) const;
0177   G4bool worthToPropagate(const G4CascadParticle& cparticle) const;
0178
0179   G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const;
0180
0181   G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const;
0182
0183   G4double absorptionCrossSection(G4double e, G4int type) const;
0184   G4double totalCrossSection(G4double ke, G4int rtype) const;
0185
0186   // Identify whether given particle can interact with dibaryons
0187   static G4bool useQuasiDeuteron(G4int ptype, G4int qdtype=0);
0188
0189 protected:
0190   G4bool passFermi(const std::vector<G4InuclElementaryParticle>& particles,
0191            G4int zone);
0192
0193   G4bool passTrailing(const G4ThreeVector& hit_position);
0194
0195   void boundaryTransition(G4CascadParticle& cparticle);
0196
0197   void choosePointAlongTraj(G4CascadParticle& cparticle);
0198
0199   G4InuclElementaryParticle generateQuasiDeuteron(G4int type1,
0200                           G4int type2,
0201                           G4int zone) const;
0202
0203   typedef std::pair<G4InuclElementaryParticle, G4double> partner;
0204
0205   std::vector<partner> thePartners;     // Buffer for output below
0206   void generateInteractionPartners(G4CascadParticle& cparticle);
0207
0208   // Function for std::sort() to use in organizing partners by path length
0209   static G4bool sortPartners(const partner& p1, const partner& p2) {
0210     return (p2.second > p1.second);
0211   }
0212
0213   // Functions used to generate model nuclear structure
0214   void fillBindingEnergies();
0215
0216   void fillZoneRadii(G4double nuclearRadius);
0217
0218   G4double fillZoneVolumes(G4double nuclearRadius);
0219
0220   void fillPotentials(G4int type, G4double tot_vol);
0221
0222   G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2,
0223                   G4double nuclearRadius) const;
0224
0225   G4double zoneIntegralGaussian(G4double ur1, G4double ur2,
0226                 G4double nuclearRadius) const;
0227
0228   G4double getRatio(G4int ip) const;    // Fraction of nucleons still available
0229
0230   // Scale nuclear density with interactions so far
0231   G4double getCurrentDensity(G4int ip, G4int izone) const;
0232
0233   // Average path length for any interaction of projectile and target
0234   //    = 1. / (density * cross-section)
0235   G4double inverseMeanFreePath(const G4CascadParticle& cparticle,
0236                    const G4InuclElementaryParticle& target,
0237                    G4int zone = -1);    // Override zone value
0238   // NOTE:  Function is non-const in order to use dummy_converter
0239
0240   // Use path-length and MFP (above) to throw random distance to collision
0241   G4double generateInteractionLength(const G4CascadParticle& cparticle,
0242                      G4double path, G4double invmfp) const;
0243
0244 private:
0245   G4int verboseLevel;
0246
0247   // Buffers for processing interactions on each cycle
0248   G4LorentzConvertor dummy_convertor;   // For getting collision frame
0249   G4CollisionOutput EPCoutput;      // For N-body inelastic collisions
0250
0251   std::vector<G4InuclElementaryParticle> qdeutrons; // For h+(NN) trials
0252   std::vector<G4double> acsecs;
0253
0254   std::vector<G4ThreeVector> coordinates;   // for initializeCascad()
0255   std::vector<G4LorentzVector> momentums;
0256   std::vector<G4InuclElementaryParticle> raw_particles;
0257
0258   std::vector<G4ThreeVector> collisionPts;
0259
0260   // Temporary buffers for computing nuclear model
0261   G4double ur[7];       // Number of skin depths for each zone
0262   G4double v[6];        // Density integrals by zone
0263   G4double v1[6];       // Pseudo-volume (delta r^3) by zone
0264   std::vector<G4double> rod;    // Nucleon density
0265   std::vector<G4double> pf; // Fermi momentum
0266   std::vector<G4double> vz; // Nucleon potential
0267
0268   // Nuclear model configuration
0269   std::vector<std::vector<G4double> > nucleon_densities;
0270   std::vector<std::vector<G4double> > zone_potentials;
0271   std::vector<std::vector<G4double> > fermi_momenta;
0272   std::vector<G4double> zone_radii;
0273   std::vector<G4double> zone_volumes;
0274   std::vector<G4double> binding_energies;
0275   G4double nuclei_radius;
0276   G4double nuclei_volume;
0277   G4int number_of_zones;
0278
0279   G4int A;
0280   G4int Z;
0281   G4InuclNuclei* theNucleus;
0282
0283   G4int neutronNumber;
0284   G4int protonNumber;
0285
0286   G4int neutronNumberCurrent;
0287   G4int protonNumberCurrent;
0288
0289   G4int current_nucl1;
0290   G4int current_nucl2;
0291
0292   G4CascadeInterpolator<30> gammaQDinterp;  // quasideuteron interpolator
0293
0294   // Symbolic names for nuclear potentials
0295   enum PotentialType { WoodsSaxon=0, Gaussian=1 };
0296
0297   // Parameters for nuclear structure
0298   const G4double crossSectionUnits; // Scale from internal to natural units
0299   const G4double radiusUnits;
0300   const G4double skinDepth;     // Fraction of radius for outer skin
0301   const G4double radiusScale;       // Coefficients for two-parameter fit
0302   const G4double radiusScale2;      //   R = 1.16*cbrt(A) - 1.3456/cbrt(A)
0303   const G4double radiusForSmall;    // Average radius of light A<5 nuclei
0304   const G4double radScaleAlpha;     // Scaling factor R_alpha/R_small
0305   const G4double fermiMomentum;
0306   const G4double R_nucleon;
0307   const G4double gammaQDscale;      // Gamma/cluster scattering rescaling
0308   const G4double potentialThickness;    // Defaulted to 1.0
0309
0310   // Cutoffs for extreme values
0311   static const G4double small;
0312   static const G4double large;
0313   static const G4double piTimes4thirds;  // FIXME:  We should not be using this!
0314
0315   static const G4double alfa3[3], alfa6[6]; // Zone boundaries in radii
0316   static const G4double pion_vp;        // Flat potentials for pi, K, Y
0317   static const G4double pion_vp_small;
0318   static const G4double kaon_vp;
0319   static const G4double hyperon_vp;
0320
0321   // Neutrons and protons, for computing trajectory placements
0322   const G4InuclElementaryParticle neutronEP;
0323   const G4InuclElementaryParticle protonEP;
0324 };
0325
0326 #endif // G4NUCLEI_MODEL_HH