5.2.  Physics Processes

Physics processes describe how particles interact with a material. Seven major categories of processes are provided by Geant4:

  1. electromagnetic ,

  2. hadronic ,

  3. decay ,

  4. photolepton-hadron ,

  5. optical ,

  6. parameterization and

  7. 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:

For each of the above DoIt methods G4VProcess provides a corresponding pure virtual GPIL method:

Other pure virtual methods in G4VProcess follow:

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 and PostStepDoIt methods.

example: transportation, ionisation(energy loss and delta ray)

G4VRestDiscreteProcess

Processes using both AtRestDoIt and PostStepDoIt methods.

example: positron annihilation, decay (both in flight and at rest)

G4VRestContinuousProcess

Processes using both AtRestDoIt and AlongStepDoIt 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,

which correspond to the three DoIt methods of G4VProcess. Each derived class should implement these methods.

5.2.1.  Electromagnetic Interactions

This section summarizes the electromagnetic physics processes which are installed in Geant4. For details on the implementation of these processes please refer to the Physics Reference Manual.

5.2.1.1.  "Standard" Electromagnetic Processes

The following is a summary of the standard electromagnetic processes available in Geant4.

  • Photon processes

    • Compton scattering (class name G4ComptonScattering)

    • Gamma conversion (also called pair production, class name G4GammaConversion)

    • Photo-electric effect (class name G4PhotoElectricEffect)

    • Muon pair production (class name G4GammaConversionToMuons)

  • Electron/positron processes

    • Ionisation and delta ray production (class name G4eIonisation)

    • Bremsstrahlung (class name G4eBremsstrahlung)

    • 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)

  • Hadron/ion processes

    • Ionisation (class name G4hIonisation)

    • Ionisation for ions (class name G4ionIonisation)

    • Ionisation for ions in low-density media (class name G4ionGasIonisation)

    • Ionisation for heavy exotic particles (class name G4hhIonisation)

    • Ionisation for classical magnetic monopole (class name G4mplIonisation)

  • Coulomb scattering processes

    • A general process in the sense that the same process/class is used to simulate the multiple scattering of the all charged particles (class name G4MultipleScattering)

    • Specialised process for more fast simulation the multiple scattering of muons and hadrons (class name G4hMultipleScattering)

    • Alternative process (beta-version) for the multiple scattering of muons (class name G4MuMultipleScattering)

    • 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:

    • Ionisation in thin absorbers (class name G4PAIModel)

An example of the registration of these processes in a physics list is given in Example 5.1, similar method is used in EM-builders of reference physics lists ($G4INSTALL/source/physics_lists/builders) and in EM examples ($G4INSTALL/examples/extended/electromagnetic).

Example 5.1.  Registration of standard electromagnetic processes

void PhysicsList::ConstructEM()

{

  theParticleIterator->reset();

  while( (*theParticleIterator)() ){

    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {

      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
      pmanager->AddDiscreteProcess(new G4ComptonScattering);
      pmanager->AddDiscreteProcess(new G4GammaConversion);

    } else if (particleName == "e-") {

      pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1);
      pmanager->AddProcess(new G4eIonisation,        -1, 2, 2);
      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3, 3);

    } else if (particleName == "e+") {

      pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1);
      pmanager->AddProcess(new G4eIonisation,        -1, 2, 2);
      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3, 3);
      pmanager->AddProcess(new G4eplusAnnihilation,   0,-1, 4);
      
    } else if( particleName == "mu+" || 
               particleName == "mu-"    ) {

      pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1);
      pmanager->AddProcess(new G4MuIonisation,       -1, 2, 2);
      pmanager->AddProcess(new G4MuBremsstrahlung,   -1, 3, 3);
               pmanager->AddProcess(new G4MuPairProduction,   -1, 4, 4);       

    } else if (particleName == "alpha" ||
               particleName == "He3" ||
               particleName == "GenericIon") {
      // ions with charge >= +2
      pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1);
      pmanager->AddProcess(new G4ionIonisation,      -1, 2, 2);
     
    } else if ((!particle->IsShortLived()) &&
               (particle->GetPDGCharge() != 0.0) && 
               (particle->GetParticleName() != "chargedgeantino")) {
      //all others charged particles except geantino and short-lived
      pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1);
      pmanager->AddProcess(new G4hIonisation,        -1, 2, 2);
            
    }
  }
}


Novice and extended electromagnetic examples illustrating the use of electromagnetic processes are available as part of the Geant4 release.

Options are available for steering the standard electromagnetic processes. These options may be invoked either by UI commands or by the interface class G4EmProcessOptions. This class has the following public methods:

  • SetLossFluctuations(G4bool)

  • SetSubCutoff(G4bool, const G4Region* r=0)

  • SetIntegral(G4bool)

  • SetMinSubRange(G4double)

  • SetMinEnergy(G4double)

  • SetMaxEnergy(G4double)

  • SetMaxEnergyForCSDARange(G4double)

  • SetMaxEnergyForMuons(G4double)

  • SetDEDXBinning(G4int)

  • SetDEDXBinningForCSDARange(G4int)

  • SetLambdaBinning(G4int)

  • SetStepFunction(G4double, G4double)

  • SetRandomStep(G4bool)

  • SetApplyCuts(G4bool)

  • SetBuildCSDARange(G4bool)

  • SetVerbose(G4int, const G4String name= "all")

  • SetLambdaFactor(G4double)

  • SetLinearLossLimit(G4double)

  • ActivateDeexcitation(G4bool val, const G4Region* r = 0)

  • SetMscStepLimitation(G4MscStepLimitType val)

  • SetMscLateralDisplacement(G4bool val)

  • SetSkin(G4double)

  • SetMscRangeFactor(G4double)

  • SetMscGeomFactor(G4double)

  • SetLPMFlag(G4bool)

  • SetBremsstrahlungTh(G4double)

The corresponding UI command can be accessed in the UI subdirectory "/process/eLoss". The following types of step limitation by multiple scattering are available:

  • fSimple - step limitation used in g4 7.1 version (used in QGSP_EMV Physics List)

  • fUseSafety - default

  • fUseDistanceToBoundary - advance method of step limitation used in EM examples, required parameter skin > 0, should be used for setup without 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,G4Region region=0)

  • GetRangeFromRestrictedDEDX(kinEnergy,particle,material,G4Region* region=0)

  • GetCSDARange(kinEnergy,particle,material,G4Region* region=0)

  • GetRange(kinEnergy,particle,material,G4Region* region=0)

  • GetKinEnergy(range,particle,material,G4Region* region=0)

  • GetCrosSectionPerVolume(kinEnergy,particle,material,G4Region* region=0)

  • GetMeanFreePath(kinEnergy,particle,material,G4Region* region=0)

  • PrintDEDXTable(particle)

  • PrintRangeTable(particle)

  • PrintInverseRangeTable(particle)

  • ComputeDEDX(kinEnergy,particle,process,material,cut=DBL_MAX)

  • ComputeElectronicDEDX(kinEnergy,particle,material,cut=DBL_MAX)

  • ComputeNuclearDEDX(kinEnergy,particle,material,cut=DBL_MAX)

  • ComputeTotalDEDX(kinEnergy,particle,material,cut=DBL_MAX)

  • ComputeCrosSectionPerVolume(kinEnergy,particle,process,material,cut=0)

  • ComputeCrosSectionPerAtom(kinEnergy,particle,process,Z,A,cut=0)

  • ComputeMeanFreePath(kinEnergy,particle,process,material,cut=0)

  • ComputeEnergyCutFromRangeCut(range,particle,material)

  • FindParticle(const G4String&)

  • FindMaterial(const G4String&)

  • FindRegion(const G4String&)

  • FindCouple(const G4Material*, const G4Region* region=0)

  • SetVerbose(G4int)

For these interfaces, particles, materials, or processes may be pointers or strings with names.

5.2.1.2.  Low Energy Electromagnetic Processes

The following is a summary of the Low Energy Electromagnetic processes available in Geant4. Further information is available in the homepage of the Geant4 Low Energy Electromagnetic Physics Working Group. The physics content of these processes is documented in Geant4 Physics Reference Manual and in other papers.

  • Photon processes

    • Compton scattering (class G4LowEnergyCompton)

    • Polarized Compton scattering (class G4LowEnergyPolarizedCompton)

    • Rayleigh scattering (class G4LowEnergyRayleigh)

    • Gamma conversion (also called pair production, class G4LowEnergyGammaConversion)

    • Photo-electric effect (classG4LowEnergyPhotoElectric)

  • Electron processes

    • Bremsstrahlung (class G4LowEnergyBremsstrahlung)

    • Ionisation and delta ray production (class G4LowEnergyIonisation)

  • Hadron and ion processes

    • Ionisation and delta ray production (class G4hLowEnergyIonisation)

An example of the registration of these processes in a physics list is given in Example 5.2.

Example 5.2.  Registration of electromagnetic low energy electron/photon processes.

void LowEnPhysicsList::ConstructEM()
{
  theParticleIterator->reset();

  while( (*theParticleIterator)() ){

    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if (particleName == "gamma") {

      theLEPhotoElectric   = new G4LowEnergyPhotoElectric();
      theLECompton         = new G4LowEnergyCompton();
      theLEGammaConversion = new G4LowEnergyGammaConversion();
      theLERayleigh        = new G4LowEnergyRayleigh();

      pmanager->AddDiscreteProcess(theLEPhotoElectric);
      pmanager->AddDiscreteProcess(theLECompton);
      pmanager->AddDiscreteProcess(theLERayleigh);
      pmanager->AddDiscreteProcess(theLEGammaConversion);

    }
    else if (particleName == "e-") {

      theLEIonisation = new G4LowEnergyIonisation();
      theLEBremsstrahlung = new G4LowEnergyBremsstrahlung();
      theeminusMultipleScattering = new G4MultipleScattering();

      pmanager->AddProcess(theeminusMultipleScattering,-1,1,1);
      pmanager->AddProcess(theLEIonisation,-1,2,2);
      pmanager->AddProcess(theLEBremsstrahlung,-1,-1,3);

    }
    else if (particleName == "e+") {

      theeplusMultipleScattering = new G4MultipleScattering();
      theeplusIonisation = new G4eIonisation();
      theeplusBremsstrahlung = new G4eBremsstrahlung();
      theeplusAnnihilation = new G4eplusAnnihilation();

      pmanager->AddProcess(theeplusMultipleScattering,-1,1,1);
      pmanager->AddProcess(theeplusIonisation,-1,2,2);
      pmanager->AddProcess(theeplusBremsstrahlung,-1,-1,3);
      pmanager->AddProcess(theeplusAnnihilation,0,-1,4);
    }
  }
}


Advanced examples illustrating the use of Low Energy Electromagnetic processes are available as part of the Geant4 release and are further documented here.

To run the Low Energy code for photon and electron electromagnetic processes, data files need to be copied by the user to his/her code repository. These files are distributed together with Geant4 release.

The user should set the environment variable G4LEDATA to the directory where he/she has copied the files.

Options are available for low energy electromagnetic processes for hadrons and ions in terms of public member functions of the G4hLowEnergyIonisation class:

  • SetHighEnergyForProtonParametrisation(G4double)

  • SetLowEnergyForProtonParametrisation(G4double)

  • SetHighEnergyForAntiProtonParametrisation(G4double)

  • SetLowEnergyForAntiProtonParametrisation(G4double)

  • SetElectronicStoppingPowerModel(const G4ParticleDefinition*,const G4String& )

  • SetNuclearStoppingPowerModel(const G4String&)

  • SetNuclearStoppingOn()

  • SetNuclearStoppingOff()

  • SetBarkasOn()

  • SetBarkasOff()

  • SetFluorescence(const G4bool)

  • ActivateAugerElectronProduction(G4bool)

  • SetCutForSecondaryPhotons(G4double)

  • SetCutForSecondaryElectrons(G4double)

The available models for ElectronicStoppingPower and NuclearStoppingPower are documented in the class diagrams.

Options are available for low energy electromagnetic processes for electrons in the G4LowEnergyIonisation class:

  • ActivateAuger(G4bool)

  • SetCutForLowEnSecPhotons(G4double)

  • SetCutForLowEnSecElectrons(G4double)

Options are available for low energy electromagnetic processes for electrons/positrons in the G4LowEnergyBremsstrahlung class, that allow the use of alternative bremsstrahlung angular generators:

  • SetAngularGenerator(G4VBremAngularDistribution* distribution);

  • SetAngularGenerator(const G4String& name);

Currently three angular generators are available: G4ModifiedTsai, 2BNGenerator and 2BSGenerator. G4ModifiedTsai is set by default, but it can be forced using the string "tsai". 2BNGenerator and 2BSGenerator can be set using the strings "2bs" and "2bn". Information regarding conditions of use, performance and energy limits of different models are available in the Physics Reference Manual and in the Geant4 Low Energy Electromagnetic Physics Working Group homepage.

Other options G4LowEnergyBremsstrahlung class are:

  • SetCutForLowEnSecPhotons(G4double)

Options can also be set in the G4LowEnergyPhotoElectric class, that allow the use of alternative photoelectron angular generators:

  • SetAngularGenerator(G4VPhotoElectricAngularDistribution* distribution);

  • SetAngularGenerator(const G4String& name);

Currently three angular generators are available: G4PhotoElectricAngularGeneratorSimple, G4PhotoElectricAngularGeneratorSauterGavrilla and G4PhotoElectricAngularGeneratorPolarized. G4PhotoElectricAngularGeneratorSimple is set by default, but it can be forced using the string "default". G4PhotoElectricAngularGeneratorSauterGavrilla and G4PhotoElectricAngularGeneratorPolarized can be set using the strings "standard" and "polarized". Information regarding conditions of use, performance and energy limits of different models are available in the Physics Reference Manual and in the Geant4 Low Energy Electromagnetic Physics Working Group homepage.

5.2.1.3.  Very Low energy Electromagnetic Processes (Geant4 DNA extension)

Geant4 low energy electromagnetic Physics processes have been extended down to energies of a few electronVolts suitable for the simulation of radiation effects in liquid water for applications at the cellular and sub-cellular level. These developments take place in the framework of the Geant4 DNA project [ http://www.ge.infn.it/geant4/dna ] and are fully described in the paper [ Chauvie2007 ].

Their implementation in Geant4 is based on the usage of innovative techniques first introduced in Monte Carlo simulation (policy-based class design), to ensure openness to future extension and evolution as well as flexibility of configuration in user applications. In this new design, a generic Geant4-DNA physics process is configured by template specialization in order to acquire physical properties (cross section, final state), using policy classes : a Cross Section policy class and a Final State policy class.

These processes apply to electrons, protons, hydrogen, alpha particles and their charge states.

Electron processes

  • Elastic scattering (two complementary models available depending on energy range)

    • Cross section policy class name, common to both models : G4CrossSectionElasticScreenedRutherford

    • Final state policy class names : G4FinalStateElasticScreenedRutherford or G4FinalStateElasticBrennerZaider

  • Excitation (one model)

    • Cross section policy class name : G4CrossSectionExcitationEmfietzoglou

    • Final state policy class name : G4FinalStateExcitationEmfietzoglou

  • Ionisation (one model)

    • Cross section policy class name : G4CrossSectionIonisationBorn

    • Final state policy class names : G4FinalStateIonisationBorn

Proton processes

  • Excitation (two complementary models available depending on energy range)

    • Cross section policy class name : G4CrossSectionExcitationMillerGreen

    • Final state policy class name : G4FinalStateExcitationMillerGreen

    • Cross section policy class name : G4CrossSectionExcitationBorn

    • Final state policy class name : G4FinalStateExcitationBorn

  • Ionisation (two complementary models available depending on energy range)

    • Cross section policy class name : G4CrossSectionIonisationRudd

    • Final state policy class name : G4FinalStateIonisationRudd

    • Cross section policy class name : G4CrossSectionIonisationBorn

    • Final state policy class name : G4FinalStateIonisationBorn

  • Charge decrease (one model)

    • Cross section policy class name : G4CrossSectionChargeDecrease

    • Final state policy class name : G4FinalStateChargeDecrease

Hydrogen processes

  • Ionisation (one model)

    • Cross section policy class name : G4CrossSectionIonisationRudd

    • Final state policy class name : G4FinalStateIonisationRudd

  • Charge increase (one model)

    • Cross section policy class name : G4CrossSectionChargeIncrease

    • Final state policy class name : G4FinalStateChargeIncrease

Helium (neutral) processes

  • Excitation (one model)

    • Cross section policy class name : G4CrossSectionExcitationMillerGreen

    • Final state policy class name : G4FinalStateExcitationMillerGreen

  • Ionisation (one model)

    • Cross section policy class name : G4CrossSectionIonisationRudd

    • Final state policy class name : G4FinalStateIonisationRudd

  • Charge increase (one model)

    • Cross section policy class name : G4CrossSectionChargeIncrease

    • Final state policy class name : G4FinalStateChargeIncrease

Helium+ (ionized once) processes

  • Excitation (one model)

    • Cross section policy class name : G4CrossSectionExcitationMillerGreen

    • Final state policy class name : G4FinalStateExcitationMillerGreen

  • Ionisation (one model)

    • Cross section policy class name : G4CrossSectionIonisationRudd

    • Final state policy class name : G4FinalStateIonisationRudd

  • Charge increase (one model)

    • Cross section policy class name : G4CrossSectionChargeIncrease

    • Final state policy class name : G4FinalStateChargeIncrease

  • Charge decrease (one model)

    • Cross section policy class name : G4CrossSectionChargeDecrease

    • Final state policy class name : G4FinalStateChargeDecrease

Helium++ (ionised twice) processes

  • Excitation (one model)

    • Cross section policy class name : G4CrossSectionExcitationMillerGreen

    • Final state policy class name : G4FinalStateExcitationMillerGreen

  • Ionisation (one model)

    • Cross section policy class name : G4CrossSectionIonisationRudd

    • Final state policy class name : G4FinalStateIonisationRudd

  • Charge decrease (one model)

    • Cross section policy class name : G4CrossSectionChargeDecrease

    • Final state policy class name : G4FinalStateChargeDecrease

An example of the registration of these processes in a physics list is given here below :

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

// Geant4 DNA header files

#include "G4DNAGenericIonsManager.hh"
#include "G4FinalStateProduct.hh"
#include "G4DNAProcess.hh"

#include "G4CrossSectionExcitationEmfietzoglou.hh"
#include "G4FinalStateExcitationEmfietzoglou.hh"

#include "G4CrossSectionElasticScreenedRutherford.hh"
#include "G4FinalStateElasticScreenedRutherford.hh"
#include "G4FinalStateElasticBrennerZaider.hh"

#include "G4CrossSectionExcitationBorn.hh"
#include "G4FinalStateExcitationBorn.hh"

#include "G4CrossSectionIonisationBorn.hh"
#include "G4FinalStateIonisationBorn.hh"

#include "G4CrossSectionIonisationRudd.hh"
#include "G4FinalStateIonisationRudd.hh"

#include "G4CrossSectionExcitationMillerGreen.hh"
#include "G4FinalStateExcitationMillerGreen.hh"

#include "G4CrossSectionChargeDecrease.hh"
#include "G4FinalStateChargeDecrease.hh"

#include "G4CrossSectionChargeIncrease.hh"
#include "G4FinalStateChargeIncrease.hh"

// Processes definition

typedef G4DNAProcess<G4CrossSectionElasticScreenedRutherford,G4FinalStateElasticScreenedRutherford> 
  ElasticScreenedRutherford;
typedef G4DNAProcess<G4CrossSectionElasticScreenedRutherford,G4FinalStateElasticBrennerZaider> 
  ElasticBrennerZaider;
typedef G4DNAProcess<G4CrossSectionExcitationEmfietzoglou,G4FinalStateExcitationEmfietzoglou> 
  ExcitationEmfietzoglou;
typedef G4DNAProcess<G4CrossSectionExcitationBorn,G4FinalStateExcitationBorn> 
  ExcitationBorn;
typedef G4DNAProcess<G4CrossSectionIonisationBorn,G4FinalStateIonisationBorn> 
  IonisationBorn;
typedef G4DNAProcess<G4CrossSectionIonisationRudd,G4FinalStateIonisationRudd> 
  IonisationRudd;
typedef G4DNAProcess<G4CrossSectionExcitationMillerGreen,G4FinalStateExcitationMillerGreen> 
  ExcitationMillerGreen;
typedef G4DNAProcess<G4CrossSectionChargeDecrease,G4FinalStateChargeDecrease> 
  ChargeDecrease;
typedef G4DNAProcess<G4CrossSectionChargeIncrease,G4FinalStateChargeIncrease> 
  ChargeIncrease;

// Processes registration

void MicrodosimetryPhysicsList::ConstructEM()
{
  theParticleIterator->reset();

  while( (*theParticleIterator)() ){

    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* processManager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();

    if (particleName == "e-") {
       processManager->AddDiscreteProcess(new ExcitationEmfietzoglou);
       processManager->AddDiscreteProcess(new ElasticScreenedRutherford);
       processManager->AddDiscreteProcess(new ElasticBrennerZaider);
       processManager->AddDiscreteProcess(new IonisationBorn);

    } else if ( particleName == "proton" ) {
       processManager->AddDiscreteProcess(new ExcitationMillerGreen);
       processManager->AddDiscreteProcess(new ExcitationBorn);
       processManager->AddDiscreteProcess(new IonisationRudd);
       processManager->AddDiscreteProcess(new IonisationBorn);
       processManager->AddDiscreteProcess(new ChargeDecrease);

    } else if ( particleName == "hydrogen" ) {
       processManager->AddDiscreteProcess(new IonisationRudd);
       processManager->AddDiscreteProcess(new ChargeIncrease);

    } else if ( particleName == "alpha" ) {
       processManager->AddDiscreteProcess(new ExcitationMillerGreen);
       processManager->AddDiscreteProcess(new IonisationRudd);
       processManager->AddDiscreteProcess(new ChargeDecrease);
    
    } else if ( particleName == "alpha+" ) {
       processManager->AddDiscreteProcess(new ExcitationMillerGreen);
       processManager->AddDiscreteProcess(new IonisationRudd);
       processManager->AddDiscreteProcess(new ChargeDecrease);
       processManager->AddDiscreteProcess(new ChargeIncrease);
    
    } else if ( particleName == "helium" ) {
       processManager->AddDiscreteProcess(new ExcitationMillerGreen);
       processManager->AddDiscreteProcess(new IonisationRudd);
       processManager->AddDiscreteProcess(new ChargeIncrease);
    }

  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Note that in the above example, "alpha" particles are helium atoms ionised twice and "helium" particles are neutral helium atoms. The definition of particles in the physics list may be for example implemented as follows :

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "G4DNAGenericIonsManager.hh"

void MicrodosimetryPhysicsList::ConstructBaryons()
{
  //  construct baryons ---

  // Geant4 DNA particles

  G4DNAGenericIonsManager * genericIonsManager;
  genericIonsManager=G4DNAGenericIonsManager::Instance();
  genericIonsManager->GetIon("alpha++");
  genericIonsManager->GetIon("alpha+");
  genericIonsManager->GetIon("helium");
  genericIonsManager->GetIon("hydrogen");

}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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.

5.2.2.  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.

5.2.2.1.  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, NeutronHPCrossSections, 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 whethee this is desirable or not for his particular problem.

5.2.2.2.  Hadrons at Rest

List of implemented "Hadron at Rest" processes

The following process classes have been implemented:

  • pi- absorption (class name G4PionMinusAbsorptionAtRest or G4PiMinusAbsorptionAtRest)

  • kaon- absorption (class name G4KaonMinusAbsorptionAtRest or G4KaonMinusAbsorption)

  • neutron capture (class name G4NeutronCaptureAtRest)

  • anti-proton annihilation (class name G4AntiProtonAnnihilationAtRest)

  • anti-neutron annihilation (class name G4AntiNeutronAnnihilationAtRest)

  • mu- capture (class name G4MuonMinusCaptureAtRest)

  • alternative CHIPS model for any negativly charged particle (class name G4QCaptureAtRest)

Obviously the last process does not, strictly speaking, deal with a ``hadron at rest''. It does, nonetheless, share common features with the others in the above list because of the implementation model chosen. The differences between the alternative implementation for kaon and pion absorption concern the fast part of the emitted particle spectrum. G4PiMinusAbsorptionAtRest, and G4KaonMinusAbsorptionAtRest focus especially on a good description of this part of the spectrum.

Implementation Interface to Geant4

All of these classes are derived from the abstract class G4VRestProcess. In addition to the constructor and destructor methods, the following public methods of the abstract class have been implemented for each of the above six processes:

  • AtRestGetPhysicalInteractionLength( const G4Track&, G4ForceCondition* )

    This method returns the time taken before the interaction actually occurs. In all processes listed above, except for muon capture, a value of zero is returned. For the muon capture process the muon capture lifetime is returned.

  • AtRestDoIt( const G4Track&, const G4Step&)

    This method generates the secondary particles produced by the process.

  • IsApplicable( const G4ParticleDefinition& )

    This method returns the result of a check to see if the process is possible for a given particle.

Example of how to use a hadron at rest process

Including a ``hadron at rest'' process for a particle, a pi- for example, into the Geant4 system is straightforward and can be done in the following way:

  • create a process:

           theProcess = new G4PionMinusAbsorptionAtRest();
        

  • register the process with the particle's process manager:

           theParticleDef = G4PionMinus::PionMinus();
           G4ProcessManager* pman = theParticleDef->GetProcessManager();
           pman->AddRestProcess( theProcess );
        

5.2.2.3.  Hadrons in Flight

What processes do you need?

For hadrons in motion, there are four physics process classes. Table 5.1 shows each process and the particles for which it is relevant.

G4HadronElasticProcess pi+, pi-, K+, K0S, K0L, K-, p, p-bar, n, n-bar, lambda, lambda-bar, Sigma+, Sigma-, Sigma+-bar, Sigma--bar, Xi0, Xi-, Xi0-bar, Xi--bar
G4HadronInelasticProcess pi+, pi-, K+, K0S, K0L, K-, p, p-bar, n, n-bar, lambda, lambda-bar, Sigma+, Sigma-, Sigma+-bar, Sigma--bar, Xi0, Xi-, Xi0-bar, Xi--bar
G4HadronFissionProcess all
G4CaptureProcess n, n-bar

Table 5.1.  Hadronic processes and relevant particles.


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:

   G4LEProtonInelastic *theProtonIE = new G4LEProtonInelastic();

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.

A sample case for the proton's inelastic interaction model class is shown in Example 5.3, where G4LEProtonInelastic.hh is the name of the include file:

Example 5.3.  An example of a proton inelastic interaction model class.

 ----------------------------- include file ------------------------------------------

#include "G4InelasticInteraction.hh"
 class G4LEProtonInelastic : public G4InelasticInteraction
 {
 public:
    G4LEProtonInelastic() : G4InelasticInteraction()
    {
      SetMinEnergy( 0.0 );
      SetMaxEnergy( 25.*GeV );
    }
    ~G4LEProtonInelastic() { }
    G4ParticleChange *ApplyYourself( const G4Track &aTrack,
                                     G4Nucleus &targetNucleus );
 private:
    void CascadeAndCalculateMomenta( required arguments );
 };

 ----------------------------- source file ------------------------------------------

 #include "G4LEProtonInelastic.hh"
 G4ParticleChange *
  G4LEProton Inelastic::ApplyYourself( const G4Track &aTrack,
                                       G4Nucleus &targetNucleus )
  {
    theParticleChange.Initialize( aTrack );
    const G4DynamicParticle *incidentParticle = aTrack.GetDynamicParticle();
    // create the target particle
    G4DynamicParticle *targetParticle = targetNucleus.ReturnTargetParticle();
    CascadeAndCalculateMomenta( required arguments )
    { ... }
    return &theParticleChange;
  }


The CascadeAndCalculateMomenta function is the bulk of the model and is to be provided by the model's creator. It should determine what secondary particles are produced in the interaction, calculate the momenta for all the particles, and put this information into the ParticleChange object which is returned.

The G4LEProtonInelastic class derives 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.

Three types of hadronic shower models have been implemented: parametrisation driven models, data driven models, and theory driven models.

  • Parametrisation driven models are used for all processes pertaining to particles coming to rest, and interacting with the nucleus. For particles in flight, two sets of models exist for inelastic scattering; low energy, and high energy models. Both sets are based originally on the GHEISHA package of Geant3.21, and the original approaches to primary interaction, nuclear excitation, intra-nuclear cascade and evaporation is kept. The models are located in the sub-directories hadronics/models/low_energy and hadronics/models/high_energy. The low energy models are targeted towards energies below 20 GeV; the high energy models cover the energy range from 20 GeV to O(TeV). Fission, capture and coherent elastic scattering are also modeled through parametrised 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.

5.2.3.  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.

5.2.3.1.  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 in G4DecayTable, and

  • using thePreAssignedDecayProducts of G4DynamicParticle

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 Section 2.5.3) is shown in Example 5.4.

Example 5.4.  Registration of the decay process to particles in the ConstructPhysics method of PhysicsList.

#include "G4Decay.hh"
void ExN02PhysicsList::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);
    }
  }
}


5.2.3.2.  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 Section 5.3.2.

5.2.3.3.  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.

5.2.4.  Photolepton-hadron Processes

To be delivered.

5.2.5.  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. Not so, if the user wishes to start optical photons as primary particles. In this case, the user must set the linear polarization using particle gun methods, the General Particle Source, or his/her 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 and Rayleigh scattering. Processes which produce optical photons include the Cerenkov effect, transition radiation 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 constants or they may be expressed as a function of the photon's wavelength. 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 composed of instantiations of the class G4MPVEntry. The G4MPVEntry is a pair of numbers, which in the case of an optical property, are the photon momentum and corresponding property value. The G4MaterialPropertyVector is implemented as a G4std::vector, with the sorting operation defined as MPVEntry1 < MPVEntry2 == photon_momentum1 < photon_momentum2. This results in all G4MaterialPropertyVectors being sorted in ascending order of photon momenta. 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 Example 5.5.

Example 5.5.  Optical properties added to a G4MaterialPropertiesTable and linked to a G4Material

const G4int NUMENTRIES = 32;

G4double ppckov[NUMENTRIES] = {2.034*eV, ......, 4.136*eV};
G4double rindex[NUMENTRIES] = {1.3435, ......, 1.3608};
G4double absorption[NUMENTRIES] = {344.8*cm, ......, 1450.0*cm];

G4MaterialPropertiesTable *MPT = new G4MaterialPropertiesTable();

MPT -> AddConstProperty("SCINTILLATIONYIELD",100./MeV);

MPT -> AddProperty("RINDEX",ppckov,rindex,NUMENTRIES};
MPT -> AddProperty("ABSLENGTH",ppckov,absorption,NUMENTRIES};

scintillator -> SetMaterialPropertiesTable(MPT);


5.2.5.1.  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 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 frequently very large number of secondaries produced in a single step (about 300/cm in water), compelled the idea in GEANT3.21 of suspending the primary particle until all its progeny have been tracked. Despite the fact that GEANT4 employs dynamic memory allocation and thus does not suffer from the limitations of GEANT3.21 with its fixed large initial ZEBRA store, GEANT4 nevertheless provides for an analogous functionality with the public method SetTrackSecondariesFirst. An example of the registration of the Cerenkov process is given in Example 5.6.

Example 5.6.  Registration of the Cerenkov process in PhysicsList.

#include "G4Cerenkov.hh"

void ExptPhysicsList::ConstructOp(){

  G4Cerenkov*   theCerenkovProcess = new G4Cerenkov("Cerenkov");

  G4int MaxNumPhotons = 300;

  theCerenkovProcess->SetTrackSecondariesFirst(true);
  theCerenkovProcess->SetMaxNumPhotonsPerStep(MaxNumPhotons);

  theParticleIterator->reset();
  while( (*theParticleIterator)() ){
    G4ParticleDefinition* particle = theParticleIterator->value();
    G4ProcessManager* pmanager = particle->GetProcessManager();
    G4String particleName = particle->GetParticleName();
    if (theCerenkovProcess->IsApplicable(*particle)) {
      pmanager->AddContinuousProcess(theCerenkovProcess);
    }
  }
}


5.2.5.2.  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. 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 Example 5.7

Example 5.7.  Specification of scintillation properties in DetectorConstruction.

  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 Example 5.8. 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.

Example 5.8.  Implementation of the scintillation process in PhysicsList.

  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π with a random linear polarization and at times characteristic for the scintillation component.

5.2.5.3.  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 momentum. WLSCOMPONENT is the relative emission spectrum of the material as a function of the photon's momentum, and WLSTIMECONSTANT accounts for any time delay which may occur between absorption and re-emission of the photon. An example is shown in Example 5.9.

Example 5.9.  Specification of WLS properties in DetectorConstruction.

 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*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.

5.2.5.4.  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 absorpted 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, σ/ω, is proportional to cos2(θ), where θ 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 attenutation length specified by the user, the program automatically calls the RayleighAttenuationLengthGenerator which calculates it for 10 degrees Celsius liquid water.

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 beteen two dielectic 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.

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 probabilites 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 [A. Levin and C. Moisan, A More Physical Approach to Model the Surface Treatment of Scintillation Counters and its Implementation into DETECT, TRIUMF Preprint TRI-PP-96-64, Oct. 1996] of the DETECT program [G.F. Knoll, T.F. Knoll and T.M. Henderson, Light Collection Scintillation Detector Composites for Neutron Detection, IEEE Trans. Nucl. Sci., 35 (1988) 872.]. 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 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.

Example 5.10.  Dielectric-dielectric surface properties defined via the G4OpticalSurface.

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. [GEANT Detector Description and Simulation Tool, Application Software Group, Computing and Networks Division, CERN, PHYS260-6 tp 260-7.].

Example 5.11.  Dielectric metal surface properties defined via the G4OpticalSurface.

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, but does not otherwise specify the model reflection probability constants, the default becomes Lambertian reflection.

5.2.6.  Parameterization

In this section we describe how to use the parameterization or "fast simulation" facilities of GEANT4. Examples are provided in the examples/novice/N05 directory.

5.2.6.1.  Generalities:

The Geant4 parameterization facilities allow you to shortcut the detailed tracking 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.

Parameterisations are bound to a G4Region object, which, in the case of fast simulation is also called an envelope. Prior to release 8.0, parameterisations were bound to a G4LogicalVolume, the root of a volume hierarchy. These root volumes are now attributes of the G4Region. Envelopes often correspond to the volumes of sub-detectors: electromagnetic calorimeters, tracking chambers, etc. With GEANT4 it is also possible to define envelopes by overlaying a parallel or "ghost" geometry as discussed in Section 5.2.6.7.

In GEANT4, parameterisations have three main features. You must specify:

  • the particle types for which your parameterisation is valid;

  • the dynamics conditions for which your parameterisation is valid and must be triggered;

  • the parameterisation itself: where the primary will be killed or moved, whether or not to create it or create secondaries, etc., and where the detector response will be computed.

GEANT4 will message your parameterisation code for each step starting in any root G4LogicalVolume (including daughters. sub-daughters, etc. of this volume) of the G4Region. It will proceed by first asking the available parameterisations for the current particle type if one of them (and only one) wants to issue a trigger. If so it will invoke its parameterisation. In this case, the tracking will not apply physics to the particle in the step. Instead, the UserSteppingAction will be invoked.

Parameterisations look like a "user stepping action" but are more advanced because:

  • parameterisation code is messaged only in the G4Region to which it is bound;

  • parameterisation code is messaged anywhere in the G4Region, that is, any volume in which the track is located;

  • GEANT4 will provide information to your parameterisation code about the current root volume of the G4Region in which the track is travelling.

5.2.6.2.  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.

G4FastSimulationManager

The G4VFastSimulationModel objects are attached to the G4Region through a G4FastSimulationManager. This object will manage the list of models and will message them at tracking time.

G4Region/Envelope

As mentioned before, an envelope in GEANT4 is a G4Region. The parameterisation is bound to the G4Region by setting a G4FastSimulationManager pointer to it.

The figure below shows how the G4VFastSimulationModel and G4FastSimulationManager objects are bound to the G4Region. Then for all root G4LogicalVolume's held by the G4Region, the fast simulation code is active.

G4FastSimulationManagerProcess

This is a G4VProcess. It provides the interface between the tracking and the parameterisation. It must be set in the process list of the particles you want to parameterise.

G4GlobalFastSimulationManager

This a singleton class which provides the management of the G4FastSimulationManager objects and some ghost facilities.

5.2.6.3.  The G4VFastSimulationModel Abstract Class

Constructors:

The G4VFastSimulationModel class has two constructors. The second one allows you to get started quickly:

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 needed G4FastSimulationManager 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 the G4VFastSimulationModel object will not keep track of the G4Region passed in the constructor. The boolean argument is there for optimization purposes: if you know that the G4Region has a unique root G4LogicalVolume, uniquely placed, you can set the boolean value to "true".

Virtual methods:

The G4VFastSimulationModel has three pure virtual methods which must be overriden 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 current G4Track, gives simple access to the current root G4LogicalVolume related features (its G4VSolid, and G4AffineTransform references between the global and the root G4LogicalVolume local coordinates systems) and simple access to the position and momentum expressed in the root G4LogicalVolume coordinate system. Using these quantities and the G4VSolid methods, you can for example easily check how far you are from the root G4LogicalVolume boundary.

G4bool IsApplicable( const G4ParticleDefinition&):

In your implementation, you must return "true" when your model is applicable to the G4ParticleDefinition passed to this method. The G4ParticleDefinition 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();
      }
      

G4bool ModelTrigger( const G4FastTrack&):

You must return "true" when the dynamic conditions to trigger your parameterisation are fulfilled. The G4FastTrack provides access to the current G4Track, gives simple access to envelope related features (G4LogicalVolume, G4VSolid, and G4AffineTransform references between the global and the envelope local coordinates systems) and simple access to the position and momentum expressed in the envelope coordinate system. Using these quantities and the G4VSolid methods, you can for example easily check how far you are from the envelope boundary.

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 the G4FastStep reference. Tracking for the final state particles is requested after your parameterisation has been invoked.

5.2.6.4.  The G4FastSimulationManager Class:

G4FastSimulationManager functionnalities regarding the use of ghost volumes are explained in Section 5.2.6.7.

Constructor:

G4FastSimulationManager( G4Region *anEnvelope, G4bool IsUnique=false):

This is the only constructor. You specify the G4Region by providing its pointer. The G4FastSimulationManager object will bind itself to this G4Region. If you know that this G4Region has a single root G4LogicalVolume, placed only once, you can set the IsUnique boolean to "true" to allow some optimization.

Note that if you choose to use the G4VFastSimulationModel(const G4String&, G4Region*, G4bool) constructor for your model, the G4FastSimulationManager will be constructed using the given G4Region* and G4bool values of the model constructor.

G4VFastSimulationModel object management:

The following two methods provide the usual management functions.

  • void AddFastSimulationModel( G4VFastSimulationModel*)

  • RemoveFastSimulationModel( G4VFastSimulationModel*)

Interface with the G4FastSimulationManagerProcess:

This is described in the User's Guide for Toolkit Developers ( section 3.9.6 )

5.2.6.5.  The G4FastSimulationManagerProcess Class

This G4VProcess serves as an interface between the tracking and the parameterisation. At tracking 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.

In the present implementation, you must set this process in the G4ProcessManager of the particles you parameterise to enable your parameterisation.

The processes ordering is:

    [n-3] ...
    [n-2] Multiple Scattering
    [n-1] G4FastSimulationManagerProcess
    [ n ] G4Transportation

This ordering is important if you use ghost geometries, since the G4FastSimulationManagerProcess will provide navigation in the ghost world to limit the step on ghost boundaries.

The G4FastSimulationManager must be added to the process list of a particle as a continuous and discrete process if you use ghost geometries for this particle. You can add it as a discrete process if you don't use ghosts.

The following code registers the G4FastSimulationManagerProcess with 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);
    }
}

5.2.6.6.  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 GlobalFastSimulationManager if you use ghost geometries.

5.2.6.7.  Parameterisation Using Ghost 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 parameterizing the interaction of charged pions. You will very likely not want electrons to see this envelope, which means that ghost geometries have to be organized 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.

In the existing implementation (temporary implementation with G4Region but before parallel geometry implementation), you may only consider ghost G4Regions with just one root G4LogicalVolume. The G4GlobalFastSimulationManager provides the construction of the ghost geometry by making first an empty "clone" of the world for tracking provided by the construct() method of your G4VUserDetectorConstruction concrete class. You provide the placement of the G4Region root G4LogicalVolume relative to the ghost world coordinates in the G4FastSimulationManager objects. A ghost G4Region is recognized by the fact that its associated G4FastSimulationManager retains a non-empty list of placements.

The G4GlobalFastSimulationManager will then use both those placements and the IsApplicable() methods of the models attached to the G4FastSimulationManager objects to build the flavour-dependant ghost geometries.

Then at the beginning of the tracking of a particle, the appropriate ghost world, if any, will be selected.

The steps required to build one ghost G4Region are:

  1. built the ghost G4Region : myGhostRegion;

  2. build the root G4LogicalVolume: myGhostLogical, set it to myGhostRegion;

  3. build a G4FastSimulationManager object, myGhostFSManager, giving myGhostRegion as argument of the constructor;

  4. give to the G4FastSimulationManager the placement of the myGhostLogical, by invoking for the G4FastSimulationManager method:

         AddGhostPlacement(G4RotationMatrix*, const G4ThreeVector&);
        

    or:

         AddGhostPlacement(G4Transform3D*);
        

    where the rotation matrix and translation vector of the 3-D transformation describe the placement relative to the ghost world coordinates.

  5. build your G4VFastSimulationModel objects and add them to the myGhostFSManager. The IsApplicable() methods of your models will be used by the G4GlobalFastSimulationManager to build the ghost geometries corresponding to a given particle type.

  6. Invoke the G4GlobalFastSimulationManager method:

         G4GlobalFastSimulationManager::getGlobalFastSimulationManager()->
    
              CloseFastSimulation();
        

This last call will cause the G4GlobalFastSimulationManager to build the flavour-dependent ghost geometries. This call must be done before the RunManager closes the geometry. (It is foreseen that the run manager in the future will invoke the CloseFastSimulation() to synchronize properly with the closing of the geometry).

Visualization facilities are provided for ghosts geometries. After the CloseFastSimulation() invocation, it is possible to ask for the drawing of ghosts in an interactive session. The basic commands are:

  • /vis/draw/Ghosts particle_name

    which makes the drawing of the ghost geometry associated with the particle specified by name in the command line.

  • /vis/draw/Ghosts

    which draws all the ghost geometries.

5.2.6.8.  Gflash Parameterization

This section describes how to use the Gflash library. Gflash is a concrete parameterization which is 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 parameterized if it has a minimum energy and the shower is expected to be contained in the calorimeter (or " parameterization envelope"). If this is fulfilled the particle is killed, as well as all secondaries, and the energy is deposited according to the Gflash equations. An example, provided in examples/extended/parametrisation/gflash/, shows 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.

5.2.6.9.  Using the Gflash Parameterisation

To use Gflash "out of the box" the following steps are necessary:

  • The user must add the fast simulation process to his process manager:

      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 envelope in which the parameterization should be performed must be specified (below: G4Region m_calo_region) and the GFlashShowerModel must be assigned to this region. Furthermore, the classes GFlashParticleBounds (which provides thresholds for the parameterization like minimal energy etc.), GflashHitMaker(a helper class to generate hits in the sensitive detector) and GFlashHomoShowerParamterisation (which does the computations) must be constructed (by the user at the moment) and assigned to the GFlashShowerModel. Please note that at the moment only homogeneous calorimeters are supported.

      m_theFastShowerModel = new GFlashShowerModel("fastShowerModel",m_calo_region);
      m_theParametrisation = new GFlashHomoShowerParamterisation(matManager->getMaterial(mat));
      m_theParticleBounds  = new GFlashParticleBounds();
      m_theHMaker          = new GFlashHitMaker();
      m_theFastShowerModel->SetParametrisation(*m_theParametrisation);
      m_theFastShowerModel->SetParticleBounds(*m_theParticleBounds) ;
      m_theFastShowerModel->SetHitMaker(*m_theHMaker);
        

    The user must also 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 ExGflashSensitiveDetector: public G4VSensitiveDetector ,public G4VGFlashSensitiveDetector 
        

    Here it is necessary to implement a separate interface, where the GFlash spots are processed.

      (ProcessHits(G4GFlashSpot*aSpot ,G4TouchableHistory* ROhist))
        

    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, some retuning might be necessary for good agreement with Geant4 showers. For experiment-specific geometries some retuning might be necessary anyway. The tuning is quite complicated since there are many parameters (some correlated) and cannot be described here (see again hep-ex/0001020). For brave users the Gflash framework already forsees the possibility of passing a class with the (users) parameters,GVFlashHomoShowerTuning, to the GFlashHomoShowerParamterisation constructor. The default parameters are the original Gflash parameters:

GFlashHomoShowerParameterisation(G4Material * aMat, GVFlashHomoShowerTuning * aPar = 0);

Now there is also a preliminary implemenation of a parameterization 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);

An implementation of some tools that should help the user to tune the parameterization is forseen.

5.2.7.  Transportation Process

To be delivered by J. Apostolakis ().