Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-03-30 07:50:10

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 /// \file G4ScreenedNuclearRecoil.hh
0027 /// \brief Definition of the G4ScreenedNuclearRecoil class
0028 
0029 // G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp
0030 // GEANT4 tag
0031 //
0032 //
0033 //
0034 // Class Description
0035 // Process for screened electromagnetic nuclear elastic scattering;
0036 // Physics comes from:
0037 // Marcus H. Mendenhall and Robert A. Weller,
0038 // "Algorithms  for  the rapid  computation  of  classical  cross  sections
0039 // for  screened  Coulomb  collisions  "
0040 // Nuclear  Instruments  and  Methods  in  Physics  Research  B58  (1991)  11-17
0041 // The only input required is a screening function phi(r/a) which is the ratio
0042 // of the actual interatomic potential for two atoms with atomic numbers Z1 and
0043 // Z2,
0044 // to the unscreened potential Z1*Z2*e^2/r where e^2 is elm_coupling in Geant4
0045 // units
0046 // the actual screening tables are computed externally in a python module
0047 // "screened_scattering.py"
0048 // to allow very specific screening functions to be added if desired, without
0049 // messing with the insides of this code.
0050 //
0051 // First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
0052 // May 1, 2008 -- Added code to allow process to have zero cross section above
0053 //  max energy, to coordinate with G4MSC.  -- mhm
0054 //
0055 // Class Description - End
0056 
0057 #ifndef G4ScreenedNuclearRecoil_h
0058 #define G4ScreenedNuclearRecoil_h 1
0059 
0060 #include "c2_function.hh"
0061 
0062 #include "G4ParticleChange.hh"
0063 #include "G4PhysicalConstants.hh"
0064 #include "G4SystemOfUnits.hh"
0065 #include "G4VDiscreteProcess.hh"
0066 #include "globals.hh"
0067 
0068 #include <map>
0069 #include <vector>
0070 
0071 class G4VNIELPartition;
0072 
0073 typedef c2_const_ptr<G4double> G4_c2_const_ptr;
0074 typedef c2_ptr<G4double> G4_c2_ptr;
0075 typedef c2_function<G4double> G4_c2_function;
0076 
0077 typedef struct G4ScreeningTables
0078 {
0079     G4double z1, z2, m1, m2, au, emin;
0080     G4_c2_const_ptr EMphiData;
0081 } G4ScreeningTables;
0082 
0083 // A class for loading ScreenedCoulombCrossSections
0084 class G4ScreenedCoulombCrossSectionInfo
0085 {
0086   public:
0087     G4ScreenedCoulombCrossSectionInfo() {}
0088     ~G4ScreenedCoulombCrossSectionInfo() {}
0089 
0090     static const char* CVSHeaderVers()
0091     {
0092       return "G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp GEANT4 tag ";
0093     }
0094     static const char* CVSFileVers();
0095 };
0096 
0097 // A class for loading ScreenedCoulombCrossSections
0098 class G4ScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSectionInfo
0099 {
0100   public:
0101     G4ScreenedCoulombCrossSection() : verbosity(1) {}
0102     G4ScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
0103       : G4ScreenedCoulombCrossSectionInfo(), verbosity(src.verbosity)
0104     {}
0105     virtual ~G4ScreenedCoulombCrossSection();
0106 
0107     typedef std::map<G4int, G4ScreeningTables> ScreeningMap;
0108 
0109     // a local, fast-access mapping of a particle's Z to its full definition
0110     typedef std::map<G4int, class G4ParticleDefinition*> ParticleCache;
0111 
0112     // LoadData is called by G4ScreenedNuclearRecoil::GetMeanFreePath
0113     // It loads the data tables, builds the elemental cross-section tables.
0114     virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff) = 0;
0115 
0116     // BuildMFPTables is called by G4ScreenedNuclearRecoil::GetMeanFreePath
0117     // to build the MFP tables for each material
0118     void BuildMFPTables(void);  // scan the MaterialsTable and construct MFP
0119                                 // tables
0120 
0121     virtual G4ScreenedCoulombCrossSection* create() = 0;
0122     // a 'virtual constructor' which clones the class
0123     const G4ScreeningTables* GetScreening(G4int Z) { return &(screeningData[Z]); }
0124     void SetVerbosity(G4int v) { verbosity = v; }
0125 
0126     // this process needs element selection weighted only by number density
0127     G4ParticleDefinition* SelectRandomUnweightedTarget(const G4MaterialCutsCouple* couple);
0128 
0129     enum
0130     {
0131       nMassMapElements = 116
0132     };
0133 
0134     G4double standardmass(G4int z1) { return z1 <= nMassMapElements ? massmap[z1] : 2.5 * z1; }
0135 
0136     // get the mean-free-path table for the indexed material
0137     const G4_c2_function* operator[](G4int materialIndex)
0138     {
0139       return MFPTables.find(materialIndex) != MFPTables.end() ? &(MFPTables[materialIndex].get())
0140                                                               : (G4_c2_function*)0;
0141     }
0142 
0143   protected:
0144     ScreeningMap screeningData;  // screening tables for each element
0145     ParticleCache targetMap;
0146     G4int verbosity;
0147     std::map<G4int, G4_c2_const_ptr> sigmaMap;
0148     // total cross section for each element
0149     std::map<G4int, G4_c2_const_ptr> MFPTables;  // MFP for each material
0150 
0151   private:
0152     static const G4double massmap[nMassMapElements + 1];
0153 };
0154 
0155 typedef struct G4CoulombKinematicsInfo
0156 {
0157     G4double impactParameter;
0158     G4ScreenedCoulombCrossSection* crossSection;
0159     G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
0160     G4ParticleDefinition* recoilIon;
0161     const G4Material* targetMaterial;
0162 } G4CoulombKinematicsInfo;
0163 
0164 class G4ScreenedCollisionStage
0165 {
0166   public:
0167     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
0168                                  const class G4Step& aStep) = 0;
0169     virtual ~G4ScreenedCollisionStage() {}
0170 };
0171 
0172 class G4ScreenedCoulombClassicalKinematics : public G4ScreenedCoulombCrossSectionInfo,
0173                                              public G4ScreenedCollisionStage
0174 {
0175   public:
0176     G4ScreenedCoulombClassicalKinematics();
0177     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
0178                                  const class G4Step& aStep);
0179 
0180     G4bool DoScreeningComputation(class G4ScreenedNuclearRecoil* master,
0181                                   const G4ScreeningTables* screen, G4double eps, G4double beta);
0182 
0183     virtual ~G4ScreenedCoulombClassicalKinematics() {}
0184 
0185   protected:
0186     // the c2_functions we need to do the work.
0187     c2_const_plugin_function_p<G4double>& phifunc;
0188     c2_linear_p<G4double>& xovereps;
0189     G4_c2_ptr diff;
0190 };
0191 
0192 class G4SingleScatter : public G4ScreenedCoulombCrossSectionInfo, public G4ScreenedCollisionStage
0193 {
0194   public:
0195     G4SingleScatter() {}
0196     virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
0197                                  const class G4Step& aStep);
0198     virtual ~G4SingleScatter() {}
0199 };
0200 
0201 /**
0202      \brief A process which handles screened Coulomb collisions between nuclei
0203 
0204 */
0205 
0206 class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo, public G4VDiscreteProcess
0207 {
0208   public:
0209     friend class G4ScreenedCollisionStage;
0210 
0211     /// \param processName the name to assign the process
0212     /// \param ScreeningKey the name of a screening function to use.
0213     /// The default functions are "zbl" (recommended for soft scattering),
0214     /// "lj" (recommended for backscattering) and "mol" (Moliere potential)
0215     /// \param GenerateRecoils if frue, ions struck by primary are converted
0216     /// into new moving particles.
0217     /// If false, energy is deposited, but no new moving ions are created.
0218     /// \param RecoilCutoff energy below which no new moving particles will be
0219     /// created, even if a GenerateRecoils is true.
0220     /// Also, a moving primary particle will be stopped if its energy falls
0221     /// below this limit.
0222     /// \param PhysicsCutoff the energy transfer to which screening tables are
0223     /// calucalted.
0224     /// There is no really
0225     /// compelling reason to change it from the 10.0 eV default.
0226     /// However, see the paper on running this
0227     /// in thin targets for further discussion, and its interaction
0228     /// with SetMFPScaling()
0229 
0230     G4ScreenedNuclearRecoil(const G4String& processName = "ScreenedElastic",
0231                             const G4String& ScreeningKey = "zbl", G4bool GenerateRecoils = 1,
0232                             G4double RecoilCutoff = 100.0 * eV, G4double PhysicsCutoff = 10.0 * eV);
0233 
0234     virtual ~G4ScreenedNuclearRecoil();
0235 
0236     virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition*);
0237 
0238     virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep);
0239     /// \param aParticleType the particle to test
0240     virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
0241     /// \param aParticleType the type of particle to build tables for
0242     virtual void BuildPhysicsTable(const G4ParticleDefinition& aParticleType);
0243     /// \param aParticleType the type of particle to build tables for
0244     virtual void DumpPhysicsTable(const G4ParticleDefinition& aParticleType);
0245     /// range of the selected nucleus
0246     /// \param A the nucleon number of the beam
0247     /// \param A1 the nucleon number of the target
0248     /// \param apsis the distance of closest approach
0249 
0250     virtual G4bool CheckNuclearCollision(G4double A, G4double A1,
0251                                          G4double apsis);  // return true if hard collision
0252 
0253     virtual G4ScreenedCoulombCrossSection* GetNewCrossSectionHandler(void);
0254 
0255     G4double GetNIEL() const { return NIEL; }
0256 
0257     void ResetTables();
0258     // clear all data tables to allow changing energy cutoff, materials, etc.
0259 
0260     /// cross section
0261     ///
0262     /// This funciton is used to coordinate this process with G4MSC.
0263     /// Typically, G4MSC should
0264     /// not be allowed to operate in a range which overlaps that of this
0265     /// process.  The criterion which is most reasonable
0266     /// is that the transition should be somewhere in the modestly
0267     /// relativistic regime (500 MeV/u for example).
0268     /// param energy energy per nucleon for the cutoff
0269 
0270     void SetMaxEnergyForScattering(G4double energy) { processMaxEnergy = energy; }
0271 
0272     std::string GetScreeningKey() const { return screeningKey; }
0273 
0274     /// \param flag if true, enable deposition of energy (the default).
0275     /// If false, disable deposition.
0276 
0277     void AllowEnergyDeposition(G4bool flag) { registerDepositedEnergy = flag; }
0278 
0279     G4bool GetAllowEnergyDeposition() const { return registerDepositedEnergy; }
0280 
0281     /// If recoils are disabled, the energy they would have received is just
0282     /// deposited.
0283     /// param flag if true, create recoil ions in cases in which the energy
0284     /// is above the recoilCutoff.
0285     /// If false, just deposit the energy.
0286 
0287     void EnableRecoils(G4bool flag) { generateRecoils = flag; }
0288 
0289     G4bool GetEnableRecoils() const { return generateRecoils; }
0290 
0291     /// \param scale the factor by which the default MFP will be scaled.
0292     /// Set to less than 1 for very thin films, typically,
0293     /// to sample multiple scattering,
0294     /// or to greater than 1 for quick simulations with a very long flight path
0295 
0296     void SetMFPScaling(G4double scale) { MFPScale = scale; }
0297 
0298     G4double GetMFPScaling() const { return MFPScale; }
0299 
0300     /// which are close enough they need hadronic phsyics.
0301     /// Default is true (skip close collisions).
0302     /// Disabling this results in excess nuclear stopping power.
0303     /// \param flag true results in hard collisions being skipped.
0304     /// false allows hard collisions.
0305 
0306     void AvoidNuclearReactions(G4bool flag) { avoidReactions = flag; }
0307 
0308     G4bool GetAvoidNuclearReactions() const { return avoidReactions; }
0309 
0310     /// be generated,
0311     /// and the energy (per nucleon) below which all ions are stopped.
0312     /// \param energy energy per nucleon
0313 
0314     void SetRecoilCutoff(G4double energy) { recoilCutoff = energy; }
0315 
0316     G4double GetRecoilCutoff() const { return recoilCutoff; }
0317 
0318     /// Typically, this is 10 eV or so, and not often changed.
0319     /// \param energy the cutoff energy
0320 
0321     void SetPhysicsCutoff(G4double energy)
0322     {
0323       physicsCutoff = energy;
0324       ResetTables();
0325     }
0326 
0327     G4double GetPhysicsCutoff() const { return physicsCutoff; }
0328 
0329 
0330     void SetNIELPartitionFunction(const G4VNIELPartition* part);
0331 
0332     /// backscattering
0333     /// \param fraction the fraction of particles to have their cross section
0334     /// boosted.
0335     /// \param HardeningFactor the factor by which to boost the scattering
0336     /// cross section.
0337 
0338     void SetCrossSectionHardening(G4double fraction, G4double HardeningFactor)
0339     {
0340       hardeningFraction = fraction;
0341       hardeningFactor = HardeningFactor;
0342     }
0343 
0344     G4double GetHardeningFraction() const { return hardeningFraction; }
0345 
0346     G4double GetHardeningFactor() const { return hardeningFactor; }
0347 
0348     G4double GetCurrentInteractionLength() const { return currentInteractionLength; }
0349 
0350     /// if the user needs non-standard behavior.
0351     /// \param cs a class which constructs the screening tables.
0352 
0353     void SetExternalCrossSectionHandler(G4ScreenedCoulombCrossSection* cs)
0354     {
0355       externalCrossSectionConstructor = cs;
0356     }
0357 
0358     G4int GetVerboseLevel() const { return verboseLevel; }
0359 
0360     std::map<G4int, G4ScreenedCoulombCrossSection*>& GetCrossSectionHandlers()
0361     {
0362       return crossSectionHandlers;
0363     }
0364 
0365     void ClearStages(void);
0366 
0367     void AddStage(G4ScreenedCollisionStage* stage) { collisionStages.push_back(stage); }
0368 
0369     G4CoulombKinematicsInfo& GetKinematics() { return kinematics; }
0370 
0371     void SetValidCollision(G4bool flag) { validCollision = flag; }
0372 
0373     G4bool GetValidCollision() const { return validCollision; }
0374 
0375     /// for internal use, primarily.
0376     class G4ParticleChange& GetParticleChange()
0377     {
0378       return static_cast<G4ParticleChange&>(*pParticleChange);
0379     }
0380 
0381     /// to partition it into NIEL and ionizing energy.
0382 
0383     void DepositEnergy(G4int z1, G4double a1, const G4Material* material, G4double energy);
0384 
0385   protected:
0386     G4double highEnergyLimit;
0387 
0388     G4double lowEnergyLimit;
0389 
0390     /// to cross over to G4MSC
0391     G4double processMaxEnergy;
0392     G4String screeningKey;
0393     G4bool generateRecoils, avoidReactions;
0394     G4double recoilCutoff, physicsCutoff;
0395     G4bool registerDepositedEnergy;
0396     G4double IonizingLoss, NIEL;
0397     G4double MFPScale;
0398     G4double hardeningFraction, hardeningFactor;
0399 
0400     G4ScreenedCoulombCrossSection* externalCrossSectionConstructor;
0401     std::vector<G4ScreenedCollisionStage*> collisionStages;
0402 
0403     std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;
0404 
0405     G4bool validCollision;
0406     G4CoulombKinematicsInfo kinematics;
0407     const G4VNIELPartition* NIELPartitionFunction;
0408 };
0409 
0410 // A customized G4CrossSectionHandler which gets its data from
0411 // an external program
0412 
0413 class G4NativeScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSection
0414 {
0415   public:
0416     G4NativeScreenedCoulombCrossSection();
0417 
0418     G4NativeScreenedCoulombCrossSection(const G4NativeScreenedCoulombCrossSection& src)
0419       : G4ScreenedCoulombCrossSection(src), phiMap(src.phiMap)
0420     {}
0421 
0422     G4NativeScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
0423       : G4ScreenedCoulombCrossSection(src)
0424     {}
0425 
0426     virtual ~G4NativeScreenedCoulombCrossSection();
0427 
0428     virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff);
0429 
0430     virtual G4ScreenedCoulombCrossSection* create()
0431     {
0432       return new G4NativeScreenedCoulombCrossSection(*this);
0433     }
0434 
0435     // get a list of available keys
0436     std::vector<G4String> GetScreeningKeys() const;
0437 
0438     typedef G4_c2_function& (*ScreeningFunc)(G4int z1, G4int z2, size_t nPoints, G4double rMax,
0439                                              G4double* au);
0440 
0441     void AddScreeningFunction(G4String name, ScreeningFunc fn) { phiMap[name] = fn; }
0442 
0443   private:
0444     // this is a map used to look up screening function generators
0445     std::map<std::string, ScreeningFunc> phiMap;
0446 };
0447 
0448 G4_c2_function& ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
0449 
0450 G4_c2_function& MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax,
0451                                  G4double* auval);
0452 
0453 G4_c2_function& LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
0454 
0455 G4_c2_function& LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
0456 
0457 #endif