Gabriele Cosmo
458 lines
17 KiB
C++
458 lines
17 KiB
C++
//
|
|
// ********************************************************************
|
|
// * License and Disclaimer *
|
|
// * *
|
|
// * The Geant4 software is copyright of the Copyright Holders of *
|
|
// * the Geant4 Collaboration. It is provided under the terms and *
|
|
// * conditions of the Geant4 Software License, included in the file *
|
|
// * LICENSE and available at http://cern.ch/geant4/license . These *
|
|
// * include a list of copyright holders. *
|
|
// * *
|
|
// * Neither the authors of this software system, nor their employing *
|
|
// * institutes,nor the agencies providing financial support for this *
|
|
// * work make any representation or warranty, express or implied, *
|
|
// * regarding this software system or assume any liability for its *
|
|
// * use. Please see the license in the file LICENSE and URL above *
|
|
// * for the full disclaimer and the limitation of liability. *
|
|
// * *
|
|
// * This code implementation is the result of the scientific and *
|
|
// * technical work of the GEANT4 collaboration. *
|
|
// * By using, copying, modifying or distributing the software (or *
|
|
// * any work based on the software) you agree to acknowledge its *
|
|
// * use in resulting scientific publications, and indicate your *
|
|
// * acceptance of all terms of the Geant4 Software license. *
|
|
// ********************************************************************
|
|
//
|
|
/// \file G4ScreenedNuclearRecoil.hh
|
|
/// \brief Definition of the G4ScreenedNuclearRecoil class
|
|
|
|
// G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp
|
|
// GEANT4 tag
|
|
//
|
|
//
|
|
//
|
|
// Class Description
|
|
// Process for screened electromagnetic nuclear elastic scattering;
|
|
// Physics comes from:
|
|
// Marcus H. Mendenhall and Robert A. Weller,
|
|
// "Algorithms for the rapid computation of classical cross sections
|
|
// for screened Coulomb collisions "
|
|
// Nuclear Instruments and Methods in Physics Research B58 (1991) 11-17
|
|
// The only input required is a screening function phi(r/a) which is the ratio
|
|
// of the actual interatomic potential for two atoms with atomic numbers Z1 and
|
|
// Z2,
|
|
// to the unscreened potential Z1*Z2*e^2/r where e^2 is elm_coupling in Geant4
|
|
// units
|
|
// the actual screening tables are computed externally in a python module
|
|
// "screened_scattering.py"
|
|
// to allow very specific screening functions to be added if desired, without
|
|
// messing with the insides of this code.
|
|
//
|
|
// First version, April 2004, Marcus H. Mendenhall, Vanderbilt University
|
|
// May 1, 2008 -- Added code to allow process to have zero cross section above
|
|
// max energy, to coordinate with G4MSC. -- mhm
|
|
//
|
|
// Class Description - End
|
|
|
|
#ifndef G4ScreenedNuclearRecoil_h
|
|
#define G4ScreenedNuclearRecoil_h 1
|
|
|
|
#include "c2_function.hh"
|
|
|
|
#include "G4ParticleChange.hh"
|
|
#include "G4PhysicalConstants.hh"
|
|
#include "G4SystemOfUnits.hh"
|
|
#include "G4VDiscreteProcess.hh"
|
|
#include "globals.hh"
|
|
|
|
#include <map>
|
|
#include <vector>
|
|
|
|
class G4VNIELPartition;
|
|
|
|
typedef c2_const_ptr<G4double> G4_c2_const_ptr;
|
|
typedef c2_ptr<G4double> G4_c2_ptr;
|
|
typedef c2_function<G4double> G4_c2_function;
|
|
|
|
typedef struct G4ScreeningTables
|
|
{
|
|
G4double z1, z2, m1, m2, au, emin;
|
|
G4_c2_const_ptr EMphiData;
|
|
} G4ScreeningTables;
|
|
|
|
// A class for loading ScreenedCoulombCrossSections
|
|
class G4ScreenedCoulombCrossSectionInfo
|
|
{
|
|
public:
|
|
G4ScreenedCoulombCrossSectionInfo() {}
|
|
~G4ScreenedCoulombCrossSectionInfo() {}
|
|
|
|
static const char* CVSHeaderVers()
|
|
{
|
|
return "G4ScreenedNuclearRecoil.hh,v 1.24 2008/05/01 19:58:59 marcus Exp GEANT4 tag ";
|
|
}
|
|
static const char* CVSFileVers();
|
|
};
|
|
|
|
// A class for loading ScreenedCoulombCrossSections
|
|
class G4ScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSectionInfo
|
|
{
|
|
public:
|
|
G4ScreenedCoulombCrossSection() : verbosity(1) {}
|
|
G4ScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
|
|
: G4ScreenedCoulombCrossSectionInfo(), verbosity(src.verbosity)
|
|
{}
|
|
virtual ~G4ScreenedCoulombCrossSection();
|
|
|
|
typedef std::map<G4int, G4ScreeningTables> ScreeningMap;
|
|
|
|
// a local, fast-access mapping of a particle's Z to its full definition
|
|
typedef std::map<G4int, class G4ParticleDefinition*> ParticleCache;
|
|
|
|
// LoadData is called by G4ScreenedNuclearRecoil::GetMeanFreePath
|
|
// It loads the data tables, builds the elemental cross-section tables.
|
|
virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff) = 0;
|
|
|
|
// BuildMFPTables is called by G4ScreenedNuclearRecoil::GetMeanFreePath
|
|
// to build the MFP tables for each material
|
|
void BuildMFPTables(void); // scan the MaterialsTable and construct MFP
|
|
// tables
|
|
|
|
virtual G4ScreenedCoulombCrossSection* create() = 0;
|
|
// a 'virtual constructor' which clones the class
|
|
const G4ScreeningTables* GetScreening(G4int Z) { return &(screeningData[Z]); }
|
|
void SetVerbosity(G4int v) { verbosity = v; }
|
|
|
|
// this process needs element selection weighted only by number density
|
|
G4ParticleDefinition* SelectRandomUnweightedTarget(const G4MaterialCutsCouple* couple);
|
|
|
|
enum
|
|
{
|
|
nMassMapElements = 116
|
|
};
|
|
|
|
G4double standardmass(G4int z1) { return z1 <= nMassMapElements ? massmap[z1] : 2.5 * z1; }
|
|
|
|
// get the mean-free-path table for the indexed material
|
|
const G4_c2_function* operator[](G4int materialIndex)
|
|
{
|
|
return MFPTables.find(materialIndex) != MFPTables.end() ? &(MFPTables[materialIndex].get())
|
|
: (G4_c2_function*)0;
|
|
}
|
|
|
|
protected:
|
|
ScreeningMap screeningData; // screening tables for each element
|
|
ParticleCache targetMap;
|
|
G4int verbosity;
|
|
std::map<G4int, G4_c2_const_ptr> sigmaMap;
|
|
// total cross section for each element
|
|
std::map<G4int, G4_c2_const_ptr> MFPTables; // MFP for each material
|
|
|
|
private:
|
|
static const G4double massmap[nMassMapElements + 1];
|
|
};
|
|
|
|
typedef struct G4CoulombKinematicsInfo
|
|
{
|
|
G4double impactParameter;
|
|
G4ScreenedCoulombCrossSection* crossSection;
|
|
G4double a1, a2, sinTheta, cosTheta, sinZeta, cosZeta, eRecoil;
|
|
G4ParticleDefinition* recoilIon;
|
|
const G4Material* targetMaterial;
|
|
} G4CoulombKinematicsInfo;
|
|
|
|
class G4ScreenedCollisionStage
|
|
{
|
|
public:
|
|
virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
|
|
const class G4Step& aStep) = 0;
|
|
virtual ~G4ScreenedCollisionStage() {}
|
|
};
|
|
|
|
class G4ScreenedCoulombClassicalKinematics : public G4ScreenedCoulombCrossSectionInfo,
|
|
public G4ScreenedCollisionStage
|
|
{
|
|
public:
|
|
G4ScreenedCoulombClassicalKinematics();
|
|
virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
|
|
const class G4Step& aStep);
|
|
|
|
G4bool DoScreeningComputation(class G4ScreenedNuclearRecoil* master,
|
|
const G4ScreeningTables* screen, G4double eps, G4double beta);
|
|
|
|
virtual ~G4ScreenedCoulombClassicalKinematics() {}
|
|
|
|
protected:
|
|
// the c2_functions we need to do the work.
|
|
c2_const_plugin_function_p<G4double>& phifunc;
|
|
c2_linear_p<G4double>& xovereps;
|
|
G4_c2_ptr diff;
|
|
};
|
|
|
|
class G4SingleScatter : public G4ScreenedCoulombCrossSectionInfo, public G4ScreenedCollisionStage
|
|
{
|
|
public:
|
|
G4SingleScatter() {}
|
|
virtual void DoCollisionStep(class G4ScreenedNuclearRecoil* master, const class G4Track& aTrack,
|
|
const class G4Step& aStep);
|
|
virtual ~G4SingleScatter() {}
|
|
};
|
|
|
|
/**
|
|
\brief A process which handles screened Coulomb collisions between nuclei
|
|
|
|
*/
|
|
|
|
class G4ScreenedNuclearRecoil : public G4ScreenedCoulombCrossSectionInfo, public G4VDiscreteProcess
|
|
{
|
|
public:
|
|
friend class G4ScreenedCollisionStage;
|
|
|
|
/// \param processName the name to assign the process
|
|
/// \param ScreeningKey the name of a screening function to use.
|
|
/// The default functions are "zbl" (recommended for soft scattering),
|
|
/// "lj" (recommended for backscattering) and "mol" (Moliere potential)
|
|
/// \param GenerateRecoils if frue, ions struck by primary are converted
|
|
/// into new moving particles.
|
|
/// If false, energy is deposited, but no new moving ions are created.
|
|
/// \param RecoilCutoff energy below which no new moving particles will be
|
|
/// created, even if a GenerateRecoils is true.
|
|
/// Also, a moving primary particle will be stopped if its energy falls
|
|
/// below this limit.
|
|
/// \param PhysicsCutoff the energy transfer to which screening tables are
|
|
/// calucalted.
|
|
/// There is no really
|
|
/// compelling reason to change it from the 10.0 eV default.
|
|
/// However, see the paper on running this
|
|
/// in thin targets for further discussion, and its interaction
|
|
/// with SetMFPScaling()
|
|
|
|
G4ScreenedNuclearRecoil(const G4String& processName = "ScreenedElastic",
|
|
const G4String& ScreeningKey = "zbl", G4bool GenerateRecoils = 1,
|
|
G4double RecoilCutoff = 100.0 * eV, G4double PhysicsCutoff = 10.0 * eV);
|
|
|
|
virtual ~G4ScreenedNuclearRecoil();
|
|
|
|
virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition*);
|
|
|
|
virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep);
|
|
/// \param aParticleType the particle to test
|
|
virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
|
|
/// \param aParticleType the type of particle to build tables for
|
|
virtual void BuildPhysicsTable(const G4ParticleDefinition& aParticleType);
|
|
/// \param aParticleType the type of particle to build tables for
|
|
virtual void DumpPhysicsTable(const G4ParticleDefinition& aParticleType);
|
|
/// range of the selected nucleus
|
|
/// \param A the nucleon number of the beam
|
|
/// \param A1 the nucleon number of the target
|
|
/// \param apsis the distance of closest approach
|
|
|
|
virtual G4bool CheckNuclearCollision(G4double A, G4double A1,
|
|
G4double apsis); // return true if hard collision
|
|
|
|
virtual G4ScreenedCoulombCrossSection* GetNewCrossSectionHandler(void);
|
|
|
|
G4double GetNIEL() const { return NIEL; }
|
|
|
|
void ResetTables();
|
|
// clear all data tables to allow changing energy cutoff, materials, etc.
|
|
|
|
/// cross section
|
|
///
|
|
/// This funciton is used to coordinate this process with G4MSC.
|
|
/// Typically, G4MSC should
|
|
/// not be allowed to operate in a range which overlaps that of this
|
|
/// process. The criterion which is most reasonable
|
|
/// is that the transition should be somewhere in the modestly
|
|
/// relativistic regime (500 MeV/u for example).
|
|
/// param energy energy per nucleon for the cutoff
|
|
|
|
void SetMaxEnergyForScattering(G4double energy) { processMaxEnergy = energy; }
|
|
|
|
std::string GetScreeningKey() const { return screeningKey; }
|
|
|
|
/// \param flag if true, enable deposition of energy (the default).
|
|
/// If false, disable deposition.
|
|
|
|
void AllowEnergyDeposition(G4bool flag) { registerDepositedEnergy = flag; }
|
|
|
|
G4bool GetAllowEnergyDeposition() const { return registerDepositedEnergy; }
|
|
|
|
/// If recoils are disabled, the energy they would have received is just
|
|
/// deposited.
|
|
/// param flag if true, create recoil ions in cases in which the energy
|
|
/// is above the recoilCutoff.
|
|
/// If false, just deposit the energy.
|
|
|
|
void EnableRecoils(G4bool flag) { generateRecoils = flag; }
|
|
|
|
G4bool GetEnableRecoils() const { return generateRecoils; }
|
|
|
|
/// \param scale the factor by which the default MFP will be scaled.
|
|
/// Set to less than 1 for very thin films, typically,
|
|
/// to sample multiple scattering,
|
|
/// or to greater than 1 for quick simulations with a very long flight path
|
|
|
|
void SetMFPScaling(G4double scale) { MFPScale = scale; }
|
|
|
|
G4double GetMFPScaling() const { return MFPScale; }
|
|
|
|
/// which are close enough they need hadronic phsyics.
|
|
/// Default is true (skip close collisions).
|
|
/// Disabling this results in excess nuclear stopping power.
|
|
/// \param flag true results in hard collisions being skipped.
|
|
/// false allows hard collisions.
|
|
|
|
void AvoidNuclearReactions(G4bool flag) { avoidReactions = flag; }
|
|
|
|
G4bool GetAvoidNuclearReactions() const { return avoidReactions; }
|
|
|
|
/// be generated,
|
|
/// and the energy (per nucleon) below which all ions are stopped.
|
|
/// \param energy energy per nucleon
|
|
|
|
void SetRecoilCutoff(G4double energy) { recoilCutoff = energy; }
|
|
|
|
G4double GetRecoilCutoff() const { return recoilCutoff; }
|
|
|
|
/// Typically, this is 10 eV or so, and not often changed.
|
|
/// \param energy the cutoff energy
|
|
|
|
void SetPhysicsCutoff(G4double energy)
|
|
{
|
|
physicsCutoff = energy;
|
|
ResetTables();
|
|
}
|
|
|
|
G4double GetPhysicsCutoff() const { return physicsCutoff; }
|
|
|
|
|
|
void SetNIELPartitionFunction(const G4VNIELPartition* part);
|
|
|
|
/// backscattering
|
|
/// \param fraction the fraction of particles to have their cross section
|
|
/// boosted.
|
|
/// \param HardeningFactor the factor by which to boost the scattering
|
|
/// cross section.
|
|
|
|
void SetCrossSectionHardening(G4double fraction, G4double HardeningFactor)
|
|
{
|
|
hardeningFraction = fraction;
|
|
hardeningFactor = HardeningFactor;
|
|
}
|
|
|
|
G4double GetHardeningFraction() const { return hardeningFraction; }
|
|
|
|
G4double GetHardeningFactor() const { return hardeningFactor; }
|
|
|
|
G4double GetCurrentInteractionLength() const { return currentInteractionLength; }
|
|
|
|
/// if the user needs non-standard behavior.
|
|
/// \param cs a class which constructs the screening tables.
|
|
|
|
void SetExternalCrossSectionHandler(G4ScreenedCoulombCrossSection* cs)
|
|
{
|
|
externalCrossSectionConstructor = cs;
|
|
}
|
|
|
|
G4int GetVerboseLevel() const { return verboseLevel; }
|
|
|
|
std::map<G4int, G4ScreenedCoulombCrossSection*>& GetCrossSectionHandlers()
|
|
{
|
|
return crossSectionHandlers;
|
|
}
|
|
|
|
void ClearStages(void);
|
|
|
|
void AddStage(G4ScreenedCollisionStage* stage) { collisionStages.push_back(stage); }
|
|
|
|
G4CoulombKinematicsInfo& GetKinematics() { return kinematics; }
|
|
|
|
void SetValidCollision(G4bool flag) { validCollision = flag; }
|
|
|
|
G4bool GetValidCollision() const { return validCollision; }
|
|
|
|
/// for internal use, primarily.
|
|
class G4ParticleChange& GetParticleChange()
|
|
{
|
|
return static_cast<G4ParticleChange&>(*pParticleChange);
|
|
}
|
|
|
|
/// to partition it into NIEL and ionizing energy.
|
|
|
|
void DepositEnergy(G4int z1, G4double a1, const G4Material* material, G4double energy);
|
|
|
|
protected:
|
|
G4double highEnergyLimit;
|
|
|
|
G4double lowEnergyLimit;
|
|
|
|
/// to cross over to G4MSC
|
|
G4double processMaxEnergy;
|
|
G4String screeningKey;
|
|
G4bool generateRecoils, avoidReactions;
|
|
G4double recoilCutoff, physicsCutoff;
|
|
G4bool registerDepositedEnergy;
|
|
G4double IonizingLoss, NIEL;
|
|
G4double MFPScale;
|
|
G4double hardeningFraction, hardeningFactor;
|
|
|
|
G4ScreenedCoulombCrossSection* externalCrossSectionConstructor;
|
|
std::vector<G4ScreenedCollisionStage*> collisionStages;
|
|
|
|
std::map<G4int, G4ScreenedCoulombCrossSection*> crossSectionHandlers;
|
|
|
|
G4bool validCollision;
|
|
G4CoulombKinematicsInfo kinematics;
|
|
const G4VNIELPartition* NIELPartitionFunction;
|
|
};
|
|
|
|
// A customized G4CrossSectionHandler which gets its data from
|
|
// an external program
|
|
|
|
class G4NativeScreenedCoulombCrossSection : public G4ScreenedCoulombCrossSection
|
|
{
|
|
public:
|
|
G4NativeScreenedCoulombCrossSection();
|
|
|
|
G4NativeScreenedCoulombCrossSection(const G4NativeScreenedCoulombCrossSection& src)
|
|
: G4ScreenedCoulombCrossSection(src), phiMap(src.phiMap)
|
|
{}
|
|
|
|
G4NativeScreenedCoulombCrossSection(const G4ScreenedCoulombCrossSection& src)
|
|
: G4ScreenedCoulombCrossSection(src)
|
|
{}
|
|
|
|
virtual ~G4NativeScreenedCoulombCrossSection();
|
|
|
|
virtual void LoadData(G4String screeningKey, G4int z1, G4double m1, G4double recoilCutoff);
|
|
|
|
virtual G4ScreenedCoulombCrossSection* create()
|
|
{
|
|
return new G4NativeScreenedCoulombCrossSection(*this);
|
|
}
|
|
|
|
// get a list of available keys
|
|
std::vector<G4String> GetScreeningKeys() const;
|
|
|
|
typedef G4_c2_function& (*ScreeningFunc)(G4int z1, G4int z2, size_t nPoints, G4double rMax,
|
|
G4double* au);
|
|
|
|
void AddScreeningFunction(G4String name, ScreeningFunc fn) { phiMap[name] = fn; }
|
|
|
|
private:
|
|
// this is a map used to look up screening function generators
|
|
std::map<std::string, ScreeningFunc> phiMap;
|
|
};
|
|
|
|
G4_c2_function& ZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
|
|
|
|
G4_c2_function& MoliereScreening(G4int z1, G4int z2, size_t npoints, G4double rMax,
|
|
G4double* auval);
|
|
|
|
G4_c2_function& LJScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
|
|
|
|
G4_c2_function& LJZBLScreening(G4int z1, G4int z2, size_t npoints, G4double rMax, G4double* auval);
|
|
|
|
#endif
|