Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:21:03

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