Physics Processes¶
Overview¶
Physics processes describe how particles interact with a material. Seven major categories of processes are provided by Geant4:
electromagnetic,
hadronic,
decay,
photolepton-hadron,
optical,
parameterization, and
transportation.
The generalization and abstraction of physics processes is a key issue in the design of Geant4. All physics processes are treated in the same manner from the tracking point of view. The Geant4 approach enables anyone to create a process and assign it to a particle type. This openness should allow the creation of processes for novel, domain-specific or customised purposes by individuals or groups of users.
Each process has two groups of methods which play an important role in
tracking, GetPhysicalInteractionLength
(GPIL) and DoIt
. The GPIL
method gives the step length from the current space-time point to the
next space-time point. It does this by calculating the probability of
interaction based on the process’s cross section information. At the end
of this step the DoIt
method should be invoked. The DoIt
method
implements the details of the interaction, changing the particle’s
energy, momentum, direction and position, and producing secondary tracks
if required. These changes are recorded as G4VParticleChange
objects (see Particle change).
G4VProcess¶
G4VProcess
is the base class for all physics processes. Each physics
process must implement virtual methods of G4VProcess
which describe
the interaction (DoIt) and determine when an interaction should occur
(GPIL). In order to accommodate various types of interactions
G4VProcess
provides three DoIt
methods:
G4VParticleChange* AlongStepDoIt( const G4Track& track, const G4Step& stepData )
This method is invoked while
G4SteppingManager
is transporting a particle through one step. The correspondingAlongStepDoIt
for each defined process is applied for every step regardless of which process produces the minimum step length. Each resulting change to the track information is recorded and accumulated inG4Step
. After all processes have been invoked, changes due toAlongStepDoIt
are applied toG4Track
, including the particle relocation and the safety update. Note that after the invocation ofAlongStepDoIt
, the endpoint of theG4Track
object is in a new volume if the step was limited by a geometric boundary. In order to obtain information about the old volume,G4Step
must be accessed, since it contains information about both endpoints of a step.G4VParticleChange* PostStepDoIt( const G4Track& track, const G4Step& stepData )
This method is invoked at the end point of a step, only if its process has produced the minimum step length, or if the process is forced to occur.
G4Track
will be updated after each invocation ofPostStepDoIt
, in contrast to theAlongStepDoIt
method.G4VParticleChange* AtRestDoIt( const G4Track& track, const G4Step& stepData )
This method is invoked only for stopped particles, and only if its process produced the minimum step length or the process is forced to occur.
For each of the above DoIt
methods G4VProcess
provides a
corresponding pure virtual GPIL method:
PostStepGetPhysicalInteractionLength
G4double PostStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4ForceCondition* condition )
This method generates the step length allowed by its process. It also provides a flag to force the interaction to occur regardless of its step length.
AlongStepGetPhysicalInteractionLength
G4double AlongStepGetPhysicalInteractionLength(const G4Track& track, G4double previousStepSize, G4double currentMinimumStep, G4double& proposedSafety, G4GPILSelection* selection )
This method generates the step length allowed by its process.
AtRestGetPhysicalInteractionLength
G4double AtRestGetPhysicalInteractionLength(const G4Track& track, G4ForceCondition* condition )
This method generates the step length in time allowed by its process. It also provides a flag to force the interaction to occur regardless of its step length.
Other pure virtual methods in G4VProcess
follow:
virtual G4bool IsApplicable(const G4ParticleDefinition&)
returns true if this process object is applicable to the particle type.
virtual void PreparePhysicsTable(const G4ParticleDefinition&)
andvirtual void BuildPhysicsTable(const G4ParticleDefinition&)
is messaged by the process manager, whenever cross section tables should be prepared and rebuilt due to changing cut-off values. It is not mandatory if the process is not affected by cut-off values.
virtual void StartTracking()
andvirtual void EndTracking()
are messaged by the tracking manager at the beginning and end of tracking the current track.
Other base classes for processes¶
Specialized processes may be derived from seven additional virtual base
classes which are themselves derived from G4VProcess
. Three of these
classes are used for simple processes:
G4VRestProcess
Processes using only the
AtRestDoIt
method.example: neutron capture
G4VDiscreteProcess
Processes using only the
PostStepDoIt
method.example: Compton scattering, hadron inelastic interaction
The other four classes are provided for rather complex processes:
G4VContinuousDiscreteProcess
Processes using both
AlongStepDoIt
andPostStepDoIt
methods.example: transportation, ionisation(energy loss and delta ray)
G4VRestDiscreteProcess
Processes using both
AtRestDoIt
andPostStepDoIt
methods.example: positron annihilation, decay (both in flight and at rest)
G4VRestContinuousProcess
Processes using both
AtRestDoIt
andAlongStepDoIt
methods.G4VRestContinuousDiscreteProcess
Processes using
AtRestDoIt
,AlongStepDoIt and
PostStepDoIt methods.
Particle change¶
G4VParticleChange
and its descendants are used to store the final
state information of the track, including secondary tracks, which has
been generated by the DoIt
methods. The instance of
G4VParticleChange
is the only object whose information is updated by
the physics processes, hence it is responsible for updating the step.
The stepping manager collects secondary tracks and only sends requests
via particle change to update G4Step
.
G4VParticleChange
is introduced as an abstract class. It has a
minimal set of methods for updating G4Step
and handling secondaries.
A physics process can therefore define its own particle change derived
from G4VParticleChange
. Three pure virtual methods are provided,
virtual G4Step* UpdateStepForAtRest( G4Step* step)
,virtual G4Step* UpdateStepForAlongStep( G4Step* step )
, andvirtual G4Step* UpdateStepForPostStep( G4Step* step)
,
which correspond to the three DoIt
methods of G4VProcess
. Each
derived class should implement these methods.
Electromagnetic Interactions¶
This section summarizes the electromagnetic (EM) physics processes which are provided with Geant4. Extended information are available at EM web pages. For details on the implementation of these processes please refer to the Physics Reference Manual.
To use the electromagnetic physics data files are needed. The user should set the environment variable G4LEDATA to the directory with this files. These files are distributed together with Geant4 and can be obtained via the Geant4 download web page.
Electromagnetic Processes¶
The following is a summary of the electromagnetic processes available in Geant4.
Photon processes
Gamma conversion (also called pair production, class name
G4GammaConversion
)Photo-electric effect (class name
G4PhotoElectricEffect
)Compton scattering (class name
G4ComptonScattering
)Rayleigh scattering (class name
G4RayleighScattering
)Muon pair production (class name
G4GammaConversionToMuons
)
Electron/positron processes
Ionisation and delta ray production (class name
G4eIonisation
)Bremsstrahlung (class name
G4eBremsstrahlung
)e+e- pair production (class name
G4ePairProduction
)Multiple scattering (class name
G4eMultipleScattering
)Positron annihilation into two gammas (class name
G4eplusAnnihilation
)Positron annihilation into two muons (class name
G4AnnihiToMuPair
)Positron annihilation into hadrons (class name
G4eeToHadrons
)
Muon processes
Ionisation and delta ray production (class name
G4MuIonisation
)Bremsstrahlung (class name
G4MuBremsstrahlung
)e+e- pair production (class name
G4MuPairProduction
)Multiple scattering (class name
G4MuMultipleScattering
)
Hadron/ion processes
Ionisation (class name
G4hIonisation
)Ionisation for ions (class name
G4ionIonisation
)Ionisation for heavy exotic particles (class name
G4hhIonisation
)Ionisation for classical magnetic monopole (class name
G4mplIonisation
)Multiple scattering (class name
G4hMultipleScattering
)Bremsstrahlung (class name
G4hBremsstrahlung
)e+e- pair production (class name
G4hPairProduction
)
Coulomb scattering processes
Alternative process for simulation of single Coulomb scattering of all charged particles (class name
G4CoulombScattering
)Alternative process for simulation of single Coulomb scattering of ions (class name
G4ScreenedNuclearRecoil
)
Processes for simulation of polarized electron and gamma beams
Compton scattering of circularly polarized gamma beam on polarized target (class name
G4PolarizedCompton
)Pair production induced by circularly polarized gamma beam (class name
G4PolarizedGammaConversion
)Photo-electric effect induced by circularly polarized gamma beam (class name
G4PolarizedPhotoElectricEffect
)Bremsstrahlung of polarized electrons and positrons (class name
G4ePolarizedBremsstrahlung
)Ionisation of polarized electron and positron beam (class name
G4ePolarizedIonisation
)Annihilation of polarized positrons (class name
G4eplusPolarizedAnnihilation
)
Processes for simulation of X-rays and optical protons production by charged particles
Synchrotron radiation (class name
G4SynchrotronRadiation
)Transition radiation (class name
G4TransitionRadiation
)Cerenkov radiation (class name
G4Cerenkov
)Scintillations (class name
G4Scintillation
)
The processes described above use physics model classes, which may be combined according to particle energy. It is possible to change the energy range over which different models are valid, and to apply other models specific to particle type, energy range, and G4Region. The following alternative models are available in the standard EM sub-library:
Ionisation in thin absorbers (class name
G4PAIModel
)Ionisation in thin absorbers (class name
G4PAIPhotModel
)Ionisation in low-density media (class name
G4BraggIonGasModel
)Ionisation in low-density media (class name
G4BetheBlochIonGasModel
)Ionisation of relativistic ions (class name
G4AtimaEnergyLossModel
)Ionisation of relativistic ions (class name
G4LindhardSorensenIonModel
)Electron/positron bremsstrahlung (class name
G4BetheHeitler5DModel
)Positron annihilation into 2 or 3 gamma (class name
G4eplusTo2GammaOKVIModel
)Multiple scattering (class name
G4UrbanMscModel
)Multiple scattering (class name
G4GoudsmitSaundersonMscModel
)Multiple scattering (class name
G4WentzelVIModel
)Multiple scattering (class name
G4LowEWentzelVIModel
)Single scattering (class name
G4eCoulombScatteringModel
)Single scattering (class name
G4eSingleCoulombScatteringModel
)Single scattering (class name
G4IonCoulombScatteringModel
)
It is recommended to use physics constructor classes provided with
reference physics lists (in subdirectory
source/physics_lists/constructors/electromagnetic
of the Geant4
source distribution):
default EM physics, multiple scattering is simulated with “UseSafety” type of step limitation by combined
G4WentzelVIModel
andG4eCoulombScatteringModel
for all particle types, for of e+- below 100 MeVG4UrbanMscModel
is used,G4LivermorePhotoElectricModel
is used for simulation of the photo-electric effect, the Rayleigh scattering process is enabled below 1 MeV, physics tables are built from 100 eV to 100 TeV, 7 bins per energy decade of physics tables are used (class nameG4EmStandardPhysics
),optional EM physics providing fast but less accurate electron transport due to “Minimal” method of step limitation by multiple scattering,
RangeFactor = 0.2
, Rayleigh scattering is disabled, photo-electric effect is usingG4PEEffectFluoModel
, and enabled “ApplyCuts” option (class nameG4EmStandardPhysics_option1
),optional EM physics providing fast but less accurate electron transport due to “Minimal” method of step limitation by multiple scattering,
RangeFactor = 0.2
, “Simple” method of step limitation by multiple scattering, Rayleigh scattering is disabled, and photo-electric effect is usingG4PEEffectFluoModel
(class nameG4EmStandardPhysics_option2
)EM physics for simulation with high accuracy due to “UseDistanceToBoundary” multiple scattering step limitation and usage of
G4UrbanMscModel
for all charged particles, reduced finalRange parameter of stepping function optimized per particle type, alternative modelG4KleinNishinaModel
for Compton scattering, enabled fluorescence, enabled nuclear stopping,G4Generator2BS
angular generator for bremsstrahlung,G4IonParameterisedLossModel
for ion ionisation,G4ePairProduction
for electron/positron, 20 bins energy decade of physics tables, and 10 eV low-energy limit for tables (class nameG4EmStandardPhysics_option3
)Combination of EM models for simulation with high accuracy includes multiple scattering with “UseSafetyPlus” type of step limitation by combined
G4WentzelVIModel
andG4eCoulombScatteringModel
for all particle types, for of e+- below 100 MeVG4GoudsmitSaundersonMscModel
is used,RangeFactor = 0.08
,Scin = 3
(error free stepping near geometry boundaries), reduced finalRange parameter of stepping function optimized per particle type, enabled fluorescence, enabled nuclear stopping, enable accurate angular generator for ionisation models,G4LowEPComptonModel
below 20 MeV andG4KleinNishinaModel
above,G4BetheHeitler5DModel
for gamma conversion,G4LivermoreIonisationModel
for electrons and positrons below 100 keV,G4IonParameterisedLossModel
for ion ionisation,G4Generator2BS
angular generator for bremsstrahlung,G4ePairProduction
for electron/positron, and 20 bins per energy decade of physics tables, (class nameG4EmStandardPhysics_option4
)Models based on Livermore data bases for electrons and gamma, enabled Rayleigh scattering, enabled fluorescence,
G4BetheHeitler5DModel
is used for gamma conversion, enabled nuclear stopping, enable accurate angular generator for ionisation models,G4IonParameterisedLossModel
for ion ionisation, for of e+- below 100 MeVG4GoudsmitSaundersonMscModel
is used with “UseSafetyPlus” multiple scattering step limitation,RangeFactor = 0.08
,Scin = 3
(error free stepping near geometry boundaries),G4Generator2BS
angular generator for bremsstrahlung,G4ePairProduction
for electron/positron, and 20 bins per energy decade of physics tables, (G4EmLivermorePhysics
);Models for simulation of linear polarized gamma based on Livermore data bases for electrons and gamma (
G4EmLivermorePolarizedPhysics
);Models based on Livermore data bases and new model for Compton scattering
G4LowEPComptonModel
,G4BetheHeitler5DModel
is used for gamma conversion, low-energy model of multiple scatteringG4LowEWenzelMscModel
, and newG4LindhardSorensenIonModel
for ions above 10 MeV,G4hBremsstrahlung
andG4hPairProduction
for ions (G4EmLowEPPhysics
);Penelope2008 models for electrons, positrons and gamma, enabled Rayleigh scattering, enabled fluorescence, enabled nuclear stopping, enable accurate angular generator for ionisation models,
G4IonParameterisedLossModel
for ion ionisation, for of e+- below 100 MeVG4GoudsmitSaundersonMscModel
is used with “UseSafetyPlus” multiple scattering step limitation,RangeFactor = 0.08
,Scin = 3
(error free stepping near geometry boundaries), and 20 bins per energy decade of physics tables, (G4EmPenelopePhysics
);Experimental physics with multiple scattering of e+- below 100 MeV simulated by
G4GoudsmitSaundersonMscModel
is done on top of the default EM physics (G4EmStandardPhysicsGS
);Experimental physics is done on top of the default EM physics with multiple scattering of e+- below 100 MeV simulated by a combination of
G4WentzelVIModel
andG4eCoulombScatteringModel
, and for ions above 2 MeVG4AtimaEnergyLossModel
(G4EmStandardPhysicsWVI
);Experimental physics with single scattering models instead of multiple scattering is done on top of the default EM physics, for all leptons and hadrons
G4eCoulombScatteringModel
is used, for ions -G4IonCoulombScatteringModel
(G4EmStandardPhysicsSS
);Low-energy Geant4-DNA physics (
G4EmDNAPhysics
).Alternative low-energy Geant4-DNA physics constructors (
G4EmDNAPhysics_optionX
, where X is 1 to 8). Refer to Geant4-DNA section.
Examples of the registration of these physics constructor and
construction of alternative combinations of options are shown in basic,
extended and advanced examples, which can be found in the subdirectories
examples/basic
, examples/extended/electromagnetic
,
examples/medical
, examples/advanced
, and of the Geant4 source
distribution. Examples illustrating the use of electromagnetic processes
are available as part of the Geant4
release.
Options are available for steering of electromagnetic processes.
These options may be invoked either by UI commands or by the new C++
interface class G4EmParameters
. The interface
G4EmParameters::Instance()
is thread safe, EM parameters are shared
between threads, and parameters are shared between all EM processes.
Parameters may be modified at G4State_PreInit or G4State_Idle states
of Geant4. Note, that when any of EM physics constructor is instantiated
a default set of EM parameters for this EM physics configuration is
defined. So, parameters modification should be applied only after. This
class has the following public methods:
Dump()
StreamInfo(std::ostream&)
SetDefaults()
SetLossFluctuations(G4bool)
SetBuildCSDARange(G4bool)
SetLPM(G4bool)
SetSpline(G4bool)
SetUseCutAsFinalRange(G4bool)
SetApplyCuts(G4bool)
SetFluo(G4bool val)
SetBeardenFluoDir(G4bool val)
SetAuger(G4bool val)
SetAugerCascade(G4bool val) - obsolete
SetPixe(G4bool val)
SetDeexcitationIgnoreCut(G4bool val)
SetLateralDisplacement(G4bool val)
SetLateralDisplacementAlg96(G4bool val)
SetMuHadLateralDisplacement(G4bool val)
SetLatDisplacementBeyondSafety(G4bool val)
ActivateAngularGeneratorForIonisation(G4bool val)
SetUseMottCorrection(G4bool val)
SetIntegral(G4bool val)
SetBirksActive(G4bool val)
SetUseICRU90Data(G4bool val)
SetDNAFast(G4bool val)
SetDNAStationary(G4bool val)
SetDNAElectronMsc(G4bool val)
SetGeneralProcessActive(G4bool val)
SetEnableSamplingTable(G4bool val)
SetEnablePolarisation(G4bool val)
SetDirectionalSplitting(G4bool val)
SetQuantumEntanglement(G4bool val)
SetEmSaturation(G4EmSaturation*)
SetOnIsolated(G4bool val)
ActivateDNA(G4bool val)
SetMinSubRange(G4double)
SetMinEnergy(G4double)
SetMaxEnergy(G4double)
SetMaxEnergyForCSDARange(G4double)
SetLowestElectronEnergy(G4double)
SetLowestMuHadEnergy(G4double)
SetLowestTripletEnergy(G4double)
SetLinearLossLimit(G4double)
SetBremsstrahlungTh(G4double)
SetLambdaFactor(G4double)
SetFactorForAngleLimit(G4double)
SetMscThetaLimit(G4double)
SetMscEnergyLimit(G4double)
SetMscRangeFactor(G4double)
SetMscMuHadRangeFactor(G4double)
SetMscGeomFactor(G4double)
SetMscSafetyFactor(G4double)
SetMscLambdaLimit(G4double)
SetMscSkin(G4double)
SetScreeningFactor(G4double)
SetMaxNIELEnergy(G4double)
SetMaxEnergyFor5DMuPair(G4double)
SetStepFunction(G4double, G4double)
SetStepFunctionMuHad(G4double, G4double)
SetDirectionalSplittingRadius(G4double)
SetDirectionalSplittingTarget(const G4ThreeVector&)
SetNumberOfBins(G4int)
SetNumberOfBinsPerDecade(G4int)
SetVerbose(G4int)
SetWorkerVerbose(G4int)
SetMscStepLimitType(G4MscStepLimitType val)
SetMscMuHadStepLimitType(G4MscStepLimitType val)
SetNuclearFormFactorType(G4NuclearFormFactorType val)
SetDNAeSolvationSubType(G4DNAModelSubType val)
SetConversionType(G4int val)
SetPIXECrossSectionModel(const G4String&)
SetPIXEElectronCrossSectionModel(const G4String&)
AddPAIModel(const G4String& particle, const G4String& region, const G4String& type)
AddMicroElecModel(const G4String& region)
AddDNA(const G4String& region, const G4String& type)
AddPhysics(const G4String& region, const G4String& physics_type)
SetSubCutoff(G4bool, const G4String& region)
SetDeexActiveRegion(const G4String& region, G4bool, G4bool, G4bool)
SetProcessBiasingFactor(const G4String& process, G4double, G4bool)
ActivateForcedInteraction(const G4String& process, const G4String& region, G4double, G4bool)
ActivateSecondaryBiasing(const G4String& process, const G4String& region, G4double, G4double)
SetDirectionalSplitting(G4int v)
SetDirectionalSplittingTarget(const G4ThreeVector& v)
SetDirectionalSplittingRadius(G4double)
The old interface class G4EmProcessOptions
is still available but
but is strongly recommended not to be used. It will be removed in the
next major release.
The corresponding UI command can be accessed in the UI subdirectories “/process/eLoss”, “/process/em”, and “/process/msc”. The following types of step limitation by multiple scattering are available:
fMinimal - simplified step limitation (used in _EMV and _EMX Physics Lists)
fUseSafety - default
fUseDistanceToBoundary - advance method of step limitation used in EM examples, required parameter skin > 0 , should be used for setup without magnetic field
fUseSafetyPlus - advanced method may be used with magnetic field
G4EmCalculator is a class which provides access to cross sections
and stopping powers. This class can be used anywhere in the user code
provided the physics list has already been initialised (G4State_Idle).
G4EmCalculator
has “Get” methods which can be applied to materials for
which physics tables are already built, and “Compute” methods which can
be applied to any material defined in the application or existing in the
Geant4 internal database. The public methods of this class are:
GetDEDX(kinEnergy,particle,material,const G4Region* r=nullptr)
GetRangeFromRestrictedDEDX(kinEnergy,particle,material,const G4Region* r=nullptr)
GetCSDARange(kinEnergy,particle,material,const G4Region* r=nullptr)
GetRange(kinEnergy,particle,material,const G4Region* r=nullptr)
GetKinEnergy(range,particle,material,const G4Region* r=nullptr)
GetCrossSectionPerVolume(kinEnergy,particle,material,const G4Region* r=nullptr)
GetShellIonisationCrossSectionPerAtom(particle,Z,shell,kinEnergy)
GetMeanFreePath(kinEnergy,particle,material,const G4Region* r=nullptr)
PrintDEDXTable(particle)
PrintRangeTable(particle)
PrintInverseRangeTable(particle)
ComputeDEDX(kinEnergy,particle,process,material,cut=DBL_MAX)
ComputeElectronicDEDX(kinEnergy,particle,material,cut=DBL_MAX)
ComputeDEDXForCutInRange(kinEnergy,particle,material,cut=DBL_MAX)
ComputeNuclearDEDX(kinEnergy,particle,material,cut=DBL_MAX)
ComputeTotalDEDX(kinEnergy,particle,material,cut=DBL_MAX)
ComputeCrossSectionPerVolume(kinEnergy,particle,process,material,cut=nullptr)
ComputeCrossSectionPerAtom(kinEnergy,particle,process,Z,A,cut=nullptr)
ComputeCrossSectionPerShell(kinEnergy,particle,process,Z,shellIdx,cut=nullptr)
ComputeGammaAttenuationLength(kinEnergy,material)
ComputeShellIonisationCrossSectionPerAtom(particle,Z,shell,kinEnergy)
ComputeMeanFreePath(kinEnergy,particle,process,material,cut=nullptr)
ComputeEnergyCutFromRangeCut(range,particle,material)
FindParticle(const G4String&)
FindIon(G4int Z, G4int A)
FindMaterial(const G4String&)
FindRegion(const G4String&)
FindCouple(const G4Material*, const const G4Region* r=nullptr)
FindProcess(particle, const G4String& processName)
SetVerbose(G4int)
For these interfaces, particles, materials, or processes may be pointers
(const G4ParticleDefinition*
, const G4Material*
)
or strings with names (const G4String&
).
G4NIELCalculator is a class which provides computation of NIEL energy loss
at a step independently on cuts and tracking.
G4NIELCalculator
has follow methods:
ComputeNIEL(const G4Step*)
RecoilEnergy(const G4Step*)
AddEmModel(G4VEmModel*)
the last method allows customisation of NIEL model.
Low Energy Electromagnetic Library¶
A physical interaction is described by a process class which can handle physics models, described by model classes. The following is a summary of the Low Energy Electromagnetic physics models available in Geant4. Further information is available in the web pages of the Geant4 Low Energy Electromagnetic Physics Working Group, accessible from the Geant4 web site, “who we are” section, then “working groups”.
The physics content of these models is documented in the Geant4 Physics Reference Manual. They are based on the Livermore data library, on the ICRU73 data tables or on the Penelope2008 Monte Carlo code. They adopt the same software design as the “standard” Geant4 electromagnetic models.
Examples of the registration of physics constructor with low-energy
electromagnetic models are shown in Geant4 extended examples
(examples/extended/electromagnetic
and examples/extended/medical
in the Geant4 source distribution). Advanced examples
(examples/advanced
in the Geant4
source distribution) illustrate alternative instantiation of these
processes. Both are available as part of the Geant4 release.
Production Cuts¶
Remember that production cuts for secondaries can be specified as range cuts, which are converted at initialisation time into energy thresholds for secondary gamma, electron, positron and proton production. The cut for proton is applied by elastic scattering processes to all recoil ions.
A range cut value is set by default to 0.7 mm in Geant4 reference physics lists. This value can be specified in the optional SetCuts() method of the user Physics list or via UI commands. For e.g. to set a range cut of 10 micrometers, one can use
/run/setCut 0.01 mm
or, for a given particle type (for e.g. electron),
/run/setCutForAGivenParticle e- 0.01 mm
If a range cut equivalent to an energy lower than 990 eV is specified, the energy cut is still set to 990 eV. In order to decrease this value (for e.g. down to 250 eV, in order to simulate low energy emission lines of the fluorescence spectrum), one may use the following UI command before the “/run/initialize” command
/cuts/setLowEdge 250 eV
or alternatively directly in the user Physics list, in the optional SetCuts() method, using:
G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(250*eV, 1*GeV);
A command is also available in order to disable usage of production threshold for fluorescence and Auger electron production
/process/em/deexcitationIgnoreCut true
Angular Generators¶
For part of EM processes it is possible to factorise out sampling of
secondary energy and direction. Using an interface G4VEmModel
base
class SetAngularDistribution(G4VEmAngularDistribution*)
it is
possible to substitute default angular generator of a model. Angular
generators in standard and lowenergy sub-packages follow the same
abstract interface.
For photoelectric models several angular generators are available:
G4SauterGavrilaAngularDistribution (default);
G4PhotoElectricAngularGeneratorSauterGavrila;
G4PhotoElectricAngularGeneratorPolarized.
For bremsstrahlung models following angular generators are available:
G4ModifiedTsai (default);
G4DipBustGenerator;
G4Generator2BS;
G4Generator2BN;
G4PenelopeBremsstrahlungAngular.
For gamma conversion models following angular generators are available:
G4ModifiedTsai (default);
G4DipBustGenerator.
For models of ionisation a new optional angular generator is available:
G4DeltaAngle.
Electromagnetics secondary biasing¶
It may be useful to create more than one secondary at an interaction. For example, electrons incident on a target in a medical linac produce photons through bremsstrahlung. The variance reduction technique of bremsstrahlung splitting involves choosing N photons from the expected distribution, and assigning each a weight of 1/N.
Similarly, if the secondaries are not important, one can kill them with a survival probability of 1/N. The weight of the survivors is increased by a factor N. This is known as Russian roulette.
Neither biasing technique is applied if the resulting daughter particles would have a weight below 1/N, in the case of brem splitting, or above 1, in the case of Russian roulette.
These techniques can be enabled in Geant4 electromagnetics with the macro commands
/process/em/setSecBiasing processName Region factor energyLimit energyUnit
where: processName is the name of the process to apply the biasing to; Region is the region in which to apply biasing; factor is the inverse of the brem splitting or Russian roulette factor (1/N); energyLimit energyUnit is the high energy limit. If the first secondary has energy above this limit, no biasing is applied.
For example,
/process/em/setSecBiasing eBrem target 10 100 MeV
will result in electrons undergoing bremsstrahlung in the target region being split 10 times (if the first photon sampled has an energy less than 100 MeV).
Note that the biasing needs to be specified for each process individually. To apply Russian Roulette to daughter electrons from interactions of photons, issue the macro command for the processes phot, compt, conv.
Directional splitting¶
This biasing may be enabled based on the direction of the outgoing particles (“directional splitting”). The user may specify a spherical volume of interest by giving the center and radius of the volume. In an interaction, if the incident particle has high weight, the outgoing particles are split. N particles from the distribution are created, each with weight 1/N. For each particle, if it is not directed towards the volume of interest, Russian Roulette is played. Typically one will want directional splitting to take place for all interactions.
For example,
/process/em/setDirectionalSplitting true
/process/em/setDirectionalSplittingTarget 1000 0 0 mm # x, y, z components of center
/process/em/setDirectionalSplittingRadius 10 cm
/process/em/setSecBiasing eBrem world 100 100 MeV
/process/em/setSecBiasing Rayl world 100 100 MeV
/process/em/setSecBiasing phot world 100 100 MeV
/process/em/setSecBiasing compt world 100 100 MeV
/process/em/setSecBiasing annihil world 100 100 MeV
Reference: BEAMnrc Users Manual, D.W.O Rogers, B. Walters, I. Kawrakow. NRCC Report PIRS-0509(A)revL, available at http://www.irs.inms.nrc.ca/inms/irs/BEAM/beamhome.html
Livermore Data Based Models¶
Photon models
Photo-electric effect (class
G4LivermorePhotoElectricModel
)Polarized Photo-electric effect (class
G4LivermorePolarizedPhotoElectricModel
)Compton scattering (class
G4LivermoreComptonModel
)Compton scattering (class
G4LowEPComptonModel
)Polarized Compton scattering (class
G4LivermorePolarizedComptonModel
)Rayleigh scattering (class
G4LivermoreRayleighModel
)Polarized Rayleigh scattering (class
G4LivermorePolarizedRayleighModel
)Gamma conversion (also called pair production, class
G4LivermoreGammaConversionModel
)Nuclear gamma conversion (class
G4LivermoreNuclearGammaConversionModel
)Radiative correction for pair production (class
G4LivermoreGammaConversionModelRC
)Polarized gamma conversion (class
G4LivermorePolarizedGammaConversionModel
)
Electron models
Bremsstrahlung (class
G4LivermoreBremsstrahlungModel
)Ionisation and delta ray production (class
G4LivermoreIonisationModel
)
ICRU73 Based Ion Model¶
Ionisation and delta ray production (class G4IonParametrisedLossModel)
The ion model uses data files initially converted from the ICRU 73 report. Later authors of the ICRU 73 report provided Geant4 recomputed tables for more combinations of projectile and target ions. In 2015 newer calculation results were provided. The algorithm of selection of ion stopping powers applying following condition: if a projectile/target combination exists in the data base and the projectile energy is below 1 GeV/nucleon then tabulated data is used, otherwise applies a Bethe-Bloch based formalism. For compounds, ICRU 73 stopping powers are employed if the material name coincides with the name of Geant4 NIST materials (e.g. G4_WATER). Elemental materials are matched to the corresponding ICRU 73 stopping powers by means of the atomic number of the material. The material name may be arbitrary in this case. For a list of applicable materials, the user is referred to the ICRU 73 report.
The model requires data files to be copied by the user to his/her code repository. These files are distributed together with the Geant4 release. The user should set the environment variable G4LEDATA to the directory where he/she has copied the files.
The model is dedicated to be used with the G4ionIonisation process and its applicability is restricted to G4GenericIon particles. The ion model is not used by default by this process and must be instantiated and registered by the user:
G4ionIonisation* ionIoni = new G4ionIonisation();
ionIoni -> SetEmModel(new G4IonParametrisedLossModel());
Penelope2008 Based Models¶
Photon models
Compton scattering (class
G4PenelopeComptonModel
)Rayleigh scattering (class
G4PenelopeRayleighModel
)Gamma conversion (also called pair production, class
G4PenelopeGammaConversionModel
)Photo-electric effect (class
G4PenelopePhotoElectricModel
)
Electron models
Bremsstrahlung (class
G4PenelopeBremsstrahlungModel
)Ionisation and delta ray production (class
G4PenelopeIonisationModel
)
Positron models
Bremsstrahlung (class
G4PenelopeBremsstrahlungModel
)Ionisation and delta ray production (class
G4PenelopeIonisationModel
)Positron annihilation (class
G4PenelopeAnnihilationModel
)
All Penelope models can be applied up to a maximum energy of 100 GeV, although it is advisable not to use them above a few hundreds of MeV.
Options are available in the all Penelope Models, allowing to set (and retrieve) the verbosity level of the model, namely the amount of information which is printed on the screen.
SetVerbosityLevel(G4int)
GetVerbosityLevel()
The default verbosity level is 0 (namely, no textual output on the screen). The default value should be used in general for normal runs. Higher verbosity levels are suggested only for testing and debugging purposes.
The verbosity scale defined for all Penelope processes is the following:
0 = no printout on the screen (default)
1 = issue warnings only in the case of energy non-conservation in the final state (should never happen)
2 = reports full details on the energy budget in the final state
3 = writes also information on cross section calculation, data file opening and sampling of atoms
4 = issues messages when entering in methods
Very Low energy Electromagnetic Processes (Geant4-DNA extension)¶
The Geant4 low energy electromagnetic Physics package has been extended down to energies of a few electron Volts suitable for the simulation of radiation effects in liquid water for applications in micro/nanodosimetry at the cellular and sub-cellular level. These developments take place in the framework of the on-going Geant4-DNA project (see more in the Geant4-DNA web pages or in the EM web pages of the Geant4 Electromagnetic Physics Working Group).
The Geant4 -DNA process and model classes apply to electrons, protons, hydrogen, alpha particles and their charge states, in liquid water (“G4_WATER” material).
Electron processes and models
Elastic scattering:
process class is G4DNAElastic
four alternative model classes are: G4DNAScreenedRutherfordElasticModel or G4DNAChampionElasticModel (default) or G4DNAUeharaScreenedRutherfordElasticModel or G4DNACPA100ElasticModel
Excitation
process class is G4DNAExcitation
model class is G4DNABornExcitationModel (default) or G4DNAEmfietzoglouExcitationModel or G4DNACPA100ExcitationModel
Ionisation
process class is G4DNAIonisation
model class is G4DNABornIonisationModel (default) or G4DNAEmfietzoglouIonisationModel or G4DNACPA100IonisationModel
Attachment
process class is G4DNAAttachment
model class is G4DNAMeltonAttachmentModel
Vibrational excitation
process class is G4DNAVibExcitation
model class is G4DNASancheExcitationModel
Proton processes and models
Elastic scattering:
process class is G4DNAElastic
G4DNAIonElasticModel
Excitation
process class is G4DNAExcitation
two complementary model classes are G4DNAMillerGreenExcitationModel (below 500 keV) and G4DNABornExcitationModel (above)
Ionisation
process class is G4DNAIonisation
two complementary model classes are G4DNARuddIonisationModel (below 500 keV) and G4DNABornIonisationModel (above)
Charge decrease
process class is G4DNAChargeDecrease
model class is G4DNADingfelderChargeDecreaseModel
Hydrogen processes and models
Elastic scattering :
process class is G4DNAElastic
G4DNAIonElasticModel
Excitation
process class is G4DNAExcitation
model class is G4DNAMillerGreenExcitationModel
Ionisation
process class is G4DNAIonisation
model class is G4DNARuddIonisationModel
Charge increase
process class is G4DNAChargeIncrease
model class is G4DNADingfelderChargeIncreaseModel
Helium (neutral) processes and models
Elastic scattering :
process class is G4DNAElastic
G4DNAIonElasticModel
Excitation
process class is G4DNAExcitation
model class is G4DNAMillerGreenExcitationModel
Ionisation
process class is G4DNAIonisation
model class is G4DNARuddIonisationModel
Charge increase
process class is G4DNAChargeIncrease
model class is G4DNADingfelderChargeIncreaseModel
Helium+ (ionized once) processes and models
Elastic scattering :
process class is G4DNAElastic
G4DNAIonElasticModel
Excitation
process class is G4DNAExcitation
model class is G4DNAMillerGreenExcitationModel
Ionisation
process class is G4DNAIonisation
model classes is G4DNARuddIonisationModel
Charge increase
process class is G4DNAChargeIncrease
model classes is G4DNADingfelderChargeIncreaseModel
Charge decrease
process class is G4DNAChargeDecrease
model classes is G4DNADingfelderChargeDecreaseModel
Helium++ (ionised twice) processes and models
Elastic scattering :
process class is G4DNAElastic
G4DNAIonElasticModel
Excitation
process class is G4DNAExcitation
model classes is G4DNAMillerGreenExcitationModel
Ionisation
process class is G4DNAIonisation
model classes is G4DNARuddIonisationModel
Charge decrease
process class is G4DNAChargeDecrease
model classes is G4DNADingfelderChargeDecreaseModel
Li, Be, B, C, N, O, Si, Fe processes and models
Ionisation
process class is G4DNAIonisation
model class is G4DNARuddIonisationExtendedModel
Examples of the registration of these processes in a physics list are
given in the G4EmDNAPhysics* constructors (in
source/physics_lists/constructors/electromagnetic
in the Geant4
source distribution). An example of the usage of these constructors in a
physics list is given in the “dnaphysics” extended example, which
explains how to extract basic information from Geant4-DNA Physics
processes.
Geant4-DNA physics constructors are described at the Geant4-DNA website.
The “microdosimetry” extended example illustrates how to combine Geant4-DNA processes with Standard electromagnetic processes (combination of discrete and condensed history Geant4 electromagnetic processes at different scales).
A set of Geant4-DNA models applicable to biological materials is available since release 10.4. These models are named G4DNAPTBElasticModel, G4DNAPTBExcitationModel, G4DNAPTBIonisationModel and G4DNAPTBAugerModel. They can be used for electrons in THF, PY, PU, TMP precursors and in backbone, cytosine, thymine, adenine, guanine materials of DNA. The G4DNAPTBIonisationModel can also be used with protons in THF, PY and TMP. Their usage is illustrated in the “icsd” extended example.
Since Geant4 release 10.1, Geant4-DNA can also be used for the modelling of water radiolysis (physico-chemistry and chemistry stages). Three extended examples, “chem1”, “chem2”, “chem3” and “chem4” illustrate this. More information is available from the Geant4-DNA website.
To run the Geant4-DNA extension, data files need to be copied by the user to his/her code repository. These files are distributed together with the Geant4 release. The user should set the environment variable G4LEDATA to the directory where he/she has copied the files.
A full list of publications regarding Geant4-DNA is directly available from the Geant4-DNA website or from the Geant4@IN2P3 web site).
Atomic Deexcitation¶
A unique interface named G4VAtomicDeexcitation is available in Geant4 for the simulation of atomic deexcitation using Standard, Low Energy and Very Low Energy electromagnetic processes. Atomic deexcitation includes fluorescence and Auger electron emission induced by photons, electrons and ions (PIXE); see more details in:
A. Mantero et al., PIXE Simulation in Geant4 , X-Ray Spec. , 40, 135-140, 2011.
It can be activated for processes producing vacancies in atomic shells. Currently these processes are the photoelectric effect, ionization and Compton scattering.
Activation of atomic deexcitation
The activation of atomic deexcitation in continuous processes in a user physics list can be done through the following G4EmParameters class methods described above or via UI commands
/process/em/deexcitation region true true true
/process/em/fluo true
/process/em/auger true
/process/em/pixe true
One can define parameters in the G4State_PreInit or G4State_Idle states. Fluorescence from photons and electrons is activated by default in Option3, Option4, Livermore and Penelope physics constructors, while Auger production and PIXE are not.
The alternative set of data by Bearden et al. (1967) for the modelling of fluorescence lines had been added to the G4LEDATA archive. This set can be selected via UI command
/process/em/fluoBearden true
Another important UI commands enable simulation of the full Auger and/or fluorescence cascade
/process/em/deexcitationIgnoreCut true
How to change ionisation cross section models ?
The user can also select which cross section model to use in order to calculate shell ionisation cross sections for generating PIXE
/process/em/pixeXSmodel name
/process/em/pixeElecXSmodel name
where the name can be “Empirical”, “ECPSSR_FormFactor” or “ECPSSR_Analytical” corresponds to different PIXE cross sections. Following shell cross sections models are available : “ECPSSR_Analytical” models derive from an analytical calculation of the ECPSSR theory (see A. Mantero et al., X-Ray Spec.40 (2011) 135-140) and it reproduces K and L shell cross sections over a wide range of energies; “ECPSSR_FormFactor” models derive from A. Taborda et al. calculations (see A. Taborda et al., X-Ray Spec. 40 (2011) 127-134) of ECPSSR values directly form Form Factors and it covers K, L shells on the range 0.1-100 MeV and M shells in the range 0.1-10 MeV; the “empirical” models are from Paul “reference values” (for protons and alphas for K-Shell) and Orlic empirical model for L shells (only for protons and ions with Z>2). The later ones are the models used by default. Out of the energy boundaries, “ECPSSR_Analytical” model is used. We recommend to use default settings if not sure what to use.
Example
The TestEm5 extended/electromagetic example shows how to simulate atomic deexcitation (see for e.g. the pixe.mac macro).
Very Low energy Electromagnetic Processes in Silicon for microelectronics application (Geant4-MuElec extension)¶
(Previously named Geant4-MuElec)
The Geant4 low energy electromagnetic Physics package has been extended down to energies of a few electron Volts suitable for the simulation of radiation effects in highly integrated microelectronic components.
The Geant4-MicroElec process and model classes apply to electrons, protons and heavy ions in silicon.
Electron processes and models
Elastic scattering :
process class is G4MicroElastic
model class is G4MicroElecElasticModel
Ionization
process class is G4MicroElecInelastic
model class is G4MicroElecInelasticModel
Proton processes and models
Ionisation
process class is G4MicroElecInelastic
model class is G4MicroElecInelasticModel
Heavy ion processes and models
Ionization
process class is G4MicroElecInelastic
model class is G4MicroElecInelasticModel
A full list of publications regarding Geant4-MicroElec is directly available from the Geant4-MicroElec website.
New Compton model by Monash U., Australia¶
A new Compton scattering model for unpolarised photons has been developed in the relativistic impulse approximation. The model was developed as an alternative to low energy electromagnetic Compton scattering models developed from Ribberfors’ Compton scattering framework (Livermore, Penelope Compton models). The model class is named named G4LowEPComptonModel.
G4LowEPComptonModel has been added to the physics constructor G4EmStandardPhysics_option4, containing the most accurate models from the Standard and Low Energy Electromagnetic physics working groups.
Multi-scale Processes¶
Hadron Impact Ionisation and PIXE¶
The G4hImpactIonisation process deals with ionisation by impact of hadrons and alpha particles, and the following generation of PIXE (Particle Induced X-ray Emission). This process and related classes can be found in source/processes/electromagnetic/pii.
Further documentation about PIXE simulation with this process is available here.
A detailed description of the related physics features can be found in:
A brief summary of the related physics features can be found in the Geant4 Physics Reference Manual.
An example of how to use this process is shown below. A more extensive example is available in the eRosita Geant4 advanced example (see examples/advanced/eRosita in your Geant4 installation source).
#include "G4hImpactIonisation.hh"
[...]
void eRositaPhysicsList::ConstructProcess()
{
[...]
theParticleIterator->reset();
while( (*theParticleIterator)() )
{
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* processManager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (particleName == "proton")
{
// Instantiate the G4hImpactIonisation process
G4hImpactIonisation* hIonisation = new G4hImpactIonisation();
// Select the cross section models to be applied for K, L and M shell vacancy creation
// (here the ECPSSR model is selected for K, L and M shell; one can mix and match
// different models for each shell)
hIonisation->SetPixeCrossSectionK("ecpssr");
hIonisation->SetPixeCrossSectionL("ecpssr");
hIonisation->SetPixeCrossSectionM("ecpssr");
// Register the process with the processManager associated with protons
processManager -> AddProcess(hIonisation, -1, 2, 2);
}
}
}
Available cross section model options
The following cross section model options are available:
protons
K shell
ecpssr
(based on the ECPSSR theory)ecpssr_hs
(based on the ECPSSR theory, with Hartree-Slater correction)ecpssr_ua
(based on the ECPSSR theory, with United Atom Hartree-Slater correction)ecpssr_he
(based on the ECPSSR theory, with high energy correction)pwba
(plane wave Born approximation)paul
(based on the empirical model by Paul and Sacher)kahoul
(based on the empirical model by Kahoul et al.)
L shell
ecpssr
ecpssr_ua
pwba
miyagawa
(based on the empirical model by Miyagawa et al.)orlic
(based on the empirical model by Orlic et al.)sow
(based on the empirical model by Sow et al.)
M shell
ecpssr
pwba
alpha particles
K shell
ecpssr
ecpssr_hs
pwba
paul
(based on the empirical model by Paul and Bolik)
L shell
ecpssr
pwba
M shell
ecpssr
pwba
PIXE data library
The G4hImpactIonisation process uses a PIXE Data Library.
The PIXE Data Library is distributed in the Geant4 G4PII data set, which must be downloaded along with Geant4 source code.
The G4PIIDATA environment variable must be defined to refer to the location of the G4PII PIXE data library in your filesystem; for instance, if you use a c-like shell
setenv G4PIIDATA path_to_where_G4PII_has_been_downloaded
Further documentation about the PIXE Data Library is available here.
Hadronic Interactions¶
This section briefly introduces the hadronic physics processes installed in Geant4. For details of the implementation of hadronic interactions available in Geant4, please refer to the Physics Reference Manual.
Treatment of Cross Sections¶
Cross section data sets¶
Each hadronic process object (derived from G4HadronicProcess
) may
have one or more cross section data sets associated with it. The term
“data set” is meant, in a broad sense, to be an object that encapsulates
methods and data for calculating total cross sections for a given
process. The methods and data may take many forms, from a simple
equation using a few hard-wired numbers to a sophisticated
parameterisation using large data tables. Cross section data sets are
derived from the abstract class G4VCrossSectionDataSet
, and are
required to implement the following methods:
G4bool IsApplicable( const G4DynamicParticle*, const G4Element* )
This method must return True
if the data set is able to calculate a
total cross section for the given particle and material, and False
otherwise.
G4double GetCrossSection( const G4DynamicParticle*, const G4Element* )
This method, which will be invoked only if True
was returned by
IsApplicable
, must return a cross section, in Geant4 default units,
for the given particle and material.
void BuildPhysicsTable( const G4ParticleDefinition& )
This method may be invoked to request the data set to recalculate its internal database or otherwise reset its state after a change in the cuts or other parameters of the given particle type.
void DumpPhysicsTable( const G4ParticleDefinition& ) = 0
This method may be invoked to request the data set to print its internal database and/or other state information, for the given particle type, to the standard output stream.
Cross section data store¶
Cross section data sets are used by the process for the calculation of
the physical interaction length. A given cross section data set may only
apply to a certain energy range, or may only be able to calculate cross
sections for a particular type of particle. The class
G4CrossSectionDataStore
has been provided to allow the user to
specify, if desired, a series of data sets for a process, and to arrange
the priority of data sets so that the appropriate one is used for a
given energy range, particle, and material. It implements the following
public methods:
G4CrossSectionDataStore()
~G4CrossSectionDataStore()
and
G4double GetCrossSection( const G4DynamicParticle*, const G4Element* )
For a given particle and material, this method returns a cross section
value provided by one of the collection of cross section data sets
listed in the data store object. If there are no known data sets, a
G4Exception
is thrown and DBL_MIN
is returned. Otherwise, each
data set in the list is queried, in reverse list order, by invoking its
IsApplicable
method for the given particle and material. The first
data set object that responds positively will then be asked to return a
cross section value via its GetCrossSection
method. If no data set
responds positively, a G4Exception
is thrown and DBL_MIN
is
returned.
void AddDataSet( G4VCrossSectionDataSet* aDataSet )
This method adds the given cross section data set to the end of the list
of data sets in the data store. For the evaluation of cross sections,
the list has a LIFO (Last In First Out) priority, meaning that data sets
added later to the list will have priority over those added earlier to
the list. Another way of saying this, is that the data store, when given
a GetCrossSection
request, does the IsApplicable
queries in the
reverse list order, starting with the last data set in the list and
proceeding to the first, and the first data set that responds positively
is used to calculate the cross section.
void BuildPhysicsTable( const G4ParticleDefinition& aParticleType )
This method may be invoked to indicate to the data store that there has
been a change in the cuts or other parameters of the given particle
type. In response, the data store will invoke the BuildPhysicsTable
of each of its data sets.
void DumpPhysicsTable( const G4ParticleDefinition& )
This method may be used to request the data store to invoke the
DumpPhysicsTable
method of each of its data sets.
Default cross sections¶
The defaults for total cross section data and calculations have been
encapsulated in the singleton class G4HadronCrossSections
. Each
hadronic process: G4HadronInelasticProcess
,
G4HadronElasticProcess
, G4HadronFissionProcess
, and
G4HadronCaptureProcess
, comes already equipped with a cross section
data store and a default cross section data set. The data set objects
are really just shells that invoke the singleton
G4HadronCrossSections
to do the real work of calculating cross
sections.
The default cross sections can be overridden in whole or in part by the
user. To this end, the base class G4HadronicProcess
has a get
method:
G4CrossSectionDataStore* GetCrossSectionDataStore()
which gives public access to the data store for each process. The user’s cross section data sets can be added to the data store according to the following framework:
G4Hadron...Process aProcess(...)
MyCrossSectionDataSet myDataSet(...)
aProcess.GetCrossSectionDataStore()->AddDataSet( &MyDataSet )
The added data set will override the default cross section data whenever
so indicated by its IsApplicable
method.
In addition to the get
method, G4HadronicProcess
also has the
method
void SetCrossSectionDataStore( G4CrossSectionDataStore* )
which allows the user to completely replace the default data store with a new data store.
It should be noted that a process does not send any information about
itself to its associated data store (and hence data set) objects. Thus,
each data set is assumed to be formulated to calculate cross sections
for one and only one type of process. Of course, this does not prevent
different data sets from sharing common data and/or calculation methods,
as in the case of the G4HadronCrossSections
class mentioned above.
Indeed, G4VCrossSectionDataSet
specifies only the abstract interface
between physics processes and their data sets, and leaves the user free
to implement whatever sort of underlying structure is appropriate.
The current implementation of the data set G4HadronCrossSections
reuses the total cross-sections for inelastic and elastic scattering,
radiative capture and fission as used with GHEISHA to provide
cross-sections for calculation of the respective mean free paths of a
given particle in a given material.
Cross-sections for low energy neutron transport¶
The cross section data for low energy neutron transport are organized in
a set of files that are read in by the corresponding data set classes at
time zero. Hereby the file system is used, in order to allow highly
granular access to the data. The ``root’’ directory of the
cross-section directory structure is accessed through an environment
variable, G4NEUTRONHPDATA
, which is to be set by the user. The
classes accessing the total cross-sections of the individual processes,
i.e., the cross-section data set classes for low energy neutron
transport, are G4NeutronHPElasticData
, G4NeutronHPCaptureData
,
G4NeutronHPFissionData
, and G4NeutronHPInelasticData
.
For detailed descriptions of the low energy neutron total cross-sections, they may be registered by the user as described above with the data stores of the corresponding processes for neutron interactions.
It should be noted that using these total cross section classes does not require that the neutron_hp models also be used. It is up to the user to decide whether this is desirable or not for his particular problem.
A prototype of the compact version of neutron and light ion
cross sections derived
from HP database are provided with classes
G4NeutronHPElasticData
, G4NeutronCaptureXS
,
G4NeutronElasticXS
, G4NeutronInelasticXS
, and G4ParticleInelasticXS
.
These cross-section directory structure is accessed through an environment
variable, G4PARTICLEXSDATA
Cross-sections for low-energy charged particle transport¶
The cross-section data for low-energy charged particle transport are
organized in a set of files that are read at initialization, similarly
to the case of low-energy neutron transport. The “root” directory of the
cross-section directory structure is accessed through an environment
variable, G4PARTICLEHPDATA
, which has to be set by the user. This
variable has to point to the directory where the low-energy charged
particle data have been installed, e.g. G4TENDL1.3
for the Geant4
release 10.3
(note that the download of this data library from the
Geant4 web site is not done automatically, i.e. it must be done manually
by the user):
export G4PARTICLEHPDATA=/your/path/G4TENDL1.3/
It is expected that the directory $G4PARTICLEHPDATA
has the
following five subdirectories, corresponding to the charged particles
that can be handled by the low-energy charged particle transport:
Proton/
, Deuteron/
, Triton/
, He3/
, Alpha/
. It is
possible for the user to overwrite the default directory structure with
individual environment variables pointing to custom data libraries for
each particle type. This is considered an advanced/expert user feature.
These directories are set by the following environment variables:
G4PROTONHPDATA
, for proton; G4DEUTERONHPDATA
, for deuteron;
G4TRITONHPDATA
, for triton; G4HE3HPDATA
, for He3;
G4ALPHAHPDATA
, for alpha. Note that if any of these variables is not
defined explicitly, e.g. G4TRITONHPDATA
, then the corresponding data
library is expected to be a subdirectory of $G4PARTICLEHPDATA/
, e.g.
$G4PARTICLEHPDATA/Triton/
. If instead all the above five
environmental variables are set, then G4PARTICLEHPDATA
does not need
to be specified; even if it is set, then its value will be ignored
(because the per-particle ones take precedence).
Hadrons at Rest¶
List of implemented “Hadron at Rest” processes¶
The following process classes have been implemented:
\(\pi^-, K^-, \sigma^-, \xi^-, \omega^-\) absorption (class name
G4HadronicAbsorptionBertini
)neutron capture (class name
G4HadronCaptureProcess
)anti-proton, anti-\(\sigma^+\), anti-deuteron, anti-triton, anti-alpha, anti-He3 annihilation (class name
G4HadronicAbsorptionFritiof
)mu- capture (class name
G4MuonMinusCapture
)
Capture of low-energy negatively charged particles is a complex process
involving formation of mesonic atoms, X-ray cascade and Auger cascade,
nuclear interaction. In the case of mu- there is also a probability to
decay from K-shell of mesonic atom. To handle this a base process class
G4HadronicStoppingProcess
is used.
For the case of neutrons, Geant4 offer simulation down to thermal energies. The capture cross section generally increases when neutron energy decreases and there are many nuclear resonances. In Geant4 neutron capture cross sections are parameterized using ENDF database.
Hadrons in Flight¶
What processes do you need?¶
For hadrons in motion, there are four physics process classes. Table 8 shows each process and the particles for which it is relevant.
|
\(\pi^+, \pi^-, K^+, K^0_S, K^0_L, K^-, p, \bar{p}, n, \bar{n}, \lambda, \bar{\lambda}, \Sigma^+, \Sigma^-, \bar{\Sigma}^+, \bar{\Sigma}^-, \Xi^0, \Xi^-, \bar{\Xi}^0, \bar{\Xi}^- \) |
|
\(\pi^+, \pi^-, K^+, K^0_S, K^0_L, K^-, p, \bar{p}, n, \bar{n}, \lambda, \bar{\lambda}, \Sigma^+, \Sigma^-, \bar{\Sigma}^+, \bar{\Sigma}^-, \Xi^0, \Xi^-, \bar{\Xi}^0, \bar{\Xi}^- \) |
|
all |
|
\(n, \bar{n}\) |
How to register Models¶
To register an inelastic process model for a particle, a proton for example, first get the pointer to the particle’s process manager:
G4ParticleDefinition *theProton = G4Proton::ProtonDefinition();
G4ProcessManager *theProtonProcMan = theProton->GetProcessManager();
Create an instance of the particle’s inelastic process:
G4ProtonInelasticProcess *theProtonIEProc = new G4ProtonInelasticProcess();
Create an instance of the model which determines the secondaries produced in the interaction, and calculates the momenta of the particles, for instance the Bertini cascade model:
G4CascadeInterface *theProtonIE = new G4CascadeInterface();
Register the model with the particle’s inelastic process:
theProtonIEProc->RegisterMe( theProtonIE );
Finally, add the particle’s inelastic process to the list of discrete processes:
theProtonProcMan->AddDiscreteProcess( theProtonIEProc );
The particle’s inelastic process class, G4ProtonInelasticProcess
in
the example above, derives from the G4HadronicInelasticProcess
class, and simply defines the process name and calls the
G4HadronicInelasticProcess
constructor. All of the specific particle
inelastic processes derive from the G4HadronicInelasticProcess
class, which calls the PostStepDoIt
function, which returns the
particle change object from the G4HadronicProcess
function
GeneralPostStepDoIt
. This class also gets the mean free path, builds
the physics table, and gets the microscopic cross section. The
G4HadronicInelasticProcess
class derives from the
G4HadronicProcess
class, which is the top level hadronic process
class. The G4HadronicProcess
class derives from the
G4VDiscreteProcess
class. The inelastic, elastic, capture, and
fission processes derive from the G4HadronicProcess
class. This pure
virtual class also provides the energy range manager object and the
RegisterMe
access function.
In-flight, final-state hadronic models derive, directly or indirectly,
from the G4InelasticInteraction
class, which is an abstract base
class since the pure virtual function ApplyYourself
is not defined
there. G4InelasticInteraction
itself derives from the
G4HadronicInteraction
abstract base class. This class is the base
class for all the model classes. It sorts out the energy range for the
models and provides class utilities. The G4HadronicInteraction
class
provides the Set/GetMinEnergy
and the Set/GetMaxEnergy
functions
which determine the minimum and maximum energy range for the model. An
energy range can be set for a specific element, a specific material, or
for general applicability:
void SetMinEnergy( G4double anEnergy, G4Element *anElement )
void SetMinEnergy( G4double anEnergy, G4Material *aMaterial )
void SetMinEnergy( const G4double anEnergy )
void SetMaxEnergy( G4double anEnergy, G4Element *anElement )
void SetMaxEnergy( G4double anEnergy, G4Material *aMaterial )
void SetMaxEnergy( const G4double anEnergy )
Which models are there, and what are the defaults¶
In Geant4, any model can be run together with any other model without the need for the implementation of a special interface, or batch suite, and the ranges of applicability for the different models can be steered at initialisation time. This way, highly specialised models (valid only for one material and particle, and applicable only in a very restricted energy range) can be used in the same application, together with more general code, in a coherent fashion.
Each model has an intrinsic range of applicability, and the model chosen for a simulation depends very much on the use-case. Consequently, there are no “defaults”. However, physics lists are provided which specify sets of models for various purposes.
Two types of hadronic shower models have been implemented: data driven models and theory driven models.
Data driven models are available for the transport of low energy neutrons in matter in sub-directory
hadronics/models/neutron_hp
. The modeling is based on the data formats of ENDF/B-VI, and all distributions of this standard data format are implemented. The data sets used are selected from data libraries that conform to these standard formats. The file system is used in order to allow granular access to, and flexibility in, the use of the cross sections for different isotopes, and channels. The energy coverage of these models is from thermal energies to 20 MeV.Theory driven models are available for inelastic scattering in a first implementation, covering the full energy range of LHC experiments. They are located in sub-directory
hadronics/models/generator
. The current philosophy implies the usage of parton string models at high energies, of intra-nuclear transport models at intermediate energies, and of statistical break-up models for de-excitation.
High-precision neutron interactions (NeutronHP)¶
Nuclear models fail (sometimes catastrophically) at predicting with reasonable accuracies the nuclear cross sections of neutrons (and other particles). For this reason, all physical quantities relevant for an accurate modeling of nuclear reactions in Monte Carlo simulations need to be provided as a database which includes, ideally:
cross sections
angular distributions of the emitted particles
energy spectra of the emitted particles
energy-angle correlated spectrum (double-differential cross sections, DDX)
neutrons per fission
fission spectra
fission product yields
photo production data
For the case of neutron induced reactions, such databases are called “evaluated data”, in the sense that they contain recommended values for different quantities that rely on compilations of experimental nuclear data and usually completed with theoretical predictions, benchmarked against available experimental data (i.e. integral and differential experiments) when possible. It should be noticed that the information available varies from isotope to isotope and can be incomplete or totally missing.
The G4NeutronHP package in Geant4 allows using evaluated nuclear data libraries in the G4NDL format. Geant4 users should know that any simulation involving neutrons with energies below 20 MeV and not using the G4NeutronHP package can lead to unreliable results. Geant4 users are therefore encouraged to use it, although they should be aware of the limitations of using evaluated nuclear data libraries.
An example about how to implement the G4NeutronHP package into physics
list in a Geant4 application can be found in the example case (among
others distributed with Geant4) extended/radioactivedecay/rdecay02
.
Three different processes are included in that example: elastic, capture
and inelastic. The inelastic reactions in G4NeutronHP are all reactions
except elastic, capture and fission, so fission should also be included
in the physics list, if needed, and it is done in the same way as it is
done for the other three.
The G4NeutronHP package must be used together with evaluated nuclear data libraries. They are available on the Geant4 download page and from the IAEA nuclear data web site where a larger set of different libraries, including isotopes with Z > 92, is available.
The evaluated nuclear data libraries do differ and thus the results of the Monte Carlo simulations will depend on the library used. It is a safe practice to perform simulations with (at least) two different libraries for estimating the uncertainties associated to the nuclear data.
Together with a good implementation of the physics list, users must be very careful with the definition of the materials performed in a Monte Carlo simulation when low energy neutron transport is relevant. In contrast to other kind of simulations, the isotopic composition of the elements which compose the different materials can strongly affect the obtained simulation results. Because of this, it is strongly recommended to define specifically the isotopic composition of each element used in the simulation, as it is described in the Geant4 user’s manual. In principle, such a practice is not mandatory if natural isotopic compositions are used, since Geant4 contains them in their databases. However, by defining them explicitly some unexpected problems may be avoided and a better control of the simulation will be achieved.
It is highly recommended or mandatory to set the following UNIX environment variables running a Geant4 application:
G4NEUTRONHPDATA
[path to the G4NDL format data libraries] (mandatory).
G4NEUTRONHP_SKIP_MISSING_ISOTOPES=1
It sets to zero the cross section of the isotopes which are not present in the neutron library. If Geant4 doesn’t find an isotope, then it looks for the natural composition data of that element. Only if the element is not found then the cross section is set to zero. On the contrary, if this variable is not defined, Geant4 looks then for the neutron data of another isotope close in Z and A, which will have completely different nuclear properties and lead to incorrect results (highly recommended).
G4NEUTRONHP_DO_NOT_ADJUST_FINAL_STATE=1
If this variable is not defined, a Geant4 model that attempts to satisfy the energy and momentum conservation in some nuclear reactions, by generating artificial gamma rays. By setting such a variable one avoids the correction and leads to the result obtained with the ENDF-6 libraries. Even though energy and momentum conservation are desirable, the ENDF-6 libraries do not provide the necessary correlations between secondary particles for satisfying them in all cases. On the contrary, ENDF-6 libraries intrinsically violate energy and momentum conservation for several processes and have been built for preserving the overall average quantities such as average energy releases, average number of secondaries… (highly recommended).
AllowForHeavyElements=1
Activates the physics for isotopes with Z>92 (recommended).
The G4NDL format libraries are based on the ENDF-6 format libraries, which contain evaluated (i.e. recommended) nuclear data prepared for their use in transport codes. These data are essentially nuclear reaction cross sections together with the distribution in energy and angle of the secondary reaction products. As a consequence of how the data is written in the ENDF files, there are some features that may be or may be not expected in the results of a Monte Carlo calculation.
The information concerning the creation of the reaction products can be incomplete and/or uncorrelated, in the sense that is described below:
Incomplete information.
This applies when there is no information about how to generate a secondary particle. As an example, it is possible to have only the cross section data of an (n,p) reaction, without any information concerning the energy and angle of the secondary proton. In this case Geant4 will produce the proton considering that it is emitted isotropically in the center of mass frame, with an energy which is deduced from assuming that the residual nucleus is in its ground state.
As a consequence of what has been presented above, some general features can be expected in the results of a Monte Carlo calculation performed with the G4NeutronHP package:
The neutron transport, which means how the neutron looses energy in the collisions, when and how it is absorbed…, is quite trustable, since the main purpose of the ENDF neutron libraries is to perform this neutron transport.
The production of neutrons due to neutron induced nuclear reactions is usually trustable, with the exception of the energy-angle correlations when several neutrons are produced in the same nuclear reaction.
The results concerning the production of charged particles have to be always questioned. A look into the ENDF format library used can indicate which results are trustable and which are not. This can be done, for example, in http://t2.lanl.gov/nis/data.html, among other websites.
The results concerning the production of \(\gamma\)-rays have to be questioned always. For example, the information on the number and energies of \(\gamma\)-rays emitted in the neutron capture process is incomplete for almost all the nuclei and is frequently also uncorrelated. When the information is available, it will be used, but one can obtain results which are quite far from reality on an event by event basis: the total energy of the cascade won’t be correct in many cases and only some specific \(\gamma\)-rays which are stored in the neutron databases will be emitted. If there isn’t any information concerning these \(\gamma\)-rays, Geant4 will use a simple a model instead which is generally missing the relevant spectroscopic information. The results concerning the generation of residual nuclei (for example, in activation calculations) are usually trustable, with the exception of libraries with MT=5 reactions, as described above (uncorrelated).
As a general conclusion, users should always be critical with the results obtained with Monte Carlo simulation codes, and this also applies to Geant4. They have to anticipate which results can be trusted and which results should be questioned. For the particular case of the a closer look into the underlying evaluated nuclear data in the ENDF format libraries will allow to check what is the information available in a certain library for some specific isotope and a certain reaction. There are several public nuclear data sites like http://t2.lanl.gov/nis/data.html.
The transport of very low energy neutrons (below 5 eV) has to be performed using the thermal neutron data libraries. At these energies, the fact that the nuclei are in atoms which form part of a certain molecule inside a material (crystal lattice, liquid, plastic…) plays an important role, since there can be a transference of momentum between the neutron and the whole structure of the material, not only with the nucleus. This is of particular importance for material used as neutron moderators, i.e., materials with low A (mass number) used to decrease the incident neutron energy in only a few collisions. Since the property is related to the nucleus in the material, as an example, there is the need for having different thermal libraries for Hydrogen in polyethylene, Hydrogen in water and so on.
If neutron collisions at these energies are relevant for the problem to be simulated, thermal libraries should be used for the materials if they are available. If they are not, the results obtained from the simulation will not be trustable in the neutron energy range below 5 eV, especially when using low mass elements in the simulation.
To use the thermal libraries the following lines should be included in the physics list:
G4HadronElasticProcess* theNeutronElasticProcess = new G4HadronElasticProcess;
// Cross Section Data set
G4NeutronHPElasticData* theHPElasticData = new G4NeutronHPElasticData;
theNeutronElasticProcess->AddDataSet(theHPElasticData);
G4NeutronHPThermalScatteringData* theHPThermalScatteringData = new G4NeutronHPThermalScatteringData;
theNeutronElasticProcess->AddDataSet(theHPThermalScatteringData);
// Models
G4NeutronHPElastic* theNeutronElasticModel = new G4NeutronHPElastic;
theNeutronElasticModel->SetMinEnergy(4.0*eV);
theNeutronElasticProcess->RegisterMe(theNeutronElasticModel);
G4NeutronHPThermalScattering* theNeutronThermalElasticModel = new G4NeutronHPThermalScattering;
theNeutronThermalElasticModel->SetMaxEnergy(4.0*eV);
theNeutronElasticProcess->RegisterMe(theNeutronThermalElasticModel);
// Apply Processes to Process Manager of Neutron
G4ProcessManager* pmanager = G4Neutron::Neutron()->GetProcessManager();
pmanager->AddDiscreteProcess(theNeutronElasticProcess);
And the materials should be defined with a specific name. For example, to use the thermal library for Hydrogen in water, the water should be defined as:
G4Element* elTSHW = new G4Element("TS_H_of_Water", "H_WATER", 1.0, 1.0079*g/mole);
G4Material* matH2O_TS = new G4Material("Water_TS", density=1.0*g/cm3, ncomponents=2);
matH2O_TS->AddElement(elTSHW,natoms=2);
matH2O_TS->AddElement(elO,natoms=1);
where the important thing is the name “TS_H_of_Water”, which is a
specific name used by G4NeutronHP. In order to see which thermal
libraries are available, they can be found in the
G4NDL4.0/ThermalScattering folder (or equivalent, for other neutron
libraries). Then, one has to look into the
G4NeutronHPThermalScatteringNames.cc
source file, under
source/processes/hadronic/models/neutron_hp/src
. There are some
lines similar to:
names.insert(std::pair<G4String,G4String>("TS_H_of_Water", "h_water"));
where “TS_H_of_Water” means Hydrogen in water. Names similar to “TS_H_of_Water” like “TS_C_of_Graphite” or “TS_H_of_Polyethylene” can be found and used in the same way as described above.
High-precision charged particle interactions (ParticleHP)¶
Due to the coupling between the configuration for neutrons and charged particles in ParticleHP, the default one is not the recommended one from the physics point of view for all particles. A consistent configuration with thorough testing will hopefully be introduced in the next release. For the time being, in order to improve the physics performance for primary charged particles the following environment variable should be set
export DO_NOT_SET_PHP_AS_HP=1
Note that this environmental variable is a configuration option which is used only at compilation, not at run time, and it affects both primary neutrons and charged particles. It is not expected to dramatically change the behaviour for neutrons.
For further improvement with projectile charged particles, it is also recommended to set the following environmental variable used at run-time
export G4PHP_DO_NOT_ADJUST_FINAL_STATE=1
which avoids the default adjustment of the final state to ensure better conservation laws (for charge, energy, momentum, baryon number).
The adjustment of the final state is recommended for realistic detector response in the case of neutron interactions. For the use-case of reactor physics and dosimetry, where average quantities are important, not adjusting the final state (i.e. setting the above environment variable) improves accuracy.
Note that, for the time being, setting
G4PHP_DO_NOT_ADJUST_FINAL_STATE
affects both primary neutrons and
charged particles, so be careful which is the use-case you are
interested in. To summarize: if you use ParticleHP for primary neutrons,
you can safely take the default; no harm is expected if you build
ParticleHP with DO_NOT_SET_PHP_AS_HP
set; be very careful instead if
you set G4PHP_DO_NOT_ADJUST_FINAL_STATE
. If you use ParticleHP for
primary charged particles, then it is recommended to build with
DO_NOT_SET_PHP_AS_HP
set, and then run with DO_NOT_SET_PHP_AS_HP
set.
Switching statistical nuclear de-excitation models¶
Nuclear reactions at intermediate energies (from a few MeV to a few GeV) are typically modelled in two stages. The first, fast reaction stage is described by a dynamical model (quantum molecular dynamics, intranuclear cascade, pre-compound, etc.) and often results in the production of one or several excited nuclei. The second reaction stage describes the de-excitation of the excited nuclei and it is usually handled by statistical de-excitation models. The models for the two reaction stages can in principle be chosen independently, but the current design of the Geant4 hadronics framework makes it difficult to do this at the physics-list level. However, another solution exists.
Geant4 provides several nuclear de-excitation modules. The default one
is G4ExcitationHandler
, which is described in detail in the
Physics Reference Manual.
The Bertini-style G4CascadeInterface
uses an internal de-excitation
model. The ABLA V3 model is also available.
Options are available for steering of the pre-compound model and the
de-excitation module. These options may be invoked by the new C++
interface class G4DeexPrecoParameters
. The interface
G4NuclearLevelData::Instance()->GetParameters()
is thread safe,
parameters are shared between threads, and parameters are shared between
all de-excitation and pre-compound classes. Parameters may be modified
at G4State_PreInit state of Geant4. This class has the following public
methods:
Dump()
StreamInfo(std::ostream&)
SetLevelDensity(G4double)
SetR0(G4double)
SetTransitionsR0(G4double)
SetFBUEnergyLimit(G4double)
SetFermiEnergy(G4double)
SetPrecoLowEnergy(G4double)
SetPrecoHighEnergy(G4double)
SetPhenoFactor(G4double)
SetMinExcitation(G4double)
SetMaxLifeTime(G4double)
SetMinExPerNucleounForMF(G4double)
SetMinEForMultiFrag(G4double)
SetMinZForPreco(G4int)
SetMinAForPreco(G4int)
SetPrecoModelType(G4int)
SetDeexModelType(G4int)
SetTwoJMAX(G4int)
SetUploadZ(G4int)
SetVerbose(G4int)
SetNeverGoBack(G4bool)
SetUseSoftCutoff(G4bool)
SetUseCEM(G4bool)
SetUseGNASH(G4bool)
SetUseHETC(G4bool)
SetUseAngularGen(G4bool)
SetPrecoDummy(G4bool)
SetCorrelatedGamma(G4bool)
SetStoreICLevelData(G4bool)
SetInternalConversionFlag(G4bool)
SetLevelDensityFlag(G4bool)
SetDiscreteExcitationFlag(G4bool)
SetDeexChannelType(G4DeexChannelType)
It is possible to replace the default de-excitation model with ABLA V3
for any intranuclear-cascade model in Geant4 except
G4CascadeInterface
. The easiest way to do this is to call the
SetDeExcitation()
method of the relevant intranuclear-cascade-model
interface. This can be done even if you are using one of the reference
physics lists. The technique is the following.
For clarity’s sake, assume you are using the FTFP_INCLXX physics
list, which uses INCL++, the Liege Intranuclear Cascade model
(G4INCLXXInterface
) at intermediate energies. You can couple
INCL++ to ABLA V3 by adding a run action
(Usage of User Actions) and adding the following code snippet to
BeginOfRunAction()
.
#include "G4HadronicInteraction.hh"
#include "G4HadronicInteractionRegistry.hh"
#include "G4INCLXXInterface.hh"
#include "G4AblaInterface.hh"
void MyRunAction::BeginOfRunAction(const G4Run*)
{
// Get hold of pointers to the INCL++ model interfaces
std::vector<G4HadronicInteraction *> interactions = G4HadronicInteractionRegistry::Instance()
->FindAllModels(G4INCLXXInterfaceStore::GetInstance()->getINCLXXVersionName());
for(std::vector<G4HadronicInteraction *>::const_iterator iInter=interactions.begin(), e=interactions.end();
iInter!=e; ++iInter) {
G4INCLXXInterface *theINCLInterface = static_cast<G4INCLXXInterface*>(*iInter);
if(theINCLInterface) {
// Instantiate the ABLA model
G4HadronicInteraction *interaction = G4HadronicInteractionRegistry::Instance()->FindModel("ABLA");
G4AblaInterface *theAblaInterface = static_cast<G4AblaInterface*>(interaction);
if(!theAblaInterface)
theAblaInterface = new G4AblaInterface;
// Couple INCL++ to ABLA
G4cout << "Coupling INCLXX to ABLA" << G4endl;
theINCLInterface->SetDeExcitation(theAblaInterface);
}
}
}
This technique may be applied to any intranuclear-cascade model (i.e.
models that inherit from G4VIntraNuclearTransportModel
), except
G4CascadeInterface
. For example, if your physics list relies on the
Binary-Cascade model (e.g. FTF_BIC), you’ll need to do
// Get hold of a pointer to the Binary-Cascade model interface
std::vector<G4HadronicInteraction *> interactions = G4HadronicInteractionRegistry::Instance()
->FindAllModels("Binary Cascade");
for(std::vector<G4HadronicInteraction *>::const_iterator iInter=interactions.begin(), e=interactions.end();
iInter!=e; ++iInter) {
G4BinaryCascade *theBICInterface = static_cast<G4BinaryCascade*>(*iInter);
if(theBICInterface) {
// Instantiate ABLA V3 as in the example above
// [...]
// Couple BIC to ABLA
theBICInterface->SetDeExcitation(theAblaInterface);
}
}
Particle Decay Process¶
This section briefly introduces decay processes installed in Geant4. For details of the implementation of particle decays, please refer to the Physics Reference Manual.
Particle Decay Class¶
Geant4 provides a G4Decay
class for both at rest
and in
flight
particle decays. G4Decay
can be applied to all particles
except:
- massless particles, i.e.,
G4ParticleDefinition::thePDGMass <= 0
- particles with “negative” life time, i.e.,
G4ParticleDefinition::thePDGLifeTime < 0
- shortlived particles, i.e.,
G4ParticleDefinition::fShortLivedFlag = True
Decay for some particles may be switched on or off by using
G4ParticleDefinition::SetPDGStable()
as well as
ActivateProcess()
and InActivateProcess()
methods of
G4ProcessManager
.
G4Decay
proposes the step length (or step time for AtRest
)
according to the lifetime of the particle unless
PreAssignedDecayProperTime
is defined in G4DynamicParticle
.
The G4Decay
class itself does not define decay modes of the
particle. Geant4 provides two ways of doing this:
using
G4DecayChannel
inG4DecayTable
, andusing
thePreAssignedDecayProducts
ofG4DynamicParticle
The G4Decay
class calculates the PhysicalInteractionLength
and
boosts decay products created by G4VDecayChannel
or event
generators. See below for information on the determination of the decay
modes.
An object of G4Decay
can be shared by particles. Registration of the
decay process to particles in the ConstructPhysics
method of
PhysicsList (see How to Specify Physics Processes) is
shown in Listing 60.
#include "G4Decay.hh"
void MyPhysicsList::ConstructGeneral()
{
// Add Decay Process
G4Decay* theDecayProcess = new G4Decay();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (theDecayProcess->IsApplicable(*particle)) {
pmanager ->AddProcess(theDecayProcess);
// set ordering for PostStepDoIt and AtRestDoIt
pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
}
}
}
Decay Table¶
Each particle has its G4DecayTable
, which stores information on the
decay modes of the particle. Each decay mode, with its branching ratio,
corresponds to an object of various “decay channel” classes derived
from G4VDecayChannel
. Default decay modes are created in the
constructors of particle classes. For example, the decay table of the
neutral pion has G4PhaseSpaceDecayChannel
and
G4DalitzDecayChannel
as follows:
// create a decay channel
G4VDecayChannel* mode;
// pi0 -> gamma + gamma
mode = new G4PhaseSpaceDecayChannel("pi0",0.988,2,"gamma","gamma");
table->Insert(mode);
// pi0 -> gamma + e+ + e-
mode = new G4DalitzDecayChannel("pi0",0.012,"e-","e+");
table->Insert(mode);
Decay modes and branching ratios defined in Geant4 are listed in Definition of a particle.
Branching ratios and life time can be set in tracking time.
// set lifetime
G4Neutron::Neutron()->SetPDGLifeTime(885.7*second);
// allow neutron decay
G4Neutron::Neutron()->SetPDGStable(false);
Branching ratios and life time can be modified by using user commands, also.
Example: Set 100% br for dalitz decay of pi0
Idle> /particle/select pi0
Idle> /particle/property/decay/select 0
Idle> /particle/property/decay/br 0
Idle> /particle/property/decay/select 1
Idle> /particle/property/decay/br 1
Idle> /particle/property/decay/dump
G4DecayTable: pi0
0: BR: 0 [Phase Space] : gamma gamma
1: BR: 1 [Dalitz Decay] : gamma e- e+
Pre-assigned Decay Modes by Event Generators¶
Decays of heavy flavor particles such as B mesons are very complex, with
many varieties of decay modes and decay mechanisms. There are many
models for heavy particle decay provided by various event generators and
it is impossible to define all the decay modes of heavy particles by
using G4VDecayChannel
. In other words, decays of heavy particles
cannot be defined by the Geant4 decay process, but should be defined by
event generators or other external packages. Geant4 provides two ways to
do this: pre-assigned decay mode
and external decayer
.
In the latter approach, the class G4VExtDecayer
is used for the
interface to an external package which defines decay modes for a
particle. If an instance of G4VExtDecayer
is attached to
G4Decay
, daughter particles will be generated by the external decay
handler.
In the former case, decays of heavy particles are simulated by an event
generator and the primary event contains the decay information.
G4VPrimaryGenerator
automatically attaches any daughter particles to
the parent particle as the PreAssignedDecayProducts member of
G4DynamicParticle
. G4Decay
adopts these pre-assigned daughter
particles instead of asking G4VDecayChannel
to generate decay
products.
In addition, the user may assign a pre-assigned
decay time for a
specific track in its rest frame (i.e. decay time is defined in the
proper time) by using the G4PrimaryParticle::SetProperTime()
method.
G4VPrimaryGenerator
sets the PreAssignedDecayProperTime member of
G4DynamicParticle
. G4Decay
uses this decay time instead of the
life time of the particle type.
Gamma-nuclear and Lepto-nuclear Processes¶
Gamma-nuclear and lepto-nuclear reactions are handled in Geant4 as hybrid processes which typically require both electromagnetic and hadronic models for their implementation. While neutrino-induced reactions are not currently provided, the Geant4 hadronic framework is general enough to include their future implementation as a hybrid of weak and hadronic models.
The general scheme followed is to factor the full interaction into an electromagnetic (or weak) vertex, in which a virtual particle is generated, and a hadronic vertex in which the virtual particle interacts with a target nucleus. In most cases the hadronic vertex is implemented by an existing Geant4 model which handles the intra-nuclear propagation.
The cross sections for these processes are parameterizations, either directly of data or of theoretical distributions determined from the integration of lepton-nucleon cross sections double differential in energy loss and momentum transfer.
For the most part gammas can be treated as hadrons and in fact they
interact that way with the nucleus when the Bertini-style cascade
G4CascadeInterface
and QGSP models are used. These models may be
assigned to G4PhotoNuclearProcess
as shown in the following partial
code:
G4TheoFSGenerator* theHEModel = new G4TheoFSGenerator;
G4QGSModel* theStringModel = new G4QGSModel<G4GammaParticipants>;
G4ExcitedStringDecay* theStringDecay =
new G4ExcitedStringDecay(theFragmentation=new G4QGSMFragmentation);
theStringModel->SetFragmentationModel(theStringDecay);
theHEModel->SetHighEnergyGenerator(theStringModel);
theHEModel->SetTransport(new G4GeneratorPrecompoundInterface);
theHEModel->SetMinEnergy(8*GeV);
G4CascadeInterface* theLEModel = new G4CascadeInterface;
theLEModel->SetMaxEnergy(10*GeV);
G4PhotoNuclearProcess* thePhotoNuclearProcess = new G4PhotoNuclearProcess;
thePhotoNuclearProcess->RegisterMe(theLEModel);
thePhotoNuclearProcess->RegisterMe(theHEModel);
G4ProcessManager* procMan = G4Gamma::Gamma()->GetProcessManager();
procMan->AddDiscreteProcess(thePhotoNuclearProcess);
Electro-nuclear reactions in Geant4 are handled by the classes
G4ElectronNuclearProcess
and G4PositronNuclearProcess
, which are
both implmented by G4ElectroVDNuclearModel
. This model consists of
three sub-models: code which generates the virtual photon from the
lepton-nucleus vertex, the Bertini-style cascade to handle the low and
medium energy photons, and the FTFP model to handle the high energy
photons.
Muon-nuclear reactions are handled similarly. The process
G4MuonNuclearProcess
can be assigned the G4MuonVDNuclearModel
which in turn is implemented by three sub-models: virtual gamma
generation code, Bertini-style cascade and the FTFP model.
Optical Photon Processes¶
A photon is considered to be optical when its wavelength is much greater than the typical atomic spacing. In Geant4 optical photons are treated as a class of particle distinct from their higher energy gamma cousins. This implementation allows the wave-like properties of electromagnetic radiation to be incorporated into the optical photon process. Because this theoretical description breaks down at higher energies, there is no smooth transition as a function of energy between the optical photon and gamma particle classes.
For the simulation of optical photons to work correctly in Geant4, they must be imputed a linear polarization. This is unlike most other particles in Geant4 but is automatically and correctly done for optical photons that are generated as secondaries by existing processes in Geant4. If the user wishes to start optical photons as primary particles, they must set the linear polarization using particle gun methods, the General Particle Source, or their PrimaryGeneratorAction. For an unpolarized source, the linear polarization should be sampled randomly for each new primary photon.
The Geant4 catalogue of processes at optical wavelengths includes refraction and reflection at medium boundaries, bulk absorption, Mie and Rayleigh scattering. Processes which produce optical photons include the Cerenkov effect and scintillation. Optical photons are generated in Geant4 without energy conservation and their energy must therefore not be tallied as part of the energy balance of an event.
The optical properties of the medium which are key to the implementation
of these types of processes are stored as entries in a
G4MaterialPropertiesTable
which is linked to the G4Material
in
question. These properties may be independent of energy (denoted “Constant”
or “Const”) or they may be expressed as
a function of the photon’s energy. This table is a private data member
of the G4Material
class. The G4MaterialPropertiesTable
is
implemented as a hash directory, in which each entry consists of a
value and a key. The key is used to quickly and efficiently retrieve
the corresponding value. All values in the dictionary are either
instantiations of G4double
or the class
G4MaterialPropertyVector
, and all keys are of type G4String
.
A G4MaterialPropertyVector
is a typedef of
G4PhysicsOrderedFreeVector. The entries are a pair of numbers, which in
the case of an optical property, are the photon energy and corresponding
property value. It is possible for the user to add as many material
(optical) properties to the material as he wishes using the methods
supplied by the G4MaterialPropertiesTable
class. An example of this
is shown in Listing 61. In this
example the interpolation of the G4MaterialPropertyVector is to be done
by a spline fit. The default is a linear interpolation.
const G4int NUMENTRIES = 3;
G4double ppckov[NUMENTRIES] = {2.034*eV, 3.*eV, 4.136*eV};
G4double rindex[NUMENTRIES] = {1.3435, 1.351, 1.3608};
G4double absorption[NUMENTRIES] = {344.8*cm, 850.*cm, 1450.0*cm};
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("SCINTILLATIONYIELD",100./MeV);
MPT->AddProperty("RINDEX",ppckov,rindex,NUMENTRIES}->SetSpline(true);
MPT->AddProperty("ABSLENGTH",ppckov,absorption,NUMENTRIES}->SetSpline(true);
scintillator -> SetMaterialPropertiesTable(MPT);
Note
It is possible to define new material properties. For this reason there is no protection against mis-spelling a property name. The user needs to carefully check the property name.
G4OpticalPhysics
constructor¶
The most straightforward way of using optical physics is to use the G4OpticalPhysics
constructor in main(), as in the extended optical examples. This automatically includes
all the optical physics processes. An example of using G4OpticalPhysics
is shown in
Listing 62.
#include "G4OpticalPhysics.hh"
...
G4VModularPhysicsList* physicsList = new FTFP_BERT; // for example
G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics();
physicsList->RegisterPhysics(opticalPhysics);
runManager->SetUserInitialization(physicsList);
There is an associated messenger class to set many of the parameters. See the sections for each process for details. Commands are:
/process/optical/verbose
/process/optical/cerenkov/setMaxPhotons
/process/optical/cerenkov/setMaxBetaChange
/process/optical/cerenkov/setStackPhotons
/process/optical/cerenkov/setTrackSecondariesFirst
/process/optical/cerenkov/verbose
/process/optical/scintillation/setYieldFactor
/process/optical/scintillation/setExcitationRatio
/process/optical/scintillation/setByParticleType
/process/optical/scintillation/setTrackInfo
/process/optical/scintillation/setFiniteRiseTime
/process/optical/scintillation/setStackPhotons
/process/optical/scintillation/setTrackSecondariesFirst
/process/optical/scintillation/verbose
/process/optical/wls/setTimeProfile
/process/optical/wls/verbose
/process/optical/boundary/setInvokeSD
/process/optical/boundary/verbose
/process/optical/absorption/verbose
/process/optical/rayleigh/verbose
/process/optical/mie/verbose
Generation of Photons in processes/electromagnetic/xrays
- Cerenkov Effect¶
The radiation of Cerenkov light occurs when a charged particle moves through a dispersive medium faster than the group velocity of light in that medium. Photons are emitted on the surface of a cone, whose opening angle with respect to the particle’s instantaneous direction decreases as the particle slows down. At the same time, the frequency of the photons emitted increases, and the number produced decreases. When the particle velocity drops below the local speed of light, the radiation ceases and the emission cone angle collapses to zero. The photons produced by this process have an inherent polarization perpendicular to the cone’s surface at production.
The flux, spectrum, polarization and emission of Cerenkov radiation in
the AlongStepDoIt
method of the class G4Cerenkov
follow
well-known formulae, with two inherent computational limitations. The
first arises from step-wise simulation, and the second comes from the
requirement that numerical integration calculate the average number of
Cerenkov photons per step. The process makes use of a G4PhysicsTable
which contains incremental integrals to expedite this calculation. The
Cerenkov process in Geant4 is limited to normally dispersive media,
i.e., \(dn(E)/dE \ge 0\).
The time and position of Cerenkov photon emission are calculated from
quantities known at the beginning of a charged particle’s step. The step
is assumed to be rectilinear even in the presence of a magnetic field.
The user may limit the step size by specifying a maximum (average)
number of Cerenkov photons created during the step, using the
SetMaxNumPhotonsPerStep(const G4int
NumPhotons)
method. The actual number generated will necessarily be
different due to the Poissonian nature of the production. In the present
implementation, the production density of photons is distributed evenly
along the particle’s track segment, even if the particle has slowed
significantly during the step. The step can also be limited with the
SetMaxBetaChangePerStep
method, where the argument is the allowed
change in percent.
The large number of optical photons that can be produced
(about 300/cm in water) can fill the available memory. Geant4 provides
the public method SetTrackSecondariesFirst
that suspends the primary
particle until all its progeny have been tracked. This is on by default
int the G4OpticalPhysics
class. An
example of the registration of the Cerenkov process is given in
Listing 63.
Note
A user may avoid this by registering the G4OpticalPhysics
constructor.
See G4OpticalPhysics constructor.
#include "G4Cerenkov.hh"
void ExptPhysicsList::ConstructOp(){
G4Cerenkov* theCerenkovProcess = new G4Cerenkov("Cerenkov");
G4int MaxNumPhotons = 300;
theCerenkovProcess->SetTrackSecondariesFirst(true);
theCerenkovProcess->SetMaxBetaChangePerStep(10.0);
theCerenkovProcess->SetMaxNumPhotonsPerStep(MaxNumPhotons);
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (theCerenkovProcess->IsApplicable(*particle)) {
pmanager->AddProcess(theCerenkovProcess);
pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
}
}
}
Generation of Photons in processes/electromagnetic/xrays
- Scintillation¶
Every scintillating material has a characteristic light yield,
SCINTILLATIONYIELD
, and an intrinsic resolution,
RESOLUTIONSCALE
, which generally broadens the statistical
distribution of generated photons. A wider intrinsic resolution is due
to impurities which are typical for doped crystals like NaI(Tl) and
CsI(Tl). On the other hand, the intrinsic resolution can also be
narrower when the Fano factor plays a role. The actual number of emitted
photons during a step fluctuates around the mean number of photons with
a width given by ResolutionScale*sqrt(MeanNumberOfPhotons)
. The
average light yield, MeanNumberOfPhotons
, has a linear dependence on
the local energy deposition, but it may be different for minimum
ionizing and non-minimum ionizing particles.
A scintillator is also characterized by its photon emission spectrum and
by the exponential decay of its time spectrum. In Geant4 the
scintillator can have a fast and a slow component. (Only one component is
required, and if there is only one component it can be either fast or slow.)
The relative strength
of the fast component as a fraction of total scintillation yield is
given by the YIELDRATIO
. Scintillation may be simulated by
specifying these empirical parameters for each material. It is
sufficient to specify in the user’s DetectorConstruction
class a
relative spectral distribution as a function of photon energy for the
scintillating material. An example of this is shown in
Listing 64. For each time constant, the spectrum must be
specified by FASTCOMPONENT
or SLOWCOMPONENT
. The time constant must
be specified by FASTTIMECONSTANT
or SLOWTIMECONSTANT
. Rise times may
optionally be specified with FASTSCINTILLATIONRISETIME
and
SLOWSCINTILLATIONRISETIME
.
const G4int NUMENTRIES = 9;
G4double Scnt_PP[NUMENTRIES] = { 6.6*eV, 6.7*eV, 6.8*eV, 6.9*eV,
7.0*eV, 7.1*eV, 7.2*eV, 7.3*eV, 7.4*eV };
G4double Scnt_FAST[NUMENTRIES] = { 0.000134, 0.004432, 0.053991, 0.241971,
0.398942, 0.000134, 0.004432, 0.053991,
0.241971 };
G4double Scnt_SLOW[NUMENTRIES] = { 0.000010, 0.000020, 0.000030, 0.004000,
0.008000, 0.005000, 0.020000, 0.001000,
0.000010 };
G4Material* Scnt;
G4MaterialPropertiesTable* Scnt_MPT = new G4MaterialPropertiesTable();
Scnt_MPT->AddProperty("FASTCOMPONENT", Scnt_PP, Scnt_FAST, NUMENTRIES);
Scnt_MPT->AddProperty("SLOWCOMPONENT", Scnt_PP, Scnt_SLOW, NUMENTRIES);
Scnt_MPT->AddConstProperty("SCINTILLATIONYIELD", 5000./MeV);
Scnt_MPT->AddConstProperty("RESOLUTIONSCALE", 2.0);
Scnt_MPT->AddConstProperty("FASTTIMECONSTANT", 1.*ns);
Scnt_MPT->AddConstProperty("SLOWTIMECONSTANT", 10.*ns);
Scnt_MPT->AddConstProperty("YIELDRATIO", 0.8);
Scnt->SetMaterialPropertiesTable(Scnt_MPT);
In cases where the scintillation yield of a scintillator depends on the
particle type, different scintillation processes may be defined for
them. How this yield scales to the one specified for the material is
expressed with the ScintillationYieldFactor
in the user’s
PhysicsList
as shown in the
Listing 65 below. In those cases where the
fast to slow excitation ratio changes with particle type, the method
SetScintillationExcitationRatio
can be called for each scintillation
process (see the advanced underground_physics example). This overwrites
the YieldRatio
obtained from the G4MaterialPropertiesTable
.
Note
A user may avoid this by registering the G4OpticalPhysics
constructor.
See G4OpticalPhysics constructor.
G4Scintillation* theMuonScintProcess = new G4Scintillation("Scintillation");
theMuonScintProcess->SetTrackSecondariesFirst(true);
theMuonScintProcess->SetScintillationYieldFactor(0.8);
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if (theMuonScintProcess->IsApplicable(*particle)) {
if (particleName == "mu+") {
pmanager->AddProcess(theMuonScintProcess);
pmanager->SetProcessOrderingToLast(theMuonScintProcess, idxAtRest);
pmanager->SetProcessOrderingToLast(theMuonScintProcess, idxPostStep);
}
}
}
A Gaussian-distributed number of photons is generated according to the
energy lost during the step. A resolution scale of 1.0 produces a
statistical fluctuation around the average yield set with
AddConstProperty("SCINTILLATIONYIELD")
, while values > 1 broaden the
fluctuation. A value of zero produces no fluctuation. Each photon’s
frequency is sampled from the empirical spectrum. The photons originate
evenly along the track segment and are emitted uniformly into \(4\pi\) with a
random linear polarization and at times characteristic for the
scintillation component. The SCINTILLATIONYIELD
specified here is
the differential yield.
When there are multiple scintillators in the simulation and/or when the scintillation yield is a non-linear function of the energy deposited, the user can also define an array of total scintillation light yields as a function of the energy deposited and particle type. The available particles are protons, electrons, deuterons, tritons, alphas, and carbon ions. These are the particles known to significantly effect the scintillation light yield, of for example, BC501A (NE213/EJ301) liquid organic scintillator and BC420 plastic scintillator as function of energy deposited.
The method works as follows:
In the user’s physics lists, the user must set a G4bool flag that allows scintillation light emission to depend on the energy deposited by particle type:
theScintProcess->SetScintillationByParticleType(true);
The
G4OpticalPhysics
constructor includes a macro command for this:/process/optical/scintillation/setByParticleType true
The user must also specify and add, via the AddProperty method of the MPT, the scintillation light yield as function of incident particle energy with new keywords,
PROTONSCINTILLATIONYIELD
,DEUTERONSCINTILLATIONYIELD
,TRITONSCINTILLATIONYIELD
,ALPHASCINTILLATIONYIELD
,IONSCINTILLATIONYIELD
, andELECTRONSCINTILLATIONYIELD
, and pairs of protonEnergy and scintLightYield. These particle-specific yields are the total yields (not differential) as a function of energy.
Generation of Photons in processes/optical
- Wavelength Shifting¶
Wavelength Shifting (WLS) fibers are used in many high-energy particle physics experiments. They absorb light at one wavelength and re-emit light at a different wavelength and are used for several reasons. For one, they tend to decrease the self-absorption of the detector so that as much light reaches the PMTs as possible. WLS fibers are also used to match the emission spectrum of the detector with the input spectrum of the PMT.
A WLS material is characterized by its photon absorption and photon
emission spectrum and by a possible time delay between the absorption
and re-emission of the photon. Wavelength Shifting may be simulated by
specifying these empirical parameters for each WLS material in the
simulation. It is sufficient to specify in the user’s
DetectorConstruction
class a relative spectral distribution as a
function of photon energy for the WLS material. WLSABSLENGTH
is the
absorption length of the material as a function of the photon’s energy.
WLSCOMPONENT
is the relative emission spectrum of the material as a
function of the photon’s energy, and WLSTIMECONSTANT
accounts for any
time delay which may occur between absorption and re-emission of the
photon. An example is shown in the
Listing 66.
const G4int nEntries = 9;
G4double PhotonEnergy[nEntries] = { 6.6*eV, 6.7*eV, 6.8*eV, 6.9*eV,
7.0*eV, 7.1*eV, 7.2*eV, 7.3*eV, 7.4*eV};
G4double RIndexFiber[nEntries] =
{1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60, 1.60};
G4double AbsFiber[nEntries] =
{0.1*mm,0.2*mm,0.3*mm,0.4*cm,1.0*cm,10.0*cm,1.0*m,10.0*m,10.0*m};
G4double EmissionFiber[nEntries] =
{0.0, 0.0, 0.0, 0.1, 0.5, 1.0, 5.0, 10.0, 10.0};
G4Material* WLSFiber;
G4MaterialPropertiesTable* MPTFiber = new G4MaterialPropertiesTable();
MPTFiber->AddProperty("RINDEX",PhotonEnergy,RIndexFiber,nEntries);
MPTFiber->AddProperty("WLSABSLENGTH",PhotonEnergy,AbsFiber,nEntries);
MPTFiber->AddProperty("WLSCOMPONENT",PhotonEnergy,EmissionFiber,nEntries);
MPTFiber->AddConstProperty("WLSTIMECONSTANT", 0.5*ns);
WLSFiber->SetMaterialPropertiesTable(MPTFiber);
The process is defined in the PhysicsList in the usual way. The process
class name is G4OpWLS
. It should be instantiated with
theWLSProcess = new G4OpWLS("OpWLS")
and attached to the process
manager of the optical
photon as a DiscreteProcess. The way the WLSTIMECONSTANT
is used depends
on the time profile method chosen by the user. If in the PhysicsList
theWLSProcess->UseTimeGenerator("exponential")
option is set, the time
delay between absorption and re-emission of the photon is sampled from
an exponential distribution, with the decay term equal to
WLSTIMECONSTANT
. If, on the other hand,
theWLSProcess->UseTimeGenerator("delta")
is chosen, the time delay is a
delta function and equal to WLSTIMECONSTANT
. The default is “delta” in
case the G4OpWLS::UseTimeGenerator(const G4String name)
method is not
used.
Tracking of Photons in processes/optical
¶
Absorption¶
The implementation of optical photon bulk absorption,
G4OpAbsorption
, is trivial in that the process merely kills the
particle. The procedure requires the user to fill the relevant
G4MaterialPropertiesTable
with empirical data for the absorption
length, using ABSLENGTH
as the property key in the public method
AddProperty
. The absorption length is the average distance traveled
by a photon before being absorbed by the medium; i.e., it is the mean
free path returned by the GetMeanFreePath
method.
Rayleigh Scattering¶
The differential cross section in Rayleigh scattering, \(d\sigma/d\Omega\), is
proportional to \(1+\cos^2 \theta\), where \(\theta\) is the polar of
the new polarization vector with respect to the old polarization vector.
The G4OpRayleigh
scattering process samples this angle accordingly
and then calculates the scattered photon’s new direction by requiring
that it be perpendicular to the photon’s new polarization in such a way
that the final direction, initial and final polarizations are all in one
plane. This process thus depends on the particle’s polarization (spin).
The photon’s polarization is a data member of the G4DynamicParticle
class.
A photon which is not assigned a polarization at production, either via
the SetPolarization
method of the G4PrimaryParticle
class, or
indirectly with the SetParticlePolarization
method of the
G4ParticleGun
class, may not be Rayleigh scattered. Optical photons
produced by the G4Cerenkov
process have inherently a polarization
perpendicular to the cone’s surface at production. Scintillation photons
have a random linear polarization perpendicular to their direction.
The process requires a G4MaterialPropertiesTable
to be filled by the
user with Rayleigh scattering length data. The Rayleigh scattering
attenuation length is the average distance traveled by a photon before
it is Rayleigh scattered in the medium and it is the distance returned
by the GetMeanFreePath
method. The G4OpRayleigh
class provides a
RayleighAttenuationLengthGenerator
method which calculates the
attenuation coefficient of a medium following the Einstein-Smoluchowski
formula whose derivation requires the use of statistical mechanics,
includes temperature, and depends on the isothermal compressibility of
the medium. This generator is convenient when the Rayleigh attenuation
length is not known from measurement but may be calculated from first
principles using the above material constants. For a medium named
Water and no Rayleigh scattering attenuation length specified by the
user, the program automatically calls the
RayleighAttenuationLengthGenerator
which calculates it for 10ºC
liquid water.
Mie Scattering¶
Mie Scattering (or Mie solution) is an analytical solution of Maxwell’s equations for scattering of optical photons by spherical particles. It is significant only when the radius of the scattering object is of order of the wave length. The analytical expressions for Mie Scattering are very complicated since they are a series sum of Bessel functions. One common approximation made is call Henyey-Greenstein (HG). The implementation in Geant4 follows the HG approximation (for details see the Physics Reference Manual. ) and the treatment of polarization and momentum are similar to that of Rayleigh scattering. We require the final polarization direction to be perpendicular to the momentum direction. We also require the final momentum, initial polarization and final polarization to be in the same plane.
The process requires a G4MaterialPropertiesTable
to be filled by the
user with Mie scattering length data (entered with the name MIEHG
)
analogous to Rayleigh scattering. The Mie scattering attenuation length
is the average distance traveled by a photon before it is Mie scattered
in the medium and it is the distance returned by the GetMeanFreePath
method. In practice, the user not only needs to provide the attenuation
length of Mie scattering, but also needs to provide the constant
parameters of the approximation: g_f, g_b, and r_f. (with
AddConstProperty
and with the names: MIEHG_FORWARD
, MIEHG_BACKWARD
,
and MIEHG_FORWARD_RATIO
, respectively; see extended example optical/OpNovice.)
Boundary Process¶
Reference: E. Hecht and A. Zajac, Optics [Hecht1974]
For the simple case of a perfectly smooth interface between two
dielectric materials, all the user needs to provide are the refractive
indices of the two materials stored in their respective
G4MaterialPropertiesTable
. In all other cases, the optical boundary
process design relies on the concept of surfaces. The information is
split into two classes. One class in the material category keeps
information about the physical properties of the surface itself, and a
second class in the geometry category holds pointers to the relevant
physical and logical volumes involved and has an association to the
physical class. Surface objects of the second type are stored in a
related table and can be retrieved by either specifying the two ordered
pairs of physical volumes touching at the surface, or by the logical
volume entirely surrounded by this surface. The former is called a
border surface while the latter is referred to as the skin surface.
This second type of surface is useful in situations where a volume is
coded with a reflector and is placed into many different mother volumes.
A limitation is that the skin surface can only have one and the same
optical property for all of the enclosed volume’s sides. The border
surface is an ordered pair of physical volumes, so in principle, the
user can choose different optical properties for photons arriving from
the reverse side of the same interface. For the optical boundary process
to use a border surface, the two volumes must have been positioned with
G4PVPlacement
. The ordered combination can exist at many places in
the simulation. When the surface concept is not needed, and a perfectly
smooth surface exists between two dielectric materials, the only relevant
property is the index of refraction, a quantity stored with the
material, and no restriction exists on how the volumes were positioned.
When an optical photon arrives at a boundary it is absorbed if the medium of the volume being left behind has no index of refraction defined. A photon is also absorbed in case of a dielectric-dielectric polished or ground surface when the medium about to be entered has no index of refraction. It is absorbed for backpainted surfaces when the surface has no index of refraction. If the geometry boundary has a border surface this surface takes precedence, otherwise the program checks for skin surfaces. The skin surface of the daughter volume is taken if a daughter volume is entered else the program checks for a skin surface of the current volume. When the optical photon leaves a volume without entering a daughter volume the skin surface of the current volume takes precedence over that of the volume about to be entered.
The physical surface object also specifies which model the boundary
process should use to simulate interactions with that surface. In
addition, the physical surface can have a material property table all
its own. The usage of this table allows all specular constants to be
wavelength dependent. In case the surface is painted or wrapped (but not
a cladding), the table may include the thin layer’s index of refraction.
This allows the simulation of boundary effects at the intersection
between the medium and the surface layer, as well as the Lambertian
reflection at the far side of the thin layer. This occurs within the
process itself and does not invoke the G4Navigator
. Combinations of
surface finish properties, such as polished or ground and front
painted or back painted, enumerate the different situations which can
be simulated.
When a photon arrives at a medium boundary its behavior depends on the nature of the two materials that join at that boundary. Medium boundaries may be formed between two dielectric materials or a dielectric and a metal. In the case of two dielectric materials, the photon can undergo total internal reflection, refraction or reflection, depending on the photon’s wavelength, angle of incidence, and the refractive indices on both sides of the boundary. Furthermore, reflection and transmission probabilities are sensitive to the state of linear polarization. In the case of an interface between a dielectric and a metal, the photon can be absorbed by the metal or reflected back into the dielectric. If the photon is absorbed it can be detected according to the photoelectron efficiency of the metal.
As expressed in Maxwell’s equations, Fresnel reflection and refraction are intertwined through their relative probabilities of occurrence. Therefore neither of these processes, nor total internal reflection, are viewed as individual processes deserving separate class implementation. Nonetheless, an attempt was made to adhere to the abstraction of having independent processes by splitting the code into different methods where practicable.
One implementation of the G4OpBoundaryProcess
class employs the
UNIFIED model
[Levin1996] of the DETECT program [Knoll1988]. It applies to
dielectric-dielectric interfaces and tries
to provide a realistic simulation, which deals with all aspects of
surface finish and reflector coating. The surface may be assumed as
smooth and covered with a metallized coating representing a specular
reflector with given reflection coefficient, or painted with a diffuse
reflecting material where Lambertian reflection occurs. The surfaces may
or may not be in optical contact with another component and most
importantly, one may consider a surface to be made up of micro-facets
with normal vectors that follow given distributions around the nominal
normal for the volume at the impact point. For very rough surfaces, it
is possible for the photon to inversely aim at the same surface again
after reflection of refraction and so multiple interactions with the
boundary are possible within the process itself and without the need for
relocation by G4Navigator
.
The UNIFIED model (Fig. 17) provides for a range of different reflection mechanisms. The specular lobe constant represents the reflection probability about the normal of a micro facet. The specular spike constant, in turn, illustrates the probability of reflection about the average surface normal. The diffuse lobe constant is for the probability of internal Lambertian reflection, and finally the back-scatter spike constant is for the case of several reflections within a deep groove with the ultimate result of exact back-scattering. The four probabilities must add up to one, with the diffuse lobe constant being implicit. The reader may consult the reference for a thorough description of the model.
G4VPhysicalVolume* volume1;
G4VPhysicalVolume* volume2;
G4OpticalSurface* OpSurface = new G4OpticalSurface("name");
G4LogicalBorderSurface* Surface = new
G4LogicalBorderSurface("name",volume1,volume2,OpSurface);
G4double sigma_alpha = 0.1;
OpSurface->SetType(dielectric_dielectric);
OpSurface->SetModel(unified);
OpSurface->SetFinish(groundbackpainted);
OpSurface->SetSigmaAlpha(sigma_alpha);
const G4int NUM = 2;
G4double pp[NUM] = {2.038*eV, 4.144*eV};
G4double specularlobe[NUM] = {0.3, 0.3};
G4double specularspike[NUM] = {0.2, 0.2};
G4double backscatter[NUM] = {0.1, 0.1};
G4double rindex[NUM] = {1.35, 1.40};
G4double reflectivity[NUM] = {0.3, 0.5};
G4double efficiency[NUM] = {0.8, 0.1};
G4MaterialPropertiesTable* SMPT = new G4MaterialPropertiesTable();
SMPT->AddProperty("RINDEX",pp,rindex,NUM);
SMPT->AddProperty("SPECULARLOBECONSTANT",pp,specularlobe,NUM);
SMPT->AddProperty("SPECULARSPIKECONSTANT",pp,specularspike,NUM);
SMPT->AddProperty("BACKSCATTERCONSTANT",pp,backscatter,NUM);
SMPT->AddProperty("REFLECTIVITY",pp,reflectivity,NUM);
SMPT->AddProperty("EFFICIENCY",pp,efficiency,NUM);
OpSurface->SetMaterialPropertiesTable(SMPT);
The original GEANT3.21 implementation of this process is also available via the GLISUR methods flag, as shown in Listing 68.
G4LogicalVolume* volume_log;
G4OpticalSurface* OpSurface = new G4OpticalSurface("name");
G4LogicalSkinSurface* Surface = new
G4LogicalSkinSurface("name",volume_log,OpSurface);
OpSurface->SetType(dielectric_metal);
OpSurface->SetFinish(ground);
OpSurface->SetModel(glisur);
G4double polish = 0.8;
G4MaterialPropertiesTable* OpSurfaceProperty = new G4MaterialPropertiesTable();
OpSurfaceProperty->AddProperty("REFLECTIVITY",pp,reflectivity,NUM);
OpSurfaceProperty->AddProperty("EFFICIENCY",pp,efficiency,NUM);
OpSurface->SetMaterialPropertiesTable(OpSurfaceProperty);
The reflectivity off a metal surface can also be calculated by way of a
complex index of refraction. Instead of storing the REFLECTIVITY
directly, the user stores the real part (REALRINDEX
) and the imaginary
part (IMAGINARYRINDEX
) as a function of photon energy separately in the
G4MaterialPropertyTable. Geant4 then calculates the
reflectivity
depending on the incident angle, photon energy, degree of TE and TM
polarization, and this complex refractive index.
The program defaults to the GLISUR model and polished surface finish
when no specific model and surface finish is specified by the user. In
the case of a dielectric-metal interface, or when the GLISUR model is
specified, the only surface finish options available are polished or
ground. For dielectric-metal surfaces, the G4OpBoundaryProcess
also defaults to unit reflectivity and zero detection efficiency. In
cases where the user specifies the UNIFIED model
(Fig. 17), but does not otherwise
specify the model reflection probability constants, the default becomes
Lambertian reflection.
Janecek and Moses [Janecek2010]
built an instrument for measuring the angular reflectivity distribution
inside of BGO crystals with common surface treatments and reflectors
applied. These results have been incorporated into the Geant4 code. A
third class of reflection type besides dielectric_metal and
dielectric_dielectric is added: dielectric_LUT. The distributions have
been converted to 21 look-up-tables (LUT); so far for 1 scintillator
material (BGO) x 3 surface treatments x 7 reflector materials. The
modified code allows the user to specify the surface treatment
(rough-cut, chemically etched, or mechanically polished), the attached
reflector (Lumirror, Teflon, ESR film, Tyvek, or TiO2 paint), and the
bonding type (air-coupled or glued). The glue used is MeltMount, and the
ESR film used is VM2000. Each LUT consists of measured angular
distributions with 4º by 5º resolution in theta and phi, respectively,
for incidence angles from 0º to 90º, in 1º-steps. The code might
in the future be updated by adding more LUTs, for instance, for other
scintillating materials (such as LSO or NaI). To use these LUT the user
has to download them from Geant4 Software
Download
and set an environment variable, G4REALSURFACEDATA
, to the directory
of geant4/data/RealSurface2.1.1
.
The enumeration G4OpticalSurfaceFinish includes:
polishedlumirrorair, // mechanically polished surface, with lumirror
polishedlumirrorglue, // mechanically polished surface, with lumirror & meltmount
polishedteflonair, // mechanically polished surface, with teflon
polishedtioair, // mechanically polished surface, with tio paint
polishedtyvekair, // mechanically polished surface, with tyvek
polishedvm2000air, // mechanically polished surface, with esr film
polishedvm2000glue, // mechanically polished surface, with esr film & meltmount
etchedlumirrorair, // chemically etched surface, with lumirror
etchedlumirrorglue, // chemically etched surface, with lumirror & meltmount
etchedteflonair, // chemically etched surface, with teflon
etchedtioair, // chemically etched surface, with tio paint
etchedtyvekair, // chemically etched surface, with tyvek
etchedvm2000air, // chemically etched surface, with esr film
etchedvm2000glue, // chemically etched surface, with esr film & meltmount
groundlumirrorair, // rough-cut surface, with lumirror
groundlumirrorglue, // rough-cut surface, with lumirror & meltmount
groundteflonair, // rough-cut surface, with teflon
groundtioair, // rough-cut surface, with tio paint
groundtyvekair, // rough-cut surface, with tyvek
groundvm2000air, // rough-cut surface, with esr film
groundvm2000glue // rough-cut surface, with esr film & meltmount
To use a look-up-table, all the user needs to specify for an
G4OpticalSurface
is: SetType(dielectric_LUT), SetModel(LUT)
and
for example, SetFinish(polishedtyvekair)
.
A newly implemented model for optical transport is called the LUT Davis model [RoncaliCherry2013], [Stockhoff2017], [Roncali2017]. The model is based on measured surface data and allows the user to choose from a list of available surface finishes. Provided are a rough and a polished L(Y)SO surface that can be used without reflector, or in combination with a specular reflector (e.g. ESR) or a Lambertian reflector (e.g. Teflon). The specular reflector can be coupled to the crystal with air or optical grease. Teflon tape is wrapped around the crystal with 4 layers.
The user can extend the list of finishes with custom measured surface data.
Bare |
Teflon |
ESR + Air |
ESR + Optical Grease |
|
---|---|---|---|---|
Rough |
Rough_LUT |
RoughTeflon_LUT |
RoughESR_LUT |
RoughESRGrease_LUT |
Polished |
Polished_LUT |
PolishedTeflon_LUT |
PolishedESR_LUT |
PolishedESRGrease_LUT |
In the LUT database, typical roughness parameters obtained from the measurements are provided to characterize the type of surface modelled:
with \(R_a\) = average roughness; \(\sigma \) = rms roughness, \(R_{pv}\) = peak-to-valley ratio.
The detector surface, called Detector_LUT
, defines a polished surface
coupled to a photodetector with optical grease or a glass interface
(similar index of refraction 1.5). To use Detector_LUT
, the surface property
EFFICIENCY
must be greater than 0. Any surface can be used as a detector
surface when the EFFICIENCY
is set to 1.
To enable the LUT Davis Model, the user needs to specify for a
G4OpticalSurface: SetType(dielectric_LUTDAVIS), SetModel(DAVIS)
and also,
for example, SetFinish(Rough_LUT)
.
Background
The crystal topography is obtained with atomic force microscopy (AFM). From the AFM data, the probability of reflection (1) and the reflection directions (2) are computationally determined, for incidence angles ranging from 0° to 90°. Each LUT is computed for a given surface and reflector configuration. The reflection probability in the LUT combines two cases: directly reflected photons from the crystal surface and photons that are transmitted to the reflector surface and later re-enter the crystal. The key operations of the reflection process are the following: The angle between the incident photon (Old Momentum) and the surface normal are calculated. The probability of reflection is extracted from the first LUT. A Bernoulli test determines whether the photon is reflected or transmitted. In case of reflection two angles are drawn from the reflection direction LUT.
Parameterisation¶
In this section we describe how to use the parameterisation or “fast simulation” facilities of Geant4. Examples are provided in the examples/extended/parameterisations directory.
Generalities:¶
The Geant4 parameterisation facilities allow you to shortcut the detailed simulation in a given volume and for given particle types in order for you to provide your own implementation of the physics and of the detector response. This allows to make an alternative modelling of the physics processes, usually approximate and faster than the detailed simulation.
Overtaking the detailed Geant4 simulation (tracking) requires the user to specify 3 main components of the parameterisation:
Where the particles are parameterised;
Which particles are parameterised:
static conditions (particle type, PDG, charge, …)
dynamic conditions (energy, direction, …)
What happens instead of the detailed simulation:
where the particle is moved
what are the created secondaries
is the primary particle killed
what (and where) energy is deposited
Note: It is user responsibility to invoke Hit method of the sensitive detector
1. Where¶
Parameterisations are bound to a G4Region
object, which, in the
case of fast simulation can be called an envelope. A G4Region
has a
G4LogicalVolume
object (or a series of G4LogicalVolume
objects) as a root,
and the G4Region
is attached to this volume and all its ancestors:
G4Region(const G4String&) // constructor also registers to G4RegionStore
void G4Region::AddRootLogicalVolume (G4LogicalVolume*) // attach root volume to region
Envelopes often correspond to the outer volumes of (sub-)detectors: electromagnetic calorimeters, tracking chambers, etc.
With Geant4 it is also possible to define envelopes by overlaying a parallel (“ghost”) geometry as discussed in Parameterisation Using Ghost/Parallel Geometries.
2. Which particles¶
Parameterisation is usually specified only for certain particles. Those particles must
have attached G4FastSimulationManagerProcess
to their list of processes. For users
of modular physics lists (G4VModularPhysicsList
) — from which reference physics
lists (FTFP_BERT
, QGSP_BIC
,…) are derived — it is enough to use the helper class
G4FastSimulationPhysics
and activate the parameterisation for each particle type:
FTFP_BERT* physicsList = new FTFP_BERT; // G4VModularPhysicsList
auto fastSimulationPhysics = new G4FastSimulationPhysics(); // helper
fastSimulationPhysics->ActivateFastSimulation("e-"); // for envelope in mass geometry
fastSimulationPhysics->ActivateFastSimulation("pi+","pionGhostWorld"); // for envelopes in parallel geometry
physicsList->RegisterPhysics( fastSimulationPhysics ); // attach to the physics list
Geant4 will take into account the parameterisation process at tracking for any
step that starts in any root G4LogicalVolume
of the G4Region
that has been declared
as an envelope. It will proceed by first asking the available parameterisation models for
that particle (models with static conditions fulfilled by that particle). Existing
models will be checked (in the order they were added) if one of them (and only one)
wants to issue a trigger (meaning the dynamic conditions are met).
If yes, it will invoke its parameterisation. In this case, the tracking and physics
(detailed simulation) will not apply be applied to the particle in the step.
In case no trigger is issued, the simulation continues as usual, with other processes
being taken into account.
Parameterisations resembles a “user stepping action” but are more advanced because:
parameterisation code is considered only in the
G4Region
to which it is bound;parameterisation code is considered anywhere in the
G4Region
, that is, any volume in which the track is located (mother, daughter, sub-daughter, …);Geant4 will provide information to your parameterisation code about the current root volume of the
G4Region
in which the track is travelling.
3. What happens¶
Implementation of the parameterisation must be made deriving from the G4VFastSimulationModel
abstract class. Models are attached to the G4Region
and will be considered only if particle meets the selection criteria and is within the
geometrical hierarchy tree of the root logical volume of the G4Region
. The G4Region
is specified in the model constructor, together with the model name:
// constructor adds this model to G4FastSimulationManager of given envelope
G4VFastSimulationModel(const G4String&, G4Region*)
Selection criteria are to be defined in G4VFastSimulationModel
implementation:
// specify the static conditions (particle type, PDG, charge, ...)
virtual G4bool G4VFastSimulationModel::IsApplicable (const G4ParticleDefinition&) = 0
// specify the dynamic conditions (momentum, direction, position, distance to boundary, ...)
virtual G4bool G4VFastSimulationModel::ModelTrigger (const G4FastTrack&) = 0
As particle is not transported by Geant4, it is therefore responsibility of the model to describe what happens to the particle in place of the standard simulation: where the particle goes to (or to kill it), how its momentum is changed, what is the deposited energy (user needs to register all deposits in the sensitive detector), and what secondary particles are created. Those details should be implemented within:
// input information: G4FastTrack
// output information: G4FastStep
virtual G4bool G4VFastSimulationModel::DoIt(const G4FastTrack&, G4FastStep&) = 0
Geant4 contains an implementation of the Gflash parameterisation model (see more in
Gflash Parameterisation) and several examples in extended/parameterisations
directory.
Overview of Parameterisation Components¶
The Geant4 components which allow the implementation and control of parameterisations are:
G4VFastSimulationModel
This is the abstract class for the implementation of parameterisations. You must inherit from it to implement your concrete parameterisation model.
G4Region
As mentioned before, an envelope of parameterisation in Geant4 is a
G4Region
. The parameterisation is bound to theG4Region
by setting aG4FastSimulationManager
pointer to it.Fig. 19 shows how the
G4VFastSimulationModel
andG4FastSimulationManager
objects are bound to theG4Region
. Then for all root G4LogicalVolume’s held by the G4Region, the fast simulation code is active.G4FastSimulationManager
The
G4VFastSimulationModel
objects are attached to theG4Region
through aG4FastSimulationManager
. This object will manage the list of models and will message them at tracking time.G4FastSimulationManagerProcess
This is an implementation of
G4VProcess
. It invokes the parameterisation models if trigger conditions are met (particle is within an envelope, of certain type, etc.). It must be set in the process list of the particles you want to parameterise (e.g. using physics constructorG4FastSimulationPhysics
on top of any modular physics list - since 10.3 release). If added manually, one must remember that in presence of the parallel world, the ordering of processes matters.G4GlobalFastSimulationManager
This a singleton class which provides the management of the
G4FastSimulationManager
objects and some ghost facilities.
The G4VFastSimulationModel
Abstract Class¶
Constructors
The G4VFastSimulationModel
class has two constructors.
G4VFastSimulationModel(const G4String& aName)
:Here
aName
identifies the parameterisation model.G4VFastSimulationModel(const G4String& aName, G4Region*, G4bool IsUnique=false):
In addition to the model name, this constructor accepts a
G4Region
pointer. The neededG4FastSimulationManager
object is constructed if necessary, passing to it the G4Region pointer and the Boolean value. If it already exists, the model is simply added to this manager. Note that theG4VFastSimulationModel
object will not keep track of theG4Region
passed in the constructor. The Boolean argument is there for optimisation purposes: if you know that theG4Region
has a unique rootG4LogicalVolume
, uniquely placed, you can set the Boolean value totrue
.
Virtual methods
The G4VFastSimulationModel
has pure virtual methods which must
be overridden in your concrete class:
G4VFastSimulationModel(const G4String& aName):
Here aName identifies the parameterisation model.
G4bool ModelTrigger( const G4FastTrack&):
You must return
true
when the dynamic conditions to trigger your parameterisation are fulfilled.G4FastTrack
provides access to the currentG4Track
, gives simple access to the current rootG4LogicalVolume
related features (itsG4VSolid
, andG4AffineTransform
references between the global and the rootG4LogicalVolume
local coordinates systems) and simple access to the position and momentum expressed in the rootG4LogicalVolume
coordinate system. Using these quantities and theG4VSolid
methods, you can for example easily check how far you are from the rootG4LogicalVolume
boundary, or if the particle is entering or escaping the volume.G4bool IsApplicable(const G4ParticleDefinition&):
In your implementation, you must return
true
when your model is applicable to theG4ParticleDefinition
passed to this method. TheG4ParticleDefinition
provides all intrinsic particle information (mass, charge, spin, name …).If you want to implement a model which is valid only for certain particle types, it is recommended for efficiency that you use the static pointer of the corresponding particle classes.
As an example, in a model valid for gammas only, the
IsApplicable()
method should take the form:#include "G4Gamma.hh" G4bool MyGammaModel::IsApplicable(const G4ParticleDefinition& partDef) { return &partDef == G4Gamma::GammaDefinition(); }
void DoIt(const G4FastTrack&, G4FastStep&):
The details of your parameterisation will be implemented in this method. The
G4FastTrack
reference provides the input information, and the final state of the particles after parameterisation must be returned through theG4FastStep
reference (what is the changed energy and position, what are the secondaries, which particles are killed). Tracking for the final state secondary particles is requested after your parameterisation has been invoked.
The G4FastSimulationManager
Class¶
G4FastSimulationManager
functionalities regarding the use of ghost
volumes are explained in
Parameterisation Using Ghost/Parallel Geometries.
Constructor
G4FastSimulationManager(G4Region *anEnvelope, G4bool IsUnique=false):
This is the only constructor. You specify the
G4Region
by providing its pointer. TheG4FastSimulationManager
object will bind itself to thisG4Region
. If you know that thisG4Region
has a single rootG4LogicalVolume
, placed only once, you can set theIsUnique
boolean totrue
to allow some optimisation.Note that if you choose to use the
G4VFastSimulationModel(const G4String&, G4Region*, G4bool)
constructor for your model, theG4FastSimulationManager
will be constructed using the givenG4Region*
andG4bool
values of the model constructor.
Management of parameterisation models
The following two methods provide the usual management functions.
void AddFastSimulationModel( G4VFastSimulationModel*)
void RemoveFastSimulationModel( G4VFastSimulationModel*)
Messenger
To ease the communication with G4FastSimulationManager
a messenger class
was introduced. List of available commands:
/param/ // Fast Simulation print/control commands.
/param/showSetup // Show fast simulation setup (for each world: fast simulation manager
// process - which particles, region hierarchy - which models)
/param/listEnvelopes <ParticleName (default:all)> // List all the envelope names for a
// given particle (or for all particles if without parameters).
/param/listModels <EnvelopeName (default:all)> // List all the Model names for a given
// envelope (or for all envelopes if without parameters).
/param/listIsApplicable <ModelName (default:all)> // List all the Particle names a given
// model is applicable (or for all models if without parameters).
/param/ActivateModel <ModelName> // Activate a given Model.
/param/InActivateModel <ModelName> // InActivate a given Model.
The G4FastSimulationManagerProcess
Class¶
This G4VProcess
serves as an interface between the tracking and the
parameterisation. You usually don’t need to set it up directly, as you
can conveniently rely on the helper tool G4FastSimulationPhysics
(see Section 2. Which particles)
for this. At tracking (stepping) time, it collaborates with the
G4FastSimulationManager
of the current volume, if any, to allow the
models to trigger. If no manager exists or if no model issues a trigger,
the tracking goes on normally.
Parallel/Ghost Geometry
In order to register G4FastSimulationManagerProcess
operating on the regions from
the parallel geometry, it must be constructed passing the world name. Otherwise,
the mass geometry is used (and its navigator).
Moreover, G4FastSimulationManagerProcess
must be registered to G4ProcessManager
as continuous and discrete process, as it provides navigation in the ghost
world to limit the step on ghost boundaries. Hence it is important to maintain
the ordering of the along-step execution:
[n-3] ...
[n-2] Multiple Scattering
[n-1] G4FastSimulationManagerProcess
[ n ] G4Transportation
Since 10.3 release it is convenient to use G4FastSimulationPhysics
on top of
a modular physics list. To use fast simulation in the parallel geometry (register
G4FastSimulationManagerProcess
) one can do (from Par01 example):
// ----------------------------------------------
// -- PhysicsList and fast simulation activation:
// ----------------------------------------------
// -- Create a physics list (note : FTFP_BERT is a G4VModularPhysicsList
// -- which allows to use the subsequent G4FastSimulationPhysics tool to
// -- activate the fast simulation):
FTFP_BERT* physicsList = new FTFP_BERT;
// -- Create helper tool, used to activate the fast simulation:
G4FastSimulationPhysics* fastSimulationPhysics = new G4FastSimulationPhysics();
fastSimulationPhysics->BeVerbose();
// -- activation of fast simulation for particles having fast simulation models
// -- attached in the mass geometry:
fastSimulationPhysics->ActivateFastSimulation("e-");
fastSimulationPhysics->ActivateFastSimulation("e+");
fastSimulationPhysics->ActivateFastSimulation("gamma");
// -- activation of fast simulation for particles having fast simulation models
// -- attached in the parallel geometry:
fastSimulationPhysics->ActivateFastSimulation("pi+","pionGhostWorld");
fastSimulationPhysics->ActivateFastSimulation("pi-","pionGhostWorld");
// -- Attach the fast simulation physics constructor to the physics list:
physicsList->RegisterPhysics( fastSimulationPhysics );
// -- Finally passes the physics list to the run manager:
runManager->SetUserInitialization(physicsList);
If you wish to register manually G4FastSimulationManagerProcess
to all
the particles as a discrete and continuous process:
void MyPhysicsList::addParameterisation()
{
G4FastSimulationManagerProcess*
theFastSimulationManagerProcess = new G4FastSimulationManagerProcess();
theParticleIterator->reset();
while( (*theParticleIterator)() )
{
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
pmanager->AddProcess(theFastSimulationManagerProcess, -1, 0, 0);
}
}
The G4GlobalFastSimulationManager
Singleton Class¶
This class is a singleton which can be accessed as follows:
#include "G4GlobalFastSimulationManager.hh"
...
...
G4GlobalFastSimulationManager* globalFSM;
globalFSM = G4GlobalFastSimulationManager::getGlobalFastSimulationManager();
...
...
Presently, you will mainly need to use the G4GlobalFastSimulationManager
if you use ghost geometries.
Parameterisation Using Ghost/Parallel Geometries¶
In some cases, volumes of the tracking geometry do not allow envelopes to be defined. This may be the case with a geometry coming from a CAD system. Since such a geometry is flat, a parallel geometry must be used to define the envelopes.
Another interesting case involves defining an envelope which groups the electromagnetic and hadronic calorimeters of a detector into one volume. This may be useful when parameterising the interaction of charged pions. You will very likely not want electrons to see this envelope, which means that ghost geometries have to be organised by particle flavours.
Using ghost geometries implies some more overhead in the
parameterisation mechanism for the particles sensitive to ghosts, since
navigation is provided in the ghost geometry by the
G4FastSimulationManagerProcess
. Usually, however, only a few volumes
will be placed in this ghost world, so that the geometry computations
will remain rather cheap.
For details on how to use fast simulation with the parallel/world geometry please consult The G4FastSimulationManagerProcess Class.
Gflash Parameterisation¶
This section describes how to use the Gflash parameterisation. Gflash is a concrete implementation based on the equations and parameters of the original Gflash package from H1(hep-ex/0001020, Grindhammer & Peters, see physics manual) and uses the fast simulation facilities of Geant4 described above. Briefly, whenever a e-/e+ particle enters the calorimeter, it is parameterised if it has a minimum energy and the shower is expected to be contained in the calorimeter (or ‘parameterisation envelope’). If this is fulfilled the particle is killed, no secondaries are created, and instead the energy is deposited according to the Gflash equations. Examples, provided in examples/extended/parametrisation/gflash/, show how to interface Gflash to your application. The simulation time is measured, so the user can immediately see the speed increase resulting from the use of Gflash.
Using the Gflash Parameterisation¶
To use Gflash ‘out of the box’ the following steps are necessary:
The user must add the fast simulation process:
G4VModularPhysicsList* physicsList = new FTFP_BERT(); G4FastSimulationPhysics* fastSimulationPhysics = new G4FastSimulationPhysics(); fastSimulationPhysics->ActivateFastSimulation("e-"); physicsList->RegisterPhysics( fastSimulationPhysics );
The envelope in which the parameterisation should be performed must be created and specified (below:
G4Region* fRegion
) and theGFlashShowerModel
must be assigned to this region. Furthermore, the classesGFlashParticleBounds
(which provides thresholds for the parameterisation like minimal energy etc.),GflashHitMaker
(a helper class to generate hits in the sensitive detector) andGFlashHomoShowerParameterisation
(which does the computations) must be constructed and attached to theGFlashShowerModel
:G4Material* pbWO4 = nistManager->FindOrBuildMaterial("G4_PbWO4"); fFastShowerModel = new GFlashShowerModel("fFastShowerModel", fRegion); fParameterisation = new GFlashHomoShowerParameterisation(pbWO4); fFastShowerModel->SetParameterisation(*fParameterisation); fParticleBounds = new GFlashParticleBounds(); fFastShowerModel->SetParticleBounds(*fParticleBounds); fHitMaker = new GFlashHitMaker(); fFastShowerModel->SetHitMaker(*fHitMaker);
The user must set the material of the calorimeter, since the computation depends on the material.
It is mandatory to use
G4VGFlashSensitiveDetector
as (additional) base class for the sensitive detector:class ExGflash1SensitiveDetector: public G4VSensitiveDetector, public G4VGFlashSensitiveDetector
Here it is necessary to implement a separate interface, where the Gflash spots are processed:
virtual G4bool ProcessHits(G4GFlashSpot*aSpot,G4TouchableHistory*);
A separate interface is used, because the Gflash spots naturally contain less information than the full simulation.
Since the parameters in the Gflash package are taken from fits to full simulations with Geant3, with limited number of materials, and on specific range of electron energy, some re-tuning might be necessary for good agreement with Geant4 showers. Such parameterisation may be moreover made independent on material (currently atomic number Z) reducing the number of parameters. The tuning procedure is described in hep-ex/0001020 and there is on-going work to automate that procedure and make available in Geant4. For the moment, the only way to alter the parameters is to provide an implementation of GVFlashHomoShowerTuning class and pass it to the class that calculates the showers profiles, etc.:
GFlashHomoShowerParameterisation(G4Material * aMat, GVFlashHomoShowerTuning * aPar = 0);
For information, there is also a preliminary (still not validated) implementation of a parameterisation for sampling calorimeters. The user must specify the active and passive material, as well as the thickness of the active and passive layer. The sampling structure of the calorimeter is taken into account by using an “effective medium” to compute the shower shape.
All material properties needed are calculated automatically. If tuning is required, the user can pass his own parameter set in the class GFlashSamplingShowerTuning. Here the user can also set his calorimeter resolution.
All in all the constructor looks the following:
GFlashSamplingShowerParamterisation(G4Material * Mat1, G4Material * Mat2,G4double d1,G4double d2,
GVFlashSamplingShowerTuning * aPar = 0);
Transportation Process¶
A Transportation process is required for every particle which is tracked in a simulation. The Geant4 transportation processes are responsible for several key functions for tracking:
polling the Geometry Modeller via
G4Navigator
to obtain the distance to the next boundary for uncharged particles or charged particles in a volume / region without an electromagnetic field;handing off the tracking of charged particles in an EM field to
G4PropagatorInField
which finds either the endpoint of integration of the equations of motion of the particle or the state of the particle at the location in which it intersects with the next volume boundary;updating the time of flight of the particle, using the full step length (not the geometrical step length, which is reduced by multiple scattering.)
killing tracks which are found to loop inside a (magnetic) field, without making adequate progress after O(thousand) integration steps. A description of this is provided below.
Transportation comes in two flavours:
the ‘standard’
G4Transportation
process, used for most applications, andthe
G4CoupledTransportation
process, which is activated when multiple geometries are active.
Multiple geometries can be created in order to cope with different use cases:
when a mass overlay geometry is used to overlap a set of ‘top’ volume onto a complex existing geometry,
when the Geant4 scoring and/or biasing capabilities are activated.
The registration of the relevant Transportation process is handled by the
G4PhysicsListHelper
, which chooses the correct type depending on whether
any of the features which require parallel geometries have been used.
In brief there is one main difference between G4Transportation
and
G4CoupledTransportation
. The G4Transportation
process uses the
G4Navigator
of the Geant4 Geometry Modeller to obtain the distance to
the next boundary along a straight line (for a neutral particle, or a charged
particle in zero field). The G4CoupledTransportation
process uses the
G4PathFinder
class to obtain the shortest length to a boundary amongst the
geometries registered for the current particle - in effect multiplexing the
different geometries.
In addition the transportation processes estimates the time of flight for the current step. For a neutral particle or a charged particle inside a pure magnetic field, this is estimated from the initial and final velocity of the particle. This taking into account roughly the effect of energy loss from ionisation. Since the full path length is used (rather than the geometrical one) the path lengthening due to multiple scattering is also taken into account.
For a charged particle in an EM field with a non-zero electric component, or a gravity field, the time of flight is calculated taking into account the change in velocity.
For the propagation in an external field, electromagnetic or other, the
Transportation processes rely on the capabilities of G4PropagatorInField
and the integration methods detailed in the subsection
ElectroMagnetic Field
Note that the integration currently is done without taking into account either energy loss along the trajectory of motion or multiple scattering, which is applied independently at the endpoint (if it is not on a boundary.)
Further details about the caveats and control of transportation within a magnetic field are given in the Appendix: Transportation in Magnetic Field - Further Details.