Physics Reference Manual

Version: geant4 10.1 (5 December 2014)

Contents

I Introduction 1
1 Introduction 2
1.1 Scope of This Manual . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Definition of Terms . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 3

2 Monte Carlo Methods 4

2.1 Status of this document . . . . . . . . . . . . . . . . . . . . . 5

3 Particle Transport 6
3.1 Transportation . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.1.1 Status of This Document . . . . . . . . . . . . . . . . . 7
3.2 True Step Length . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2.1 The Interaction Length or Mean Free Path . . . . . . . 8
3.2.2 Determination of the Interaction Point . . . . . . . . . 9
3.2.3 Step Limitations . . . . . . . . . . . . . . . . . . . . . 9
3.2.4 Updating the Particle Time . . . . . . . . . . . . . . . 10
3.2.5 Status of This Document . . . . . . . . . . . . . . . . . 10

II Particle Decay 11
4 Decay 12
4.1 Mean Free Path for Decay in Flight . . . . . . . . . . . . . . . 12
4.2 Branching Ratios and Decay Channels . . . . . . . . . . . . . 12
4.2.1 G4PhaseSpaceDecayChannel . . . . . . . . . . . . . . . 13
4.2.2 G4DalitzDecayChannel . . . . . . . . . . . . . . . . . . 13
4.2.3 Muon Decay . . . . . . . . . . . . . . . . . . . . . . . . 14
4.2.4 Leptonic Tau Decay . . . . . . . . . . . . . . . . . . . 15
4.2.5 Kaon Decay . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 16

-10
III Electromagnetic Interactions 17
5 Gamma Incident 18
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.1.1 General Interfaces . . . . . . . . . . . . . . . . . . . . . 19
5.1.2 Status of This Document . . . . . . . . . . . . . . . . . 19
5.2 Photoelectric Effect . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.1 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.2 Final State . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.3 Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2.4 Status of this document . . . . . . . . . . . . . . . . . 23
5.3 Compton scattering . . . . . . . . . . . . . . . . . . . . . . . . 24
5.3.1 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 24
5.3.2 Sampling the Final State . . . . . . . . . . . . . . . . . 25
5.3.3 Atomic shell effects . . . . . . . . . . . . . . . . . . . . 26
5.3.4 Status of This Document . . . . . . . . . . . . . . . . . 27
5.4 Gamma Conversion into an Electron - Positron Pair . . . . . . 28
5.4.1 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 28
5.4.2 Final State . . . . . . . . . . . . . . . . . . . . . . . . 32
5.4.3 Ultra-Relativistic Model . . . . . . . . . . . . . . . . . 33
5.4.4 Status of This Document . . . . . . . . . . . . . . . . . 33
5.5 Gamma Conversion into a Muon - Anti-mu Pair . . . . . . . . 35
5.5.1 Cross Section and Energy Sharing . . . . . . . . . . . . 35
5.5.2 Parameterization of the Total Cross Section . . . . . . 38
5.5.3 Multi-differential Cross Section and Angular Variables 40
5.5.4 Procedure for the Generation of + Pairs . . . . . . 42
5.5.5 Status of this document . . . . . . . . . . . . . . . . . 49

6 Elastic scattering 50
6.1 Multiple Scattering . . . . . . . . . . . . . . . . . . . . . . . . 51
6.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 51
6.1.2 Definition of Terms . . . . . . . . . . . . . . . . . . . . 52
6.1.3 Path Length Correction . . . . . . . . . . . . . . . . . 54
6.1.4 Angular Distribution . . . . . . . . . . . . . . . . . . . 56
6.1.5 Determination of the Model Parameters . . . . . . . . 56
6.1.6 Step Limitation Algorithm . . . . . . . . . . . . . . . . 58
6.1.7 Boundary Crossing Algorithm . . . . . . . . . . . . . . 60
6.1.8 Implementation Details . . . . . . . . . . . . . . . . . . 61
6.1.9 Status of this document . . . . . . . . . . . . . . . . . 63
6.2 Discrete Processes for Charged Particles . . . . . . . . . . . . 65
6.2.1 Status of This Document . . . . . . . . . . . . . . . . . 66
 6.3 Single Scattering . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.3.1 Coulomb Scattering . . . . . . . . . . . . . . . . . . . . 67
6.3.2 Implementation Details . . . . . . . . . . . . . . . . . . 68
6.3.3 Status of This Document . . . . . . . . . . . . . . . . . 69
6.4 Ion Scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.4.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.4.2 Implementation Details . . . . . . . . . . . . . . . . . . 74
6.4.3 Status of this document . . . . . . . . . . . . . . . . . 74
6.5 Single Scattering, Screened Coulomb Potential and NIEL . . . 76
6.5.1 NucleusNucleus Interactions . . . . . . . . . . . . . . 76
6.5.2 Nuclear Stopping Power . . . . . . . . . . . . . . . . . 78
6.5.3 Non-Ionizing Energy Loss due to Coulomb Scattering . 81
6.5.4 G4IonCoulombScatteringModel . . . . . . . . . . . . . 82
6.5.5 The Method . . . . . . . . . . . . . . . . . . . . . . . . 82
6.5.6 Implementation Details . . . . . . . . . . . . . . . . . . 83
6.5.7 Status of This Document . . . . . . . . . . . . . . . . . 83
6.6 Electron Screened Single Scattering and NIEL . . . . . . . . . 85
6.6.1 Scattering Cross Section of Electrons on Nuclei . . . . 85
6.6.2 Nuclear Stopping Power of Electrons . . . . . . . . . . 94
6.6.3 Non-Ionizing Energy-Loss of Electrons . . . . . . . . . 95
6.7 G4eSingleScatteringModel . . . . . . . . . . . . . . . . . . . . 96
6.7.1 The method . . . . . . . . . . . . . . . . . . . . . . . . 97
6.7.2 Implementation Details . . . . . . . . . . . . . . . . . . 99
6.8 Status of this Document . . . . . . . . . . . . . . . . . . . . . 99

7 Energy loss of Charged Particles 101

7.1 Mean Energy Loss . . . . . . . . . . . . . . . . . . . . . . . . 102
7.1.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.1.2 General Interfaces . . . . . . . . . . . . . . . . . . . . . 103
7.1.3 Step-size Limit . . . . . . . . . . . . . . . . . . . . . . 103
7.1.4 Run Time Energy Loss Computation . . . . . . . . . . 105
7.1.5 Energy Loss by Heavy Charged Particles . . . . . . . . 106
7.1.6 Status of This Document . . . . . . . . . . . . . . . . . 107
7.2 Energy Loss Fluctuations . . . . . . . . . . . . . . . . . . . . . 109
7.2.1 Fluctuations in Thick Absorbers . . . . . . . . . . . . . 109
7.2.2 Fluctuations in Thin Absorbers . . . . . . . . . . . . . 110
7.2.3 Width Correction Algorithm . . . . . . . . . . . . . . . 112
7.2.4 Sampling of Energy Loss . . . . . . . . . . . . . . . . . 112
7.2.5 Status of This Document . . . . . . . . . . . . . . . . . 113
7.3 Correcting the Cross Section for Energy Variation . . . . . . 114
7.3.1 Status of This Document . . . . . . . . . . . . . . . . . 115
 7.4 Conversion from Cut in Range to Energy Threshold . . . . . . 116
7.4.1 Status of This Document . . . . . . . . . . . . . . . . . 118
7.5 Photoabsorption ionization model . . . . . . . . . . . . . . . . 119
7.5.1 Cross Section for Ionizing Collisions . . . . . . . . . . . 119
7.5.2 Energy Loss Simulation . . . . . . . . . . . . . . . . . 121
7.5.3 Photoabsorption Cross Section at Low Energies . . . . 122
7.5.4 Status of this document . . . . . . . . . . . . . . . . . 123

8 Electron and Positron Incident 124

8.1 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.1.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.1.2 Continuous Energy Loss . . . . . . . . . . . . . . . . . 125
8.1.3 Total Cross Section per Atom and Mean Free Path . . 127
8.1.4 Simulation of Delta-ray Production . . . . . . . . . . . 128
8.1.5 Status of this document . . . . . . . . . . . . . . . . . 129
8.2 Bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . . . . . 130
8.2.1 Seltzer-Berger bremsstrahlung model . . . . . . . . . . 130
8.2.2 Bremsstrahlung of high-energy electrons . . . . . . . . 133
8.2.3 Status of this document . . . . . . . . . . . . . . . . . 136
8.3 Positron - Electron Annihilation . . . . . . . . . . . . . . . . . 138
8.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 138
8.3.2 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 138
8.3.3 Sampling the final state . . . . . . . . . . . . . . . . . 138
8.3.4 Sampling the Gamma Energy . . . . . . . . . . . . . . 139
8.3.5 Status of This Document . . . . . . . . . . . . . . . . . 140
8.4 Positron Annihilation into + Pair in Media . . . . . . . . . 141
8.4.1 Total Cross Section . . . . . . . . . . . . . . . . . . . . 141
8.4.2 Sampling of Energies and Angles . . . . . . . . . . . . 141
8.4.3 Status of this document . . . . . . . . . . . . . . . . . 144
8.5 Positron Annihilation into Hadrons . . . . . . . . . . . . . . . 146
8.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 146
8.5.2 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 146
8.5.3 Sampling the final state . . . . . . . . . . . . . . . . . 146
8.5.4 Status of this document . . . . . . . . . . . . . . . . . 147

9 Low Energy Livermore 148

9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
9.1.1 Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
9.1.2 Data Sources . . . . . . . . . . . . . . . . . . . . . . . 149
9.1.3 Distribution of the Data Sets . . . . . . . . . . . . . . 150
9.1.4 Calculation of Total Cross Sections . . . . . . . . . . . 151
 9.1.5 Status of the document . . . . . . . . . . . . . . . . . . 151
9.2 Compton Scattering . . . . . . . . . . . . . . . . . . . . . . . 153
9.2.1 Total Cross Section . . . . . . . . . . . . . . . . . . . . 153
9.2.2 Sampling of the Final State . . . . . . . . . . . . . . . 153
9.2.3 Status of the document . . . . . . . . . . . . . . . . . . 154
9.3 Compton Scattering by Linearly Polarized Gamma Rays . . . 155
9.3.1 The Cross Section . . . . . . . . . . . . . . . . . . . . . 155
9.3.2 Angular Distribution . . . . . . . . . . . . . . . . . . . 155
9.3.3 Polarization Vector . . . . . . . . . . . . . . . . . . . . 155
9.3.4 Unpolarized Photons . . . . . . . . . . . . . . . . . . . 156
9.3.5 Status of this document . . . . . . . . . . . . . . . . . 156
9.4 Rayleigh Scattering . . . . . . . . . . . . . . . . . . . . . . . . 157
9.4.1 Total Cross Section . . . . . . . . . . . . . . . . . . . . 157
9.4.2 Sampling of the Final State . . . . . . . . . . . . . . . 157
9.4.3 Status of this document . . . . . . . . . . . . . . . . . 157
9.5 Gamma Conversion . . . . . . . . . . . . . . . . . . . . . . . . 159
9.5.1 Total cross-section . . . . . . . . . . . . . . . . . . . . 159
9.5.2 Sampling of the final state . . . . . . . . . . . . . . . . 159
9.5.3 Status of the document . . . . . . . . . . . . . . . . . . 160
9.6 Triple Gamma Conversion . . . . . . . . . . . . . . . . . . . . 161
9.6.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 161
9.6.2 Azimuthal Distribution for Electron Recoil . . . . . . . 161
9.6.3 Monte Carlo Simulation of the Asymptotic Expression 161
9.6.4 Algorithm for Non Polarized Radiation . . . . . . . . . 162
9.6.5 Algorithm for Polarized Radiation . . . . . . . . . . . . 164
9.6.6 Sampling of Energy . . . . . . . . . . . . . . . . . . . . 166
9.6.7 Status of This Document . . . . . . . . . . . . . . . . . 167
9.7 Photoelectric effect . . . . . . . . . . . . . . . . . . . . . . . . 168
9.7.1 Cross sections . . . . . . . . . . . . . . . . . . . . . . . 168
9.7.2 Sampling of the final state . . . . . . . . . . . . . . . . 168
9.7.3 Angular distribution of the emmited photoelectron . . 168
9.7.4 Status of the document . . . . . . . . . . . . . . . . . . 170
9.8 Electron ionisation . . . . . . . . . . . . . . . . . . . . . . . . 171
9.8.1 Status of the document . . . . . . . . . . . . . . . . . . 172
9.9 Bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . . . . . 173
9.9.1 Bremsstrahlung angular distributions . . . . . . . . . . 174
9.9.2 Status of the document . . . . . . . . . . . . . . . . . . 178
10 Low Energy Penelope 180
10.1 Penelope physics . . . . . . . . . . . . . . . . . . . . . . . . . 181
10.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 181
10.1.2 Compton scattering . . . . . . . . . . . . . . . . . . . . 181
10.1.3 Rayleigh scattering . . . . . . . . . . . . . . . . . . . . 183
10.1.4 Gamma conversion . . . . . . . . . . . . . . . . . . . . 184
10.1.5 Photoelectric effect . . . . . . . . . . . . . . . . . . . . 186
10.1.6 Bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . 187
10.1.7 Ionisation . . . . . . . . . . . . . . . . . . . . . . . . . 189
10.1.8 Positron Annihilation . . . . . . . . . . . . . . . . . . . 195
10.1.9 Status of the document . . . . . . . . . . . . . . . . . . 196

11 Monash University low energy photon processes 198

11.1 Monash Low Energy Model . . . . . . . . . . . . . . . . . . . 199
11.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 199
11.1.2 Physics and Simulation . . . . . . . . . . . . . . . . . . 199
11.1.3 Status of the document . . . . . . . . . . . . . . . . . . 201

12 Charged Hadron Incident 202

12.1 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
12.1.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 203
12.1.2 Continuous Energy Loss . . . . . . . . . . . . . . . . . 203
12.1.3 Nuclear Stopping . . . . . . . . . . . . . . . . . . . . . 208
12.1.4 Total Cross Section per Atom . . . . . . . . . . . . . . 208
12.1.5 Simulating Delta-ray Production . . . . . . . . . . . . 209
12.1.6 Ion Effective Charge . . . . . . . . . . . . . . . . . . . 210
12.1.7 Status of this document . . . . . . . . . . . . . . . . . 211
12.2 Low energy extentions . . . . . . . . . . . . . . . . . . . . . . 213
12.2.1 Energy losses of slow negative particles . . . . . . . . . 213
12.2.2 Energy losses of hadrons in compounds . . . . . . . . . 213
12.2.3 Fluctuations of energy losses of hadrons . . . . . . . . 214
12.2.4 ICRU 73-based energy loss model . . . . . . . . . . . . 216
12.2.5 Status of this document . . . . . . . . . . . . . . . . . 217

13 Muon Incident 219

13.1 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
13.1.1 Status of this document . . . . . . . . . . . . . . . . . 221
13.2 Bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . . . . . 222
13.2.1 Differential Cross Section . . . . . . . . . . . . . . . . . 222
13.2.2 Continuous Energy Loss . . . . . . . . . . . . . . . . . 223
13.2.3 Total Cross Section . . . . . . . . . . . . . . . . . . . . 223
 13.2.4 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . 224
13.2.5 Status of this document . . . . . . . . . . . . . . . . . 225
13.3 Positron - Electron Pair Production by Muons . . . . . . . . . 227
13.3.1 Differential Cross Section . . . . . . . . . . . . . . . . . 227
13.3.2 Total Cross Section and Restricted Energy Loss . . . . 230
13.3.3 Sampling of Positron - Electron Pair Production . . . . 231
13.3.4 Status of this document . . . . . . . . . . . . . . . . . 232
13.4 Muon Photonuclear Interaction . . . . . . . . . . . . . . . . . 234
13.4.1 Differential Cross Section . . . . . . . . . . . . . . . . . 234
13.4.2 Sampling . . . . . . . . . . . . . . . . . . . . . . . . . 235
13.4.3 Status of this document . . . . . . . . . . . . . . . . . 237

14 Atomic Relaxation 239

14.1 Atomic relaxation . . . . . . . . . . . . . . . . . . . . . . . . . 240
14.1.1 Fluorescence . . . . . . . . . . . . . . . . . . . . . . . . 240
14.1.2 Auger process . . . . . . . . . . . . . . . . . . . . . . . 241
14.1.3 PIXE . . . . . . . . . . . . . . . . . . . . . . . . . . . 241
14.1.4 Status of the document . . . . . . . . . . . . . . . . . . 242

15 Geant4-DNA 243
15.1 Geant4-DNA processes and models . . . . . . . . . . . . . . . 244
15.1.1 Status of the document . . . . . . . . . . . . . . . . . . 244

16 Microelectronics 245
16.1 The MicroElec extension for microelectronics applications . . . 246
16.1.1 Status of the document . . . . . . . . . . . . . . . . . . 247

17 Polarized Electron/Positron/Gamma Incident 248

17.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249
17.1.1 Stokes vector . . . . . . . . . . . . . . . . . . . . . . . 249
17.1.2 Transfer matrix . . . . . . . . . . . . . . . . . . . . . . 251
17.1.3 Coordinate transformations . . . . . . . . . . . . . . . 252
17.1.4 Polarized beam and material . . . . . . . . . . . . . . . 253
17.1.5 Status of this document . . . . . . . . . . . . . . . . . 255
17.2 Ionization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256
17.2.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 256
17.2.2 Total cross section and mean free path . . . . . . . . . 256
17.2.3 Sampling the final state . . . . . . . . . . . . . . . . . 258
17.2.4 Status of this document . . . . . . . . . . . . . . . . . 261
17.3 Positron - Electron Annihilation . . . . . . . . . . . . . . . . . 263
17.3.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 263
 17.3.2 Total cross section and mean free path . . . . . . . . . 263
17.3.3 Sampling the final state . . . . . . . . . . . . . . . . . 265
17.3.4 Annihilation at Rest . . . . . . . . . . . . . . . . . . . 267
17.3.5 Status of this document . . . . . . . . . . . . . . . . . 268
17.4 Polarized Compton scattering . . . . . . . . . . . . . . . . . . 269
17.4.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 269
17.4.2 Total cross section and mean free path . . . . . . . . . 269
17.4.3 Sampling the final state . . . . . . . . . . . . . . . . . 270
17.4.4 Status of this document . . . . . . . . . . . . . . . . . 273
17.5 Polarized Bremsstrahlung for electron and positron . . . . . . 274
17.5.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 274
17.5.2 Polarization in gamma conversion and bremsstrahlung 274
17.5.3 Polarization transfer to the photon . . . . . . . . . . . 275
17.5.4 Polarization transfer to the lepton . . . . . . . . . . . . 276
17.5.5 Status of this document . . . . . . . . . . . . . . . . . 278
17.6 Polarized Gamma conversion into an electronpositron pair . . 280
17.6.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 280
17.6.2 Polarization transfer . . . . . . . . . . . . . . . . . . . 280
17.6.3 Status of this document . . . . . . . . . . . . . . . . . 281
17.7 Polarized Photoelectric Effect . . . . . . . . . . . . . . . . . . 282
17.7.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 282
17.7.2 Polarization transfer . . . . . . . . . . . . . . . . . . . 282
17.7.3 Status of this document . . . . . . . . . . . . . . . . . 284

18 X-Ray Production 285

18.1 Transition radiation . . . . . . . . . . . . . . . . . . . . . . . . 286
18.1.1 Relationship of Transition Rad to Cherenkov Rad . . . 286
18.1.2 Calculating the X-ray Transition Radiation Yield . . . 287
18.1.3 Simulating X-ray Transition Radiation Production . . . 289
18.1.4 Status of this document . . . . . . . . . . . . . . . . . 292
18.2 Scintillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293
18.2.1 Status of this document . . . . . . . . . . . . . . . . . 293

18.3 Cerenkov Effect . . . . . . . . . . . . . . . . . . . . . . . . . . 294
18.3.1 Status of this document . . . . . . . . . . . . . . . . . 295
18.4 Synchrotron Radiation . . . . . . . . . . . . . . . . . . . . . . 296
18.4.1 Photon spectrum . . . . . . . . . . . . . . . . . . . . . 296
18.4.2 Validity . . . . . . . . . . . . . . . . . . . . . . . . . . 297
18.4.3 Direct inversion/generation of photon energy spectrum 298
18.4.4 Properties of the Power Spectra . . . . . . . . . . . . . 301
18.4.5 Status of This Document . . . . . . . . . . . . . . . . . 302
19 Optical Photons 304
19.1 Interactions of optical photons . . . . . . . . . . . . . . . . . . 305
19.1.1 Physics processes for optical photons . . . . . . . . . . 305
19.1.2 Photon polarization . . . . . . . . . . . . . . . . . . . . 306
19.1.3 Tracking of the photons . . . . . . . . . . . . . . . . . 307
19.1.4 Mie Scattering in Henyey-Greensterin Approximation . 310

20 Phonon-Lattice Interactions 313

20.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314
20.2 Phonon Propagation . . . . . . . . . . . . . . . . . . . . . . . 314
20.3 Lattice Parameters . . . . . . . . . . . . . . . . . . . . . . . . 315
20.4 Scattering and Mode Mixing . . . . . . . . . . . . . . . . . . . 315
20.5 Anharmonic Downconversion . . . . . . . . . . . . . . . . . . . 316
20.6 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316

21 Precision multi-scale modeling 318

21.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 319
21.2 Impact ionisation by hadrons and PIXE . . . . . . . . . . . . 319
21.3 Status of the document . . . . . . . . . . . . . . . . . . . . . . 327

22 Shower Parameterizations 328

22.1 Gflash Shower Parameterizations . . . . . . . . . . . . . . . . 329
22.1.1 Parameterization Ansatz . . . . . . . . . . . . . . . . . 329
22.1.2 Longitudinal Shower Profiles . . . . . . . . . . . . . . 329
22.1.3 Radial Shower Profiles . . . . . . . . . . . . . . . . . . 330
22.1.4 Gflash Performance . . . . . . . . . . . . . . . . . . . . 331
22.1.5 Status of this document . . . . . . . . . . . . . . . . . 332

IV Hadronic Interactions 333

23 Total Reaction Cross Section in Nucleus-nucleus Reactions 334
23.1 Sihver Formula . . . . . . . . . . . . . . . . . . . . . . . . . . 334
23.2 Kox and Shen Formulae . . . . . . . . . . . . . . . . . . . . . 335
23.3 Tripathi formula . . . . . . . . . . . . . . . . . . . . . . . . . 337
23.4 Representative Cross Sections . . . . . . . . . . . . . . . . . . 339
23.5 Tripathi Formula for light Systems . . . . . . . . . . . . . . 339
23.6 Status of this document . . . . . . . . . . . . . . . . . . . . . 340

24 Coherent elastic scattering 344

24.1 Nucleon-Nucleon elastic Scattering . . . . . . . . . . . . . . . 344
25 Hadron-nucleus Elastic Scattering at Medium/High Energy345
25.1 Method of Calculation . . . . . . . . . . . . . . . . . . . . . . 345
25.2 Status of this document . . . . . . . . . . . . . . . . . . . . . 348

26 Interactions of Stopping Particles 362

26.1 Complementary parameterised and theoretical treatment . . . 362
26.1.1 Pion absorption at rest . . . . . . . . . . . . . . . . . . 363

27 Parametrization Driven Models 365

27.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365
27.2 Low Energy Model . . . . . . . . . . . . . . . . . . . . . . . . 366
27.3 High Energy Model . . . . . . . . . . . . . . . . . . . . . . . . 367
27.3.1 Initial Interaction . . . . . . . . . . . . . . . . . . . . . 367
27.3.2 Intra-nuclear Cascade . . . . . . . . . . . . . . . . . . . 367
27.3.3 High Energy Cascading . . . . . . . . . . . . . . . . . . 368
27.3.4 High Energy Cluster Production . . . . . . . . . . . . . 373
27.3.5 Medium Energy Cascading . . . . . . . . . . . . . . . . 375
27.3.6 Medium Energy Cluster Production . . . . . . . . . . . 375
27.3.7 Elastic and Quasi-elastic Scattering . . . . . . . . . . . 375
27.4 Status of this document . . . . . . . . . . . . . . . . . . . . . 376

28 Parton string model. 377

28.1 Reaction initial state simulation. . . . . . . . . . . . . . . . . . 377
28.1.1 Allowed projectiles and bombarding energy range . . . 377
28.1.2 MC initialization procedure for nucleus. . . . . . . . . 377
28.1.3 Random choice of the impact parameter. . . . . . . . . 379
28.2 Sample of collision participants in nuclear collisions. . . . . . . 379
28.2.1 MC procedure to define collision participants. . . . . . 379
28.2.2 Separation of hadron diffraction excitation. . . . . . . . 380
28.3 Longitudinal string excitation . . . . . . . . . . . . . . . . . . 381
28.3.1 Hadronnucleon inelastic collision . . . . . . . . . . . . 381
28.3.2 The diffractive string excitation . . . . . . . . . . . . . 381
28.3.3 The string excitation by parton exchange . . . . . . . . 381
28.3.4 Transverse momentum sampling . . . . . . . . . . . . . 382
28.3.5 Sampling x-plus and x-minus . . . . . . . . . . . . . . 382
28.3.6 The diffractive string excitation . . . . . . . . . . . . . 382
28.3.7 The string excitation by parton rearrangement . . . . . 383
28.4 Longitudinal string decay. . . . . . . . . . . . . . . . . . . . . 384
28.4.1 Hadron production by string fragmentation. . . . . . . 384
28.4.2 The hadron formation time and coordinate. . . . . . . 385
28.5 Status of this document . . . . . . . . . . . . . . . . . . . . . 385
29 Fritiof (FTF) Model 387
29.1 Main assumptions of the FTF model . . . . . . . . . . . . . . 388
29.2 General properties of hadron-nucleon interactions . . . . . . . 391
29.2.1 p-interactions . . . . . . . . . . . . . . . . . . . . . . 391
29.2.2 + p-interactions . . . . . . . . . . . . . . . . . . . . . . 393
29.2.3 pp-interactions . . . . . . . . . . . . . . . . . . . . . . 394
29.2.4 K + p- and K p-interactions . . . . . . . . . . . . . . . 395
29.2.5 Protonanti-proton interactions . . . . . . . . . . . . . 397
29.3 Hadron-nucleon process cross section . . . . . . . . . . . . . . 399
29.3.1 Total, elastic and inelastic hadron-nucleon cross sections399
29.3.2 Cross sections of quark exchange processes . . . . . . . 401
29.3.3 Cross sections of anti-proton processes . . . . . . . . . 401
29.3.4 Cross sections of diffractive and non-diffractive processes402
29.4 Simulation of hadron-nucleon interactions . . . . . . . . . . . 405
29.4.1 Simulation of meson-nucleon and nucleon-nucleon interactions405
29.4.2 Simulation of anti-baryon-nucleon interactions . . . . . 408
29.5 Flowchart of the FTF model . . . . . . . . . . . . . . . . . . . 409
29.6 Simulation of nuclear interactions . . . . . . . . . . . . . . . . 411
29.6.1 Sampling of intra-nuclear collisions . . . . . . . . . . . 411
29.6.2 Reggeon cascading . . . . . . . . . . . . . . . . . . . . 417
29.6.3 Fermi motion of nuclear nucleons . . . . . . . . . . . 423
29.6.4 Excitation energy of nuclear residuals . . . . . . . . . . 427
29.7 Status of this document . . . . . . . . . . . . . . . . . . . . . 427

30 Chiral Invariant Phase Space Decay 431

30.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 431
30.2 Fundamental Concepts . . . . . . . . . . . . . . . . . . . . . . 434
30.3 Code Development . . . . . . . . . . . . . . . . . . . . . . . . 435
30.4 Nucleon-Antinucleon Annihilation at Rest . . . . . . . . . . . 436
30.4.1 Meson Production . . . . . . . . . . . . . . . . . . . . 437
30.4.2 Baryon Production . . . . . . . . . . . . . . . . . . . . 441
30.5 Nuclear Pion Capture Below Delta(3,3) . . . . . . . . . . . . . 447
30.6 Modeling of real and virtual photon interactions . . . . . . . . 464
30.7 Chiral invariant phase-space decay . . . . . . . . . . . . . . . 473
30.8 Neutrino-nuclear interactions . . . . . . . . . . . . . . . . . . . 473
30.9 Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 478
30.10Status of this document . . . . . . . . . . . . . . . . . . . . . 479
31 Bertini Intranuclear Cascade Model in Geant4 484
31.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484
31.2 The Geant4 Cascade Model . . . . . . . . . . . . . . . . . . . 485
31.2.1 Model Limits . . . . . . . . . . . . . . . . . . . . . . . 485
31.2.2 Intranuclear Cascade Model . . . . . . . . . . . . . . . 485
31.2.3 Nuclear Model . . . . . . . . . . . . . . . . . . . . . . 486
31.2.4 Pre-equilibrium Model . . . . . . . . . . . . . . . . . . 488
31.2.5 Break-up models . . . . . . . . . . . . . . . . . . . . . 488
31.2.6 Evaporation Model . . . . . . . . . . . . . . . . . . . . 489
31.3 Interfacing Bertini implementation . . . . . . . . . . . . . . . 489
31.4 Status of this document . . . . . . . . . . . . . . . . . . . . . 490

32 The Geant4 Binary Cascade 493

32.1 Modeling overview . . . . . . . . . . . . . . . . . . . . . . . . 493
32.1.1 The transport algorithm . . . . . . . . . . . . . . . . . 493
32.1.2 The description of the target nucleus and fermi motion 494
32.1.3 Optical and phenomenological potentials . . . . . . . . 495
32.1.4 Pauli blocking simulation . . . . . . . . . . . . . . . . . 496
32.1.5 The scattering term . . . . . . . . . . . . . . . . . . . 496
32.1.6 Total inclusive cross-sections . . . . . . . . . . . . . . 497
32.1.7 Channel cross-sections . . . . . . . . . . . . . . . . . . 497
32.1.8 Mass dependent resonance width and partial width . . 498
32.1.9 Resonance production cross-section in the t-channel . . 498
32.1.10 Nucleon Nucleon elastic collisions . . . . . . . . . . . . 499
32.1.11 Generation of transverse momentum . . . . . . . . . . 499
32.1.12 Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . 500
32.1.13 The escaping particle and coherent effects . . . . . . . 500
32.1.14 Light ion reactions . . . . . . . . . . . . . . . . . . . . 501
32.1.15 Transition to pre-compound modeling . . . . . . . . . . 501
32.1.16 Calculation of excitation energies and residuals . . . . 502
32.2 Comparison with experiments . . . . . . . . . . . . . . . . . . 502

33 Quantum Molecular Dynamics for Heavy Ions 511

33.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . 512
33.2 Ion-ion Implementation . . . . . . . . . . . . . . . . . . . . . . 514
33.3 Cross Sections . . . . . . . . . . . . . . . . . . . . . . . . . . . 515
33.4 Status of this document . . . . . . . . . . . . . . . . . . . . . 515

34 Abrasion-ablation Model 517

34.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 517
34.2 Initial nuclear dynamics and impact parameter . . . . . . . . . 518
 34.3 Abrasion process . . . . . . . . . . . . . . . . . . . . . . . . . 519
34.4 Abraded nucleon spectrum . . . . . . . . . . . . . . . . . . . . 521
34.5 De-excitation of nuclear pre-fragments by standard G4 . . . . 522
34.6 De-excitation of nuclear pre-fragments by nuclear ablation . . 523
34.7 Definition of the functions P and F used in the abrasion model 524
34.8 Status of this document . . . . . . . . . . . . . . . . . . . . . 525

35 Electromagnetic Dissociation Model 528

35.1 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 528
35.2 Status of this document . . . . . . . . . . . . . . . . . . . . . 531

36 Precompound model. 532

36.1 Reaction initial state. . . . . . . . . . . . . . . . . . . . . . . . 532
36.2 Simulation of pre-compound reaction . . . . . . . . . . . . . . 532
36.2.1 Statistical equilibrium condition . . . . . . . . . . . . . 533
36.2.2 Level density of excited (n-exciton) states . . . . . . . 533
36.2.3 Transition probabilities . . . . . . . . . . . . . . . . . . 533
36.2.4 Emission probabilities for nucleons . . . . . . . . . . . 535
36.2.5 Emission probabilities for complex fragments . . . . . . 535
36.2.6 The total probability . . . . . . . . . . . . . . . . . . . 536
36.2.7 Calculation of kinetic energies for emitted particle . . . 536
36.2.8 Parameters of residual nucleus. . . . . . . . . . . . . . 536
36.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 536

37 Evaporation Model 538

37.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 538
37.2 Evaporation model . . . . . . . . . . . . . . . . . . . . . . . . 538
37.2.1 Cross sections for inverse reactions . . . . . . . . . . . 539
37.2.2 Coulomb barriers . . . . . . . . . . . . . . . . . . . . . 539
37.2.3 Level densities . . . . . . . . . . . . . . . . . . . . . . . 540
37.2.4 Maximum energy available for evaporation . . . . . . . 540
37.2.5 Total decay width . . . . . . . . . . . . . . . . . . . . . 541
37.3 GEM model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 541
37.4 Nuclear fission . . . . . . . . . . . . . . . . . . . . . . . . . . . 543
37.4.1 The fission total probability . . . . . . . . . . . . . . . 543
37.4.2 The fission barrier . . . . . . . . . . . . . . . . . . . . 543
37.5 Photon evaporation . . . . . . . . . . . . . . . . . . . . . . . . 544
37.5.1 Computation of probability . . . . . . . . . . . . . . . 544
37.5.2 Discrete photon evaporation . . . . . . . . . . . . . . . 544
37.5.3 Internal conversion electron emission . . . . . . . . . . 545
37.6 Sampling procedure . . . . . . . . . . . . . . . . . . . . . . . . 546
38 Fission model. 549
38.1 Reaction initial state. . . . . . . . . . . . . . . . . . . . . . . . 549
38.2 Fission process simulation. . . . . . . . . . . . . . . . . . . . . 549
38.2.1 Atomic number distribution of fission products. . . . . 549
38.2.2 Charge distribution of fission products. . . . . . . . . . 551
38.2.3 Kinetic energy distribution of fission products. . . . . . 551
38.2.4 Calculation of the excitation energy of fission products. 552
38.2.5 Excited fragment momenta. . . . . . . . . . . . . . . . 552

39 Fermi break-up model. 554

39.1 Fermi break-up simulation for light nuclei . . . . . . . . . . . . 554
39.1.1 Allowed channels . . . . . . . . . . . . . . . . . . . . . 554
39.1.2 Break-up probability . . . . . . . . . . . . . . . . . . . 555
39.1.3 Fragment characteristics . . . . . . . . . . . . . . . . . 556
39.1.4 Sampling procedure . . . . . . . . . . . . . . . . . . . . 556

40 Multifragmentation model. 558

40.1 Multifragmentation process simulation. . . . . . . . . . . . . . 558
40.1.1 Multifragmentation probability. . . . . . . . . . . . . . 558
40.1.2 Direct simulation of low multiplicity disintegration . . 560
40.1.3 Fragment multiplicity distribution. . . . . . . . . . . . 561
40.1.4 Atomic number distribution of fragments. . . . . . . . 561
40.1.5 Charge distribution of fragments. . . . . . . . . . . . . 562
40.1.6 Kinetic energy distribution of fragments. . . . . . . . . 562
40.1.7 Calculation of the fragment excitation energies. . . . . 562

41 INCL++: the Liege Intranuclear Cascade model 564

41.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564
41.1.1 Suitable application fields . . . . . . . . . . . . . . . . 565
41.2 Generalities of the INCL++ cascade . . . . . . . . . . . . . . . 566
41.2.1 Model limits . . . . . . . . . . . . . . . . . . . . . . . . 566
41.3 Physics ingredients . . . . . . . . . . . . . . . . . . . . . . . . 567
41.3.1 Emission of composite particles . . . . . . . . . . . . . 568
41.3.2 Cascade stopping time . . . . . . . . . . . . . . . . . . 568
41.3.3 Conservation laws . . . . . . . . . . . . . . . . . . . . . 569
41.3.4 Initialisation of composite projectiles . . . . . . . . . . 569
41.3.5 De-excitation phase . . . . . . . . . . . . . . . . . . . . 569
41.4 Physics performance . . . . . . . . . . . . . . . . . . . . . . . 569
41.5 Status of this document . . . . . . . . . . . . . . . . . . . . . 571
42 ABLA V3 evaporation/fission model 574
42.1 Level densities . . . . . . . . . . . . . . . . . . . . . . . . . . . 575
42.2 Fission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575
42.3 External data file required . . . . . . . . . . . . . . . . . . . . 576
42.4 How to use ABLA V3 . . . . . . . . . . . . . . . . . . . . . . 576
42.5 Status of this document . . . . . . . . . . . . . . . . . . . . . 576

43 Low Energy Neutron Interactions 578

43.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 578
43.2 Physics and Verification . . . . . . . . . . . . . . . . . . . . . 578
43.2.1 Inclusive Cross-sections . . . . . . . . . . . . . . . . . . 578
43.2.2 Elastic Scattering . . . . . . . . . . . . . . . . . . . . . 579
43.2.3 Radiative Capture . . . . . . . . . . . . . . . . . . . . 581
43.2.4 Fission . . . . . . . . . . . . . . . . . . . . . . . . . . . 581
43.2.5 Inelastic Scattering . . . . . . . . . . . . . . . . . . . . 584
43.3 Neutron Data Library (G4NDL) Format . . . . . . . . . . . . 586
43.3.1 Cross Section . . . . . . . . . . . . . . . . . . . . . . . 586
43.3.2 Final State . . . . . . . . . . . . . . . . . . . . . . . . 587
43.3.3 Thermal Scattering Cross Section . . . . . . . . . . . . 587
43.3.4 Coherent Final State . . . . . . . . . . . . . . . . . . . 588
43.3.5 Incoherent Final State . . . . . . . . . . . . . . . . . . 589
43.3.6 Inelastic Final State . . . . . . . . . . . . . . . . . . . 590
43.4 High Precision Models and Low Energy Parameterized Models 592
43.5 Summary and Important Remark . . . . . . . . . . . . . . . . 592
43.6 Status of this document . . . . . . . . . . . . . . . . . . . . . 593

44 Radioactive Decay 595

44.1 The Radioactive Decay Module . . . . . . . . . . . . . . . . . 595
44.2 Alpha Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . 595
44.3 Beta Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . 596
44.4 Electron Capture . . . . . . . . . . . . . . . . . . . . . . . . . 596
44.5 Recoil Nucleus Correction . . . . . . . . . . . . . . . . . . . . 597
44.6 Biasing Methods . . . . . . . . . . . . . . . . . . . . . . . . . 597
44.7 Status of this document . . . . . . . . . . . . . . . . . . . . . 598

V Gamma- and Lepto-Nuclear Interactions 600

45 Introduction 601
45.1 Status of this document . . . . . . . . . . . . . . . . . . . . . 601
46 Cross Sections in Photonuclear/Electronuclear Reactions 602
46.1 Approximation of Photonuclear Cross Sections. . . . . . . . . 602
46.2 Electronuclear Cross Sections and Reactions . . . . . . . . . . 605
46.3 Common Notation for Electronuclear Reactions . . . . . . . . 605
46.4 Status of this document . . . . . . . . . . . . . . . . . . . . . 612

47 Gamma-nuclear Interactions 613

47.1 Process and Cross Section . . . . . . . . . . . . . . . . . . . . 613
47.2 Final State Generation . . . . . . . . . . . . . . . . . . . . . . 613
47.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 614

48 Electro-nuclear Interactions 615

48.1 Process and Cross Section . . . . . . . . . . . . . . . . . . . . 615
48.2 Final State Generation . . . . . . . . . . . . . . . . . . . . . . 615
48.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 616

49 Muon-nuclear Interactions 617

49.1 Process and Cross Section . . . . . . . . . . . . . . . . . . . . 617
49.2 Final State Generation . . . . . . . . . . . . . . . . . . . . . . 617
49.3 Status of this document . . . . . . . . . . . . . . . . . . . . . 618

0
 Part I

Introduction

1
Chapter 1

Introduction

1.1 Scope of This Manual

The Physics Reference Manual provides detailed explanations of the physics
implemented in the Geant4 toolkit. The manuals purpose is threefold:

to present the theoretical formulation, model, or parameterization of

the physics interactions included in Geant4,

to describe the probability of the occurrence of an interaction and the

sampling mechanisms required to simulate it, and

to serve as a reference for toolkit users and developers who wish to

consult the underlying physics of an interaction.

This manual does not discuss code implementation or how to use the
implemented physics interactions in a simulation. These topics are discussed
in the Users Guide for Application Developers. Details of the object-oriented
design and functionality of the Geant4 toolkit are given in the Users Guide
for Toolkit Developers. The Installation Guide for Setting up Geant4 in
Your Computing Environment describes how to get the Geant4 code, install
it, and run it.

1.2 Definition of Terms

Several terms used throughout the Physics Reference Manual have specific
meaning within Geant4, but are not well-defined in general usage. The defi-
nitions of these terms are given here.

2
 process - a C++ class which describes how and when a specific kind
of physical interaction takes place along a particle track. A given par-
ticle type typically has several processes assigned to it. Occaisionally
process refers to the interaction which the process class describes.

model - a C++ class whose methods implement the details of an in-

teraction, such as its kinematics. One or more models may be assigned
to each process. In sections discussing the theory of an interaction,
model may refer to the formulae or parameterization on which the
model class is based.

Geant3 - a physics simulation tool written in Fortran, and the prede-

cessor of Geant4. Although many references are made to Geant3, no
knowledge of it is required to understand this manual.

1.3 Status of this document

4.12.01 created by D.H. Wright

3
Chapter 2

Monte Carlo Methods

The Geant4 toolkit uses a combination of the composition and rejection

Monte Carlo methods. Only the basic formalism of these methods is outlined
here. For a complete account of the Monte Carlo methods, the interested user
is referred to the publications of Butcher and Messel, Messel and Crawford,
or Ford and Nelson [1, 2, 3].
Suppose we wish to sample x in the interval [x1 , x2 ] from the distribution
f (x) and the normalised probability density function can be written as :
n
X
f (x) = Ni fi (x)gi (x) (2.1)
i=1

where Ni > 0, fi (x) are normalised density functions on [x1 , x2 ] , and 0

gi (x) 1.
According to this method, x can sampled in the following way:

1. select a random integer i {1, 2, n} with probability proportional

to Ni

2. select a value x0 from the distribution fi (x)

3. calculate gi (x0 ) and accept x = x0 with probability gi (x0 );

4. if x0 is rejected restart from step 1.

It can be shown that

P this scheme is correct and the mean number of tries to
accept a value is i Ni .
In practice, a good method of sampling from the distribution f (x) has the
following properties:

all the subdistributions fi (x) can be sampled easily;

4
 the rejection functions gi (x) can be evaluated easily/quickly;

the mean number of tries is not too large.

Thus the different possible decompositions of the distribution f (x) are not
equivalent from the practical point of view (e.g. they can be very different
in computational speed) and it can be useful to optimise the decomposition.
A remark of practical importance : if our distribution is not normalised
Z x2
f (x)dx = C > 0
x1

the methodP can be used in the same manner; the mean number of tries in
this case is i Ni /C.

2.1 Status of this document

18.01.02 created by [Link].

Bibliography
[1] J.C. Butcher and H. Messel. Nucl. Phys. 20 15 (1960)

[2] H. Messel and D. Crawford. Electron-Photon shower distribution, Perg-

amon Press (1970)

[3] R. Ford and W. Nelson. SLAC-265, UC-32 (1985)

[4] Particle Data Group. Rev. of Particle Properties. Eur. Phys. J. C15.
(2000) 1. [Link]

5
Chapter 3

Particle Transport

6
3.1 Transportation
The transportation process is responsible for determining the geometrical
limits of a step. It calculates the length of step with which a track will cross
into another volume. When the track actually arrives at a boundary, the
transportation process locates the next volume that it enters.
If the particle is charged and there is an electromagnetic (or potentially
other) field, it is responsible for propagating the particle in this field. It does
this according to an equation of motion. This equation can be provided by
Geant4, for the case a magnetic or EM field, or can be provided by the user
for other fields.
The transportation updates the time of flight of a particle, utilising its
initial velocity.
Some additional details on motion in fields:
In order to intersect the model Geant4 geometry of a detector or setup,
the curved trajectory followed by a charged particle is split into chords seg-
ments. A chord is a straight line segment between two trajectory points.
Chords are created utilizing a criterion for the maximum estimated distance
between a curve point and the chord. This distance is also known as the
sagitta.
The equations of motions are solved utilising Runge Kutta methods.
Runge Kutta methods of different can be utilised for fields depending on the
numerical method utilised for approximating the field. Specialised methods
for near-constant magnetic fields are under development.

3.1.1 Status of This Document

17.11.11 minor revisions by V. Ivanchenko

7
3.2 True Step Length
Geant4 simulation of particle transport is performed step by step [1]. A
true step length for a next physics interaction is randomly sampled using the
mean free path of the interaction or by various step limitations established by
different Geant4 components. The smallest step limit defines the new true
step length.

3.2.1 The Interaction Length or Mean Free Path

Computation of mean free path of a particle in a media is performed in
Geant4 using cross section of a particular physics process and density of
atoms. In a simple material the number of atoms per volume is:
N
n=
A
where:
N Avogadros number
density of the medium
A mass of a mole
In a compound material the number of atoms per volume of the ith ele-
ment is:
N wi
ni =
Ai
where:
wi proportion by mass of the ith element
Ai mass of a mole of the ith element
The mean free path of a process, , also called the interaction length,
can be given in terms of the total cross section :
!1
X
(E) = [ni (Zi , E)]
i
P
where (Z, E) is the total cross section per atom of the process and i runs
overPall elements composing the material.
[ni (Zi , E)] is also called the macroscopic cross section. The mean free
i
path is the inverse of the macroscopic cross section.
Cross sections per atom and mean free path values may be tabulated during
initialisation.

8
3.2.2 Determination of the Interaction Point
The mean free path, , of a particle for a given process depends on the
medium and cannot be used directly to sample the probability of an inter-
action in a heterogeneous detector. The number of mean free paths which a
particle travels is:
Z x2
dx
n = , (3.1)
x1 (x)

which is independent of the material traversed. If nr is a random variable

denoting the number of mean free paths from a given point to the point of
interaction, it can be shown that nr has the distribution function:

P (nr < n ) = 1 en (3.2)

The total number of mean free paths the particle travels before reaching the
interaction point, n , is sampled at the beginning of the trajectory as:

n = log () (3.3)

where is a random number uniformly distributed in the range (0, 1). n is

updated after each step x according the formula:
x
n = n (3.4)
(x)

until the step originating from s(x) = n (x) is the shortest and this trig-
gers the specific process.

3.2.3 Step Limitations

The short description given above is the differential approach to particle
transport, which is used in the most popular simulation codes EGS and
Geant3. In this approach besides the other (discrete) processes the contin-
uous energy loss imposes a limit on the step-size too [2], because the cross
section of different processes depend of the energy of the particle. Then it
is assumed that the step is small enough so that the particle cross sections
remain approximately constant during the step. In principle one must use
very small steps in order to insure an accurate simulation, but computing
time increases as the step-size decreases. A good compromise depends on
required accuracy of a concrete simulation. For electromagnetic physics the

9
problem is reduced using integral approach, which is described below in sub-
chapter 7.3. However, this only provides effectively correct cross sections but
step limitation is needed also for more precise tracking. Thus, in Geant4 any
process may establish additional step limitation, the most important limits
see below in sub-chapters 7.1.3 and 6.1.6).

3.2.4 Updating the Particle Time

The laboratory time of a particle should be updated after each step:
1 1
tlab = 0.5x( + ), (3.5)
v1 v2
where x is a true step length traveled by the particle, v1 and v2 are particle
velocities at the beginning and at the end of the step correspondingly.

3.2.5 Status of This Document

09.10.98 created by L. Urban.
27.07.01 minor revisions by M. Maire
01.12.03 integral method subsection added by V. Ivanchenko
12.08.04 splitted and partly moved in introduction by M. Maire
25.12.06 minor revision by V. Ivanchenko
15.12.08 minor revision by J. Apostolakis
08.12.10 revisions by V. Ivanchenko
17.11.11 moved to transportation chapter by V. Ivanchenko

Bibliography
[1] S. Agostinelli et al., Geant4 a simulation toolkit Nucl. Instr. Meth.
A506 (2003) 250.

[2] J. Apostolakis et al., Geometry and physics of the Geant4 toolkit for high
and medium energy applications. Rad. Phys. Chem. 78 (2009) 859.

10
 Part II

Particle Decay

11
Chapter 4

Decay

The decay of particles in flight and at rest is simulated by the G4Decay class.

4.1 Mean Free Path for Decay in Flight

The mean free path is calculated for each step using

= c

where is the lifetime of the particle and

1
=p .
1 2

and are calculated using the momentum at the beginning of the step.
The decay time in the rest frame of the particle (proper time) is then sampled
and converted to a decay length using .

4.2 Branching Ratios and Decay Channels

G4Decay selects a decay mode for the particle according to branching ratios
defined in the G4DecayTable class, which is a member of the G4ParticleDefinition
class. Each mode is implemented as a class derived from G4VDecayChannel
and is responsible for generating the secondaries and the kinematics of the
decay. In a given decay channel the daughter particle momenta are calcu-
lated in the rest frame of the parent and then boosted into the laboratory
frame. Polarization is not currently taken into account for either the parent
or its daughters.

12
 A large number of specific decay channels may be required to simulate
an experiment, ranging from two-body to many-body decays and V A to
semi-leptonic decays. Most of these are covered by the five decay channel
classes provided by Geant4:
G4PhaseSpaceDecayChannel : phase space decay
G4DalitzDecayChannel : dalitz decay
G4MuonDecayChannel : muon decay
G4TauLeptonicDecayChannel : tau leptonic decay
G4KL3DecayChannel : semi-leptonic decays of kaon .

4.2.1 G4PhaseSpaceDecayChannel
The majority of decays in Geant4 are implemented using the G4PhaseSpaceDecayChannel
class. It simulates phase space decays with isotropic angular distributions in
the center-of-mass system. Three private methods of G4PhaseSpaceDecayChannel
are provided to handle two-, three- and N-body decays:
TwoBodyDecayIt()
ThreeBodyDecayIt()
ManyBodyDecayIt()
Some examples of decays handled by this class are:

0 ,

and

K 0L 0+.

4.2.2 G4DalitzDecayChannel
The Dalitz decay

0 + e+ + e

and other Dalitz-like decays, such as

K 0 L + e+ + e

and

K 0 L + + +

13
are simulated by the G4DalitzDecayChannel class. In general, it handles any
decay of the form

P 0 + l+ + l ,

where P 0 is a spin-0 meson of mass M and l are leptons of mass m. The

angular distribution of the is isotropic in the center-of-mass system of the
parent particle and the leptons are generated isotropically and back-to-back
in their center-of-mass frame. The magnitude of the leptons momentum is
sampled from the distribution function

r
t 3 2m2 4m2
f (t) = (1 2 ) (1 + ) 1 ,
M t t
where t is the square of the sum of the leptons energy in their center-of-mass
frame.

4.2.3 Muon Decay

G4MuonDecayChannel simulates muon decay according to V A theory. The
electron energy is sampled from the following distribution:

GF 2 m 5 2
d = 2 (3 2)
192 3
where: : decay rate
: = Ee /Emax
Ee : electron energy
Emax : maximum electron energy = m /2

The magnitudes of the two neutrino momenta are also sampled from the
V A distribution and constrained by energy conservation. The direction of
the electron neutrino is sampled using

cos() = 1 2/Ee 2/Ee + 2/Ee /Ee

and the muon anti-neutrino momentum is chosen to conserve momentum.

Currently, neither the polarization of the muon nor the electron is considered
in this class.

14
4.2.4 Leptonic Tau Decay
G4TauLeptonicDecayChannel simulates leptonic tau decays according to V
A theory. This class is valid for both
e + + e
and
+ +
modes.
The energy spectrum is calculated without neglecting lepton mass as
follows:
GF 2 m 3
d = pl El (3El m 2 4El 2 m 2m ml 2 )
24 3
where: : decay rate
El : daughter lepton energy (total energy)
pl : daughter lepton momentum
ml : daughter lepton mass

As in the case of muon decay, the energies of the two neutrinos are not
sampled from their V A spectra, but are calculated so that energy and
momentum are conserved. Polarization of the and final state leptons is not
taken into account in this class.

4.2.5 Kaon Decay

The class G4KL3DecayChannel simulates the following four semi-leptonic de-
cay modes of the kaon:

K e3 : K 0 + e +
K 3 : K 0 + +
K 0 e3 : KL0 + e +
K 0 3 : KL0 + +

Assuming that only the vector current contributes to K l decays, the

matrix element can be described by using two dimensionless form factors, f+
and f , which depend only on the momentum transfer t = (PK P )2 .
The Dalitz plot density used in this class is as follows [1]:

(E , E ) f+2 (t)[A + B (t) + C (t)2 ]

15
where: A = mK (2E E mK E ) + m 2 ( 41 E E )
B = m 2 (E 12 E )
C = 41 m 2 E
E = E max E

Here (t) is the ratio of the two form factors

(t) = f (t)/f+ (t).

f+ (t) is assumed to depend linearly on t, i.e.

f+ (t) = f+ (0)[1 + + (t/m 2 )]

and f (t) is assumed to be constant due to time reversal invariance.

Two parameters, + and (0) are then used for describing the Dalitz plot
density in this class. The values of these parameters are taken to be the
world average values given by the Particle Data Group [2].

4.3 Status of this document

05.07.12 updated muon decay section - D.H. Wright
10.04.02 re-written by D.H. Wright
02.04.02 editing by Hisaya Kurashige
14.11.01 editing by Hisaya Kurashige

Bibliography
[1] L.M. Chounet, J.M. Gaillard, and M.K. Gaillard, Phys. Reports 4C, 199
(1972).

[2] Review of Particle Physics The European Physical Journal C, 15 (2000).

16
 Part III

Electromagnetic Interactions

17
Chapter 5

Gamma Incident

18
5.1 Introduction
All processes of gamma interaction with media in Geant4 are happen at the
end of the step, so these interactions are discrete and corresponding processes
are following G4V DiscreteP rocess interface.

5.1.1 General Interfaces

There are a number of similar functions for discrete electromagnetic pro-
cesses and for electromagnetic (EM) packages an additional base classes were
designed to provide common computations [1]. Common calculations for
discrete EM processes are performed in the class G4V EmP rocess. Derived
classes (5.1) are concrete processes providing initialisation. The physics mod-
els are implemented using the G4V EmM odel interface. Each process may
have one or many models defined to be active over a given energy range
and set of G4Regions. Models are implementing computation of energy loss,
cross section and sampling of final state. The list of EM processes and models
for gamma incident is shown in Table 5.1.

5.1.2 Status of This Document

06.12.07 created by V. Ivanchenko
11.12.08 extended list of models by V. Ivanchenko
08.12.10 cleaned up by V. Ivanchenko
20.11.11 updated list of processes/models by V. Ivanchenko
29.11.13 updated list of processes/models by V. Ivanchenko

Bibliography
[1] J. Apostolakis et al., Geometry and physics of the Geant4 toolkit for high
an dmedium energy applications. Rad. Phys. Chem. 78 (2009) 859.

19
 Table 5.1: List of process and model classes for gamma.
EM process EM model Ref.
G4PhotoElectricEffect G4PEEffectFluoModel 5.2
G4LivermorePhotoElectricModel 9.7
G4LivermorePolarizedPhotoElectricModel
G4PenelopePhotoElectricModel 10.1.5
G4PolarizedPhotoElectricEffect G4PolarizedPEEffectModel 17.1
G4ComptonScattering G4KleinNishinaCompton 5.3
G4KleinNishinaModel 5.3
G4LivermoreComptonModel 9.2
G4LivermoreComptonModelRC
G4LivermorePolarizedComptonModel 9.3
G4LowEPComptonModel 11.1
G4PenelopeComptonModel 10.1.2
G4PolarizedCompton G4PolarizedComptonModel 17.1
G4GammaConversion G4BetheHeitlerModel 5.4
G4PairProductionRelModel
G4LivermoreGammaConversionModel 9.5
G4BoldyshevTripletModel 9.6
G4LivermoreNuclearGammaConversionModel
G4LivermorePolarizedGammaConversionModel
G4PenelopeGammaConvertion 10.1.4
G4PolarizedGammaConversion G4PolarizedGammaConversionModel 17.1
G4RayleighScattering G4LivermoreRayleighModel 9.4
G4LivermorePolarizedRayleighModel
G4PenelopeRayleighModel 10.1.3
G4GammaConversionToMuons 5.5

20
5.2 PhotoElectric effect
The photoelectric effect is the ejection of an electron from a material af-
ter a photon has been absorbed by that material. In the standard model
G4PEEffectFluoModel it is simulated by using a parameterized photon ab-
sorption cross section to determine the mean free path, atomic shell data to
determine the energy of the ejected electron, and the K-shell angular distri-
bution to sample the direction of the electron.

5.2.1 Cross Section

The parameterization of the photoabsorption cross section proposed by Biggs
et al. [1] was used :

a(Z, E ) b(Z, E ) c(Z, E ) d(Z, E )

(Z, E ) = + + + (5.1)
E E2 E3 E4

Using the least-squares method, a separate fit of each of the coefficients

a, b, c, d to the experimental data was performed in several energy intervals
[2]. As a rule, the boundaries of these intervals were equal to the correspond-
ing photoabsorption edges. The cross section (and correspondingly mean free
path) are discontinuous and must be computed on the fly from the formula
5.1. Coefficients are defined to each Sandia table energy interval.
If photon energy is below the lowest Sandia energy for the material the
cross section is computed for this lowest energy, so gamma is absorbed by
photoabsorption at any energy. This approach is implemented coherently for
models of photoelectric effect of Geant4. As a result, any media become not
transparant for low-energy gammas.

5.2.2 Final State

Choosing an Element
The binding energies of the shells depend on the atomic number Z of the ma-
terial. In compound materials the ith element is chosen randomly according
to the probability:

nati (Zi , E )
P rob(Zi , E ) = P .
i [nati i (E )]

21
Shell
A quantum can be absorbed if E > Bshell where the shell energies are taken
from G4AtomicShells data: the closest available atomic shell is chosen. The
photoelectron is emitted with kinetic energy :

Tphotoelectron = E Bshell (Zi ) (5.2)

Theta Distribution of the Photoelectron

The polar angle of the photoelectron is sampled from the Sauter-Gavrila
distribution (for K-shell) [3], which is correct only to zero order in Z :

sin2
 
d 1
1 + ( 1)( 2)(1 cos ) (5.3)
d(cos ) (1 cos )4 2
where and are the Lorentz factors of the photoelectron.
cos is sampled from the probability density function :
1 2 1 (1 2r) +
f (cos ) = = cos = (5.4)
2 (1 cos )2 (1 2r) + 1
The rejection function is :

1 cos2
g(cos ) = [1 + b(1 cos )] (5.5)
(1 cos )2

with b = ( 1)( 2)/2

It can be shown that g(cos ) is positive cos [1, +1], and can be
majored by :

gsup = 2 [1 + b(1 )] if ]1, 2] (5.6)

= 2 [1 + b(1 + )] if > 2

The efficiency of this method is 50% if < 2, 25% if [2, 3].

5.2.3 Relaxation
Atomic relaxations can be sampled using the de-excitation module of the low-
energy sub-package 14.1. For that atomic de-excitation option should be acti-
vated. In the physics list sub-library this activation is done automatically for
G4EmLivermorePhysics, G4EmPenelopePhysics, G4EmStandardPhysics option3
and G4EmStandardPhysics option4. For other standard physics constructors
the de-excitation module is already added but is disabled. The simulation of

22
fluorescence and Auger electron emmision may be enabled for all geometry
via UI commands:

/process/em/fluo true
/process/em/auger true

There is a possiblity to enable atomic deexcitation only for G4Region by

its name:

/process/em/deexcitation myregion true true false

where three boolean arguments enable/disable fluorescence, Auger electron

production and PIXE (deexcitation induced by ionisation).

5.2.4 Status of this document

09.10.98 created by M. Maire
08.01.02 updated by M. Maire
22.04.02 re-worded by D.H. Wright
02.05.02 modifs in total cross section and final state (M. Maire)
15.11.02 introduction added by D.H. Wright
08.12.10 revision by V. Ivanchenko
20.11.11 revision by V. Ivanchenko
20.12.12 revision by V. Ivanchenko

Bibliography
[1] Biggs F., and Lighthill R., Preprint Sandia Laboratory, SAND 87-0070
(1990)

[2] Grichine V.M., Kostin A.P., Kotelnikov S.K. et al., Bulletin of the Lebe-
dev Institute no. 2-3, 34 (1994).

[3] Gavrila M. [Link]. 113, 514 (1959).

23
5.3 Compton scattering
The Compton scattering is an inelastic gamma scattering on atom with the
ejection of an electron. In the standard sub-package two model G4KleinNishinaCompton
and G4KleinNishinaModel are available. The first model is the fastest, in the
second model atomic shell effects are taken into account.

5.3.1 Cross Section

When simulating the Compton scattering of a photon from an atomic elec-
tron, an empirical cross section formula is used, which reproduces the cross
section data down to 10 keV:
log(1 + 2X) P2 (Z) + P3 (Z)X + P4 (Z)X 2
 
(Z, E ) = P1 (Z) + . (5.7)
X 1 + aX + bX 2 + cX 3

Z = atomic number of the medium

E = energy of the photon
X = E /mc2
m = electron mass
Pi (Z) = Z(di + ei Z + fi Z 2 ).

The values of the parameters can be found within the method which computes
the cross section per atom. A fit of the parameters was made to over 511
data points [1, 2] chosen from the intervals

1 Z 100

E [10 keV, 100 GeV].

The accuracy of the fit was estimated to be


10% for E 10 keV 20 keV
=
5 6% for E > 20 keV

To avoid sampling problems in the Compton process the cross section is set
to zero at low-energy limit of cross section table, which is 100eV in majority
of EM Phyiscs Lists.

24
5.3.2 Sampling the Final State
The Klein-Nishina differential cross section per atom is [3]:
2
sin2
  
d 2 me c 1
= re Z + 1 (5.8)
d E0 1 + 2

where re = classical electron radius Assuming an elastic col-

me c2 = electron mass
E0 = energy of the incident photon
E1 = energy of the scattered photon
= E1 /E0 .
lision, the scattering angle is defined by the Compton formula:

me c2
E1 = E0 . (5.9)
me c2 + E0 (1 cos )

Sampling the Photon Energy

The value of corresponding to the minimum photon energy (backward scat-
tering) is given by
me c2
0 = , (5.10)
me c2 + 2E0
hence [0 , 1]. Using the combined composition and rejection Monte Carlo
methods described in [4, 5, 6] one may set

sin2
  
1
() + 1 = f ()g() = [1 f1 () + 2 f2 ()]g(), (5.11)
1 + 2

1 = ln(1/0 ) ; f1 () = 1/(1 )
2
2 = (1 0 )/2 ; f2 () = /2 .
f1 and f2 are probability density functions defined on the interval [0 , 1], and
 
2
g() = 1 sin
1 + 2

is the rejection function [0 , 1] = 0 < g() 1. Given a set of

3 random numbers r, r , r uniformly distributed on the interval [0,1], the
sampling procedure for is the following:

1. decide whether to sample from f1 () or f2 ():

if r < 1 /(1 + 2 ) select f1 (), otherwise select f2 ()

25
 2. sample from the distributions corresponding to f1 or f2 :

for f1 : = r0 ( exp(r 1 ))
for f2 : = 0 + (1 20 )r
2 2

3. calculate sin2 = t(2 t) where t (1 cos ) = me c2 (1 )/(E0 )

4. test the rejection function:

if g() r accept , otherwise go to step 1.

Compute the Final State Kinematics

After the successful sampling of , the polar angles of the scattered photon
with respect to the direction of the parent photon are generated. The az-
imuthal angle, , is generated isotropically and is as defined in the previous

section. The momentum vector of the scattered photon, P1 , is then trans-
formed into the World coordinate system. The kinetic energy and momentum
of the recoil electron are then

Tel = E0 E1


Pel = P0 P1 .

Doppler broading of final electron momentum due to electron motion is

implemented only in G4KleinNishinaModel. For that emphirical electron
density profile function is used.

5.3.3 Atomic shell effects

The differential cross-section described above is valid only for those collisions
in which the energy of the recoil electron is large compared to its binding
energy (which is ignored). In the alternative model (G4KleinNishinaModel)
atomic shell effects are taken into account. For that a sampling of a shell is
performed with the weight proportional to number of shell electrons. Electron
energy distribution function is approximated via simplified form

F (T ) = exp (T /Eb )/Eb , (5.12)

where Eb is shell bound energy, T - kinetic energy of the electron.

The value T is sampled and scattering is sampled in the rest frame of
the electron according the algorithm described in the previous sub-chapter.
After sampling an inverse Lorentz transformation to the laboratory frame is
performed. Potential energy (Eb + T ) is subtracted from the scattered elec-
tron kinetic energy. If final electron energy become negative then sampling is

26
repeated. Atomic relaxation are sampled if deexcitation module is enabled.
Enabling of atomic relaxation for Compton scattering is performed in the
same way as for photoelectric effect 5.2.3.

5.3.4 Status of This Document

09.10.98 created by M. Maire
14.01.02 minor revision by M. Maire
22.04.02 reworded by D.H. Wright
18.03.04 include references for total cross section (M. Maire)
10.12.10 revised by V. Ivanchenko
20.11.12 revised by V. Ivanchenko
29.11.13 revised by V. Ivanchenko

Bibliography
[1] Hubbell, Gimm and Overbo. J. Phys. Chem. Ref. Data 9 (1980) 1023.

[2] H. Storm and H.I. Israel Nucl. Data Tables A7 (1970) 565.

[3] O. Klein and Y. Nishina. Z. Physik 52 (1929) 853.

[4] J.C. Butcher and H. Messel. Nucl. Phys. 20 (1960) 15.

[5] H. Messel and D. Crawford. Electron-Photon shower distribution, Perg-

amon Press (1970)

[6] R. Ford and W. Nelson. SLAC-265, UC-32 (1985).

[7] B. Rossi. High energy particles, Prentice-Hall 77-79 (1952)

27
5.4 Gamma Conversion into e+e Pair
In the standard sub-package two models are available. The first model is
implemented in the class G4BetheHeitlerModel, it was derived from Geant3
and is applicable below 100GeV . In the second (G4PairProductionRelModel)
Landau-Pomenrachuk-Migdal (LPM) effect is taken into account and this
model can be applied for high energy gammas (above 100M eV ).

5.4.1 Cross Section

According [1], [2] the total cross-section per atom for the conversion of a
gamma into an (e+ , e ) pair has been parameterized as
 
F3 (X)
(Z, E ) = Z(Z + 1) F1 (X) + F2 (X) Z + , (5.13)
Z

where E is the incident gamma energy and X = ln(E /me c2 ) . The functions
Fn are given by

F1 (X) = a0 + a1 X + a2 X 2 + a3 X 3 + a4 X 4 + a5 X 5 (5.14)
F2 (X) = b0 + b1 X + b2 X 2 + b3 X 3 + b4 X 4 + b5 X 5
F3 (X) = c0 + c1 X + c2 X 2 + c3 X 3 + c4 X 4 + c5 X 5 ,

with the parameters ai , bi , ci taken from a least-squares fit to the data [1].
Their values can be found in the function which computes formula 5.13.
This parameterization describes the data in the range

1 Z 100

and

E [1.5 MeV, 100 GeV].

The accuracy of the fit was estimated to be 5% with a mean value of

2.2%. Above 100 GeV the cross section is constant. Below Elow = 1.5 MeV
the extrapolation
2
E 2me c2

(E) = (Elow ) (5.15)
Elow 2me c2

is used.

28
 In a given material the mean free path, , for a photon to convert into
an (e+ , e ) pair is
!1
X
(E ) = nati (Zi , E ) (5.16)
i

where nati is the number of atoms per volume of the ith element of the
material.

Corrected Bethe-Heitler Cross Section

As written in [2], the Bethe-Heitler formula corrected for various effects is
  
d(Z, ) 2 2 2 F (Z)
= re Z[Z + (Z)] [ + (1 ) ] 1 (())
d 2
 
2 F (Z)
+ (1 ) 2 (()) (5.17)
3 2

where is the fine-structure constant and re the classical electron radius.

Here = E/E , E is the energy of the photon and E is the total energy
carried by one particle of the (e+ , e ) pair. The kinematical limits of are
therefore
me c2
= 0 1 0 . (5.18)
E

Screening Effect The screening variable, , is a function of

136 0
() = 1/3
, (5.19)
Z (1 )

and measures the impact parameter of the projectile. Two screening func-
tions are introduced in the Bethe-Heitler formula :

for 1 1 () = 20.867 3.242 + 0.625 2 (5.20)

2 () = 20.209 1.930 0.086 2
for > 1 1 () = 2 () = 21.12 4.184 ln( + 0.952).

Because the formula 5.17 is symmetric under the exchange (1 ), the

range of can be restricted to

[0 , 1/2]. (5.21)

29
Born Approximation The Bethe-Heitler formula is calculated with plane
waves, but Coulomb waves should be used instead. To correct for this, a
Coulomb correction function is introduced in the Bethe-Heitler formula :

for E < 50 MeV : F (z) = 8/3 ln Z (5.22)

for E 50 MeV : F (z) = 8/3 ln Z + 8fc (Z)

with

2 1
fc (Z) = (Z) (5.23)
1 + (Z)2
+0.20206 0.0369(Z)2 + 0.0083(Z)4 0.0020(Z)6 + .


It should be mentioned that, after these additions, the cross section becomes
negative if  
42.24 F (Z)
> max (1 ) = exp 0.952. (5.24)
8.368
This gives an additional constraint on :
r
1 1 min
max = 1 = 1 (5.25)
2 2 max
where  
1 136
min = = = 40 (5.26)
2 Z 1/3
has been introduced. Finally the range of becomes

[min = max(0 , 1 ), 1/2]. (5.27)

30
 ()

d max

d min

0 0 1 1/2 1

Gamma Conversion in the Electron Field The electron cloud gives an

additional contribution to pair creation, proportional to Z (instead of Z 2 ).
This is taken into account through the expression
ln(1440/Z 2/3 )
(Z) = . (5.28)
ln(183/Z 1/3 ) fc (Z)

Factorization of the Cross Section is sampled using the techniques of

composition+rejection, as treated in [3, 4, 5]. First, two auxiliary screening
functions should be introduced:
F1 () = 31 () 2 () F (Z)
3 1
F2 () = 1 () 2 () F (Z) (5.29)
2 2
It can be seen that F1 () and F2 () are decreasing functions of ,
[min , max ]. They reach their maximum for min = ( = 1/2) :
F10 = max F1 () = F1 (min )
F20 = max F2 () = F2 (min ). (5.30)
After some algebraic manipulations the formula 5.17 can be written :
 
d(Z, ) 2 2 1
= re Z[Z + (Z)] min
d 9 2
[N1 f1 () g1 () + N2 f2 () g2 ()] , (5.31)

31
where
 2
1 3
1 2 F1 ()
N1 = min F10 f1 () = 3 2
g1 () =
2 [ 1
2
min] F10
3 1 F2 ()
N2 = F20 f2 () = const = g2 () = .
2 [ 21 min ] F20
f1 () and f2 () are probability density functions on the interval [min , 1/2]
such that Z 1/2
fi () d = 1
min

, and g1 () and g2 () are valid rejection functions: 0 < gi () 1 .

5.4.2 Final State

The differential cross section depends on the atomic number Z of the material
in which the interaction occurs. In a compound material the element i in
which the interaction occurs is chosen randomly according to the probability

nati (Zi , E )
P rob(Zi , E ) = P . (5.32)
i [nati i (E )]

Sampling the Energy Given a triplet of uniformly distributed random

numbers (ra , rb , rc ) :

1. use ra to choose which decomposition term in 5.31 to use:

if ra < N1 /(N1 + N2 ) f1 () g1 () otherwise f2 () g2 () (5.33)

2. sample from f1 () or f2 () with rb :

   
1 1 1/3 1
= min rb or = min + min rb (5.34)
2 2 2

3. reject if g1 ()or g2 () < rc

note : below E = 2 MeV it is enough to sample uniformly on [0 , 1/2],

without rejection.

Charge The charge of each particle of the pair is fixed randomly.

32
Polar Angle of the Electron or Positron
The polar angle of the electron (or positron) is defined with respect to the
direction of the parent photon. The energy-angle distribution given by Tsai
[6] is quite complicated to sample and can be approximated by a density
function suggested by Urban [7] :

9a2
u [0, [ f (u) = [u exp(au) + d u exp(3au)] (5.35)
9+d
with
5 mc2
a= d = 27 and = u. (5.36)
8 E
A sampling of the distribution 5.35 requires a triplet of random numbers such
that
9 ln(r2 r3 ) ln(r2 r3 )
if r1 < u= otherwise u = . (5.37)
9+d a 3a
The azimuthal angle is generated isotropically. The e+ and e momenta are
assumed to be coplanar with the parent photon. This information, together
with energy conservation, is used to calculate the momentum vectors of the
(e+ , e ) pair and to rotate them to the global reference system.

5.4.3 Ultra-Relativistic Model

It is implemented in the class G4PairProductionRelModel and is configured
above 80GeV in all reference Physics lists. The cross section is computed
using direct integration of differential cross section [6] and not its parameter-
isation described in 5.4.1. LPM effect is taken into account in the same way
as for bremsstrahlung 8.2.2. Secondary generation algorithm is the same as
in the standard Bethe-Haitler model.

5.4.4 Status of This Document

12.01.02 created by [Link].
21.03.02 corrections in angular distribution (mma)
22.04.02 re-worded by D.H. Wright
10.12.10 revision by V. Ivanchenko
20.11.12 revision by V. Ivanchenko

33
Bibliography
[1] [Link], [Link], [Link] Jou. Phys. Chem. Ref. Data 9:1023
(1980)

[2] W. Heitler The Quantum Theory of Radiation, Oxford University Press

(1957)

[3] R. Ford and W. Nelson. SLAC-210, UC-32 (1978)

[4] J.C. Butcher and H. Messel. Nucl. Phys. 20 15 (1960)

[5] H. Messel and D. Crawford. Electron-Photon shower distribution, Perg-

amon Press (1970)

[6] Y. S. Tsai, Rev. Mod. Phys. 46 815 (1974), Y. S. Tsai, Rev. Mod. Phys.
49 421 (1977)

[7] [Link] in Geant3 writeup, section PHYS-211. Cern Program Library

(1993)

34
5.5 Gamma Conversion into + Pair
The class G4GammaConversionToMuons simulates the process of gamma
conversion into muon pairs. Given the photon energy and Z and A of the
material in which the photon converts, the probability for the conversions
to take place is calculated according to a parameterized total cross section.
Next, the sharing of the photon energy between the + and is deter-
mined. Finally, the directions of the muons are generated. Details of the
implementation are given below and can be also found in [1].

5.5.1 Cross Section and Energy Sharing

Muon pair production on atomic electrons, + e e + + + , has a
threshold of 2m (m + me )/me 43.9 GeV . Up to several hundred GeV
this process has a much lower cross section than the corresponding process
on the nucleus. At higher energies, the cross section on atomic electrons
represents a correction of 1/Z to the total cross section.
For the approximately elastic scattering considered here, momentum, but
no energy, is transferred to the nucleon. The photon energy is fully shared
by the two muons according to

E = E+ + E (5.38)

or in terms of energy fractions

E+ E
x+ = , x = , x+ + x = 1 .
E E
The differential cross section for electromagnetic pair creation of muons in
terms of the energy fractions of the muons is
 
d 2 2 4
= 4 Z rc 1 x+ x log(W ) , (5.39)
dx+ 3
where Z is the charge of the nucleus, rc is the classical radius of the particles
which are pair produced (here muons) and

1 + (Dn e 2) /m
W = W (5.40)
1 + B Z 1/3 e /me
where
B Z 1/3 m m2
W = = e = 1.6487 . . . .
Dn m e 2 E x+ x

35
 For hydrogen B = 202.4 Dn = 1.49
and for all other nuclei B = 183 Dn = 1.54 A0.27 . (5.41)

These formulae are obtained from the differential cross section for muon
bremsstrahlung [2] by means of crossing relations. The formulae take into
account the screening of the field of the nucleus by the atomic electrons in
the Thomas-Fermi model, as well as the finite size of the nucleus, which is
essential for the problem under consideration. The above parameterization
gives good results for E m . The fact that it is approximate close
to threshold is of little practical importance. Close to threshold, the cross
section is small and the few low energy muons produced will not travel very
far. The cross section calculated from Eq. (5.39) is positive for E > 4m
and
s s
1 1 m 1 1 m
xmin x xmax with xmin = xmax = + ,
2 4 E 2 4 E
(5.42)
except for very asymmetric pair-production, close to threshold, which can
easily be taken care of by explicitly setting = 0 whenever < 0.
Note that the differential cross section is symmetric in x+ and x and
that
x+ x = x x2
where x stands for either x+ or x . By defining a constant

0 = 4 Z 2 rc2 log(W ) (5.43)

the differential cross section Eq. (5.39) can be rewritten as a normalized and
symmetric as function of x:
 
1 d 4 2 log W
= 1 (x x ) . (5.44)
0 dx 3 log W

This is shown in Fig. 5.1 for several elements and a wide range of photon
energies. The asymptotic differential cross section for E

1 d 4
= 1 (x x2 )
0 dx 3
is also shown.

36
 1
E
H Z=1 A=1.00794
0.9 Be Z=4 A=9.01218
Pb Z=82 A=207.2
100 TeV

0.8

0.7 10 TeV

1 TeV
0.6

0.5
0 dx
d

0.4 100 GeV

0.3

10 GeV
0.2

0.1
E = 1 GeV

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Figure 5.1: Normalized differential cross section for pair production as a

function of x, the energy fraction of the photon energy carried by one of
the leptons in the pair. The function is shown for three different elements,
hydrogen, beryllium and lead, and for a wide range of photon energies.

37
5.5.2 Parameterization of the Total Cross Section
The total cross section is obtained by integration of the differential cross
section Eq. (5.39), that is
Z xmax Z xmax  
d 2 2 4
tot (E ) = dx+ = 4 Z rc 1 x+ x log(W ) dx+ .
xmin dx+ xmin 3
(5.45)
W is a function of (x+ , E ) and (Z, A) of the element (see Eq. (5.40)). Nu-
merical values of W are given in Table 5.2.

Table 5.2: Numerical values of W for x+ = 0.5 for different elements.

E W for H W for Be W for Cu W for Pb
GeV
1 2.11 1.594 1.3505 5.212
10 19.4 10.85 6.803 43.53
100 191.5 102.3 60.10 332.7
1000 1803 919.3 493.3 1476.1
10000 11427 4671 1824 1028.1
28087 8549 2607 1339.8

Values of the total cross section obtained by numerical integration are

listed in Table 5.3 for four different elements. Units are in barn , where
1 barn = 1034 m2 .

Table 5.3: Numerical values for the total cross section

E tot , H tot , Be tot , Cu tot , Pb
GeV barn barn barn barn
1 0.01559 0.1515 5.047 30.22
10 0.09720 1.209 49.56 334.6
100 0.1921 2.660 121.7 886.4
1000 0.2873 4.155 197.6 1476
10000 0.3715 5.392 253.7 1880
0.4319 6.108 279.0 2042

Well above threshold, the total cross section rises about linearly in log(E )
with the slope
1
WM = (5.46)
4 Dn e m

38
 1

0.9

0.8

0.7

/
0.6

0.5

0.4

0.3

0.2 H Pb
0.1

0
2 3 4 5 6 7 8
1 10 10 10 10 10 10 10 10

E in GeV

Figure 5.2: Total cross section for the Bethe-Heitler process + as a

function of the photon energy E in hydrogen and lead, normalized to the
asymptotic cross section .

until it saturates due to screening at . Fig. 5.2 shows the normalized cross
section where
7
= 0 and 0 = 4 Z 2 rc2 log(W ) . (5.47)
9
Numerical values of WM are listed in Table 5.4.

Table 5.4: Numerical values of WM .

Element WM
1/GeV
H 0.963169
Be 0.514712
Cu 0.303763
Pb 0.220771

The total cross section can be parameterized as

28 Z 2 rc2
par = log(1 + WM Cf Eg ) , (5.48)
9
with  t
4m s

1/s
Eg = 1 Wsat + Es . (5.49)
E
and
W 1/3
4 e m2
Wsat = =BZ .
WM me

39
The threshold behavior in the cross section was found to be well approxi-
mated by t = 1.479 + 0.00799Dn and the saturation by s = 0.88. The
agreement at lower energies is improved using an empirical correction factor,
applied to the slope WM , of the form
  
Ec
Cf = 1 + 0.04 log 1 + ,
E
where  
4347.
Ec = 18. + GeV .
B Z 1/3
A comparison of the parameterized cross section with the numerical integra-
tion of the exact cross section shows that the accuracy of the parametrization
is better than 2%, as seen in Fig. 5.3.
1.02
Pb
1.01
/ par

Cu
1
Be
0.99 H
0.98
2 3 4 5 6 7 8
1 10 10 10 10 10 10 10 10
E in GeV

Figure 5.3: Ratio of numerically integrated and parametrized total cross

sections as a function of E for hydrogen, beryllium, copper and lead.

5.5.3 Multi-differential Cross Section and Angular Vari-

ables
The angular distributions are based on the multi-differential cross section for
lepton pair production in the field of the Coulomb center

d 4 Z 2 3 m2
= u+ u
dx+ du+ du d q4

u2+ + u2

2x+ x
(1 + u2+ ) (1 + u2 )
u2+ u2
  
2u+ u (1 2x+ x ) cos
+ . (5.50)
(1 + u2+ )2 (1 + u2 )2 (1 + u2+ ) (1 + u2 )
Here
E
u = , = , q 2 = qk2 + q
2
, (5.51)
m

40
where

qk2 = qmin
2
(1 + x u2+ + x+ u2 )2 ,
2 2
(u+ u )2 + 2 u+ u (1 cos ) .
 
q = m (5.52)

q 2 is the square of the momentum q transferred to the target and qk2 and q
2

are the squares of the components of the vector q, which are parallel and
perpendicular to the initial photon momentum, respectively. The minimum
momentum transfer is qmin = m2 /(2E x+ x ).

The muon vectors have the components

p+ = p+ ( sin + cos(0 + /2) , sin + sin(0 + /2) , cos + ) ,

p = p ( sin cos(0 /2) , sin sin(0 /2) , cos ) ,
q (5.53)
where p = E2 m2 . The initial photon direction is taken as the z-axis.
The cross section of Eq. (5.50) does not depend on 0 . Because of azimuthal
symmetry, 0 can simply be sampled at random in the interval (0, 2 ).
Eq. (5.50) is too complicated for efficient Monte Carlo generation. To
simplify, the cross section is rewritten to be symmetric in u+ , u using a
new variable u and small parameters , , where u = u /2 and = u .
When higher powers in small parameters are dropped, the differential cross
section in terms of u, , becomes

d 4 Z 2 3 m2
= 2 (5.54)
dx+ d d udu

qk2 + m2 ( 2 + 2 )
(1 u2 )2 2 (1 2x+ x )
   
2 1
2 x+ x + ,
(1 + u2 )2 (1 + u2 )4 (1 + u2 )2

where, in this approximation,

qk2 = qmin
2
(1 + u2 )2 .

For Monte Carlo generation, it is convenient to replace (, ) by the polar

coordinates (, ) with = cos and = sin . Integrating Eq. 5.54
over and using symbolically du2 where du2 = 2u du yields

4Z 2 3 3 1 x+ x x+ x (1 u2 )2
 
d
= .
dx+ d du2 m2 (qk2 /m2 + 2 )2 (1 + u2 )2 (1 + u2 )4
(5.55)

41
Integration with logarithmic accuracy over gives

Z1
3 d
 
d m
Z
= log . (5.56)
(qk2 /m2 + 2 )2 qk
qk /m

Within the logarithmic accuracy, log(m /qk ) can be replaced by log(m /qmin ),
so that
4 Z 2 3 1 x+ x x+ x (1 u2 )2
   
d m
= log . (5.57)
dx+ du2 m2 (1 + u2 )2 (1 + u2 )4 qmin

Making the substitution u2 = 1/t 1, du2 = dt /t2 gives

4 Z 2 3
 
d m
= [1 2 x+ x + 4 x+ x t (1 t)] log . (5.58)
dx+ dt m2 qmin

Atomic screening and the finite nuclear radius may be taken into account by
multiplying the differential cross section determined by Eq. (5.55) with the
factor
(Fa (q) Fn (q) )2 , (5.59)
where Fa and Fn are atomic and nuclear form factors. Please note that after
integrating Eq. 5.55 over , the q-dependence is lost.

5.5.4 Procedure for the Generation of + Pairs

Given the photon energy E and Z and A of the material in which the
converts, the probability for the conversions to take place is calculated ac-
cording to the parametrized total cross section Eq. (5.48). The next step,
determining how the photon energy is shared between the + and , is
done by generating x+ according to Eq. (5.39). The directions of the muons
are then generated via the auxilliary variables t, , . In more detail, the
final state is generated by the following five steps, in which R1,2,3,4,... are ran-
dom numbers with a flat distribution in the interval [0,1]. The generation
proceeds as follows.

1) Sampling of the positive muon energy E+ = x+ E .

This is done using the rejection technique. x+ is first sampled from a flat
distribution within kinematic limits using

x+ = xmin + R1 (xmax xmin )

42
and then brought to the shape of Eq. (5.39) by keeping all x+ which satisfy
 
4 log(W )
1 x+ x < R2 .
3 log(Wmax )

Here Wmax = W (x+ = 1/2) is the maximum value of W , obtained for sym-
metric pair production at x+ = 1/2. About 60% of the events are kept in this
step. Results of a Monte Carlo generation of x+ are illustrated in Fig. 5.4.
The shape of the histograms agrees with the differential cross section illus-
trated in Fig. 5.1.

30000
E

25000 10 GeV
100 GeV
20000
1000 GeV
15000

10000

5000

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x+

Figure 5.4: Histogram of generated x+ distributions for beryllium at three

different photon energies. The total number of entries at each energy is 106 .

2) Generate t(= 2 12 +1 ) .
The distribution in t is obtained from Eq.(5.58) as

1 2 x+ x + 4 x+ x t (1 t)
f1 (t) dt = dt , 0 < t 1. (5.60)
1 + C1 /t2

with form factors taken into account by

(0.35 A0.27 )2
C1 = . (5.61)
x+ x E /m

In the interval considered, the function f1 (t) will always be bounded from
above by
1 x+ x
max[f1 (t)] = .
1 + C1
For small x+ and large E , f1 (t) approaches unity, as shown in Fig. 5.5.

43
 1 1 E = 1 TeV 0.01
E = 10 GeV
0.1
0.1
0.8 0.8 0.25
0.25

0.6 0.5 0.6

0.5

f1(t)
f1(t)

0.4 0.4 x+

0.01
0.2 0.2
x+
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t t

Figure 5.5: The function f1 (t) at E = 10 GeV (left) and E = 1 TeV (right)
in beryllium for different values of x+ .

30000
E = 10 GeV
25000

20000
E = 1 TeV
15000

10000

5000

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t

Figure 5.6: Histograms of generated t distributions for E = 10 GeV (solid

line) and E = 100 GeV (dashed line) with 106 events each.

30000

25000

20000

15000
1 GeV
10000 10 GeV
100 GeV
5000 1000 GeV

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
/

Figure 5.7: Histograms of generated distributions for beryllium at four

different photon energies.

44
The Monte Carlo generation is done using the rejection technique. About
70% of the generated numbers are kept in this step. Generated t-distributions
are shown in Fig. 5.6.

3) Generate by the rejection technique using t generated in the previous

step for the frequency distribution
h i
f2 () = 12 x+ x +4 x+ x t (1t) (1+cos(2)) , 0 2 . (5.62)

The maximum of f2 () is

max[f2 ()] = 1 2 x+ x [1 4 t (1 t)] . (5.63)

Generated distributions in are shown in Fig. 5.7.

4) Generate .
The distribution in has the form

3 d
f3 () d = , 0 max , (5.64)
4 + C 2
where  
1.9 1
2max = 0.27 1 , (5.65)
A t
and
" 2  2 # 2
4 m me
C2 = + . (5.66)
x+ x 2E x+ x t 183 Z 1/3 m

The distribution is obtained by a direct transformation applied to uniform

random numbers Ri according to

= [C2 (exp( Ri ) 1)]1/4 , (5.67)

where
C2 + 4max
 
= log . (5.68)
C2
Generated distributions of are shown in Fig. 5.8
5) Calculate + , and from t, , with

E
r
1
= and u= 1. (5.69)
m t

45
 x 100
2000
1800
1600
1400
1200
1000 1 TeV
800
600 E = 10 GeV

400
200
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 5.8: Histograms of generated distributions for beryllium at two

different photon energies. The total number of entries at each energy is 106 .

x 100
5000
4500 1TeV 1 GeV
10 GeV
4000
100 GeV
3500 1000 GeV
3000 100 GeV
2500
2000
1500
10 GeV 1 GeV
1000
500
0
0 .01 .02 .03 .04 .05 .06 .07 .08 .09 0.1
+

Figure 5.9: Histograms of generated + distributions at different photon

energies.

46
according to
1   1  
+ = u + cos , = u cos sin .and =
+ 2 2 u
(5.70)
The muon vectors can now be constructed from Eq. (5.53), where 0 is chosen
randomly between 0 and 2. Fig. 5.9 shows distributions of + at different
photon energies (in beryllium). The spectra peak around 1/ as expected.
The most probable values are + m /E+ = 1/+ . In the small angle
approximation used here, the values of + and can in principle be any
positive value from 0 to . In the simulation, this may lead (with a very
small probability, of the order of m /E ) to unphysical events in which + or
is greater than . To avoid this, a limiting angle cut = is introduced,
and the angular sampling repeated, whenever max(+ , ) > cut .
1.4
simulated
1.2 exact

1.0
Coulomb centre
0.8

0.6

0.4

0.2

0.0
0 0.2 0.4 0.6 0.8 1
1 / ( 1 + +2 +2 )

Figure 5.10: Angular distribution of positive (or negative) muons. The solid
curve represents the results of the exact calculations. The histogram is the
simulated distribution. The angular distribution for pairs created in the field
of the Coulomb centre (point-like target) is shown by the dashed curve for
comparison.

Figs. 5.10,5.11 and 5.12 show distributions of the simulated angular char-
acteristics of muon pairs in comparison with results of exact calculations.
The latter were obtained by means of numerical integration of the squared
matrix elements with respective nuclear and atomic form factors. All these
calculations were made for iron, with E = 10 GeV and x+ = 0.3. As seen
from Fig. 5.10, wide angle pairs (at low values of the argument in the fig-
ure) are suppressed in comparison with the Coulomb center approximation.
This is due to the influence of the finite nuclear size which is comparable
to the inverse mass of the muon. Typical angles of particle emission are of

47
 7

0
10 -1 1
+ +

Figure 5.11: Angular distribution in logarithmic scale. The curve corre-

sponds to the exact calculations and the histogram is the simulated distri-
bution.

0 -3 -2 -1
10 10 10 1
| + + - - - |

Figure 5.12: Distribution of the difference of transverse momenta of positive

and negative muons (with logarithmic x-scale).

48
the order of 1/ = m /E (Fig. 5.11). Fig. 5.12 illustrates the influence of
the momentum transferred to the target on the angular characteristics of the
produced pair. In the frame of the often used model which neglects target
recoil, the pair particles would be symmetric in transverse momenta, and
coplanar with the initial photon.

5.5.5 Status of this document

28.05.02 created by H. Burkhardt.
01.12.02 re-worded by D.H. Wright

Bibliography
[1] H. Burkhardt, S. Kelner, and R. Kokoulin, Monte Carlo Generator for
Muon Pair Production. CERN-SL-2002-016 (AP) and CLIC Note 511,
May 2002.

[2] S.R. Kelner, R.P. Kokoulin, and A.A. Petrukhin, About cross section
for high energy muon bremsstrahlung. Moscow Phys. Eng. Inst. 024-95,
1995. 31pp.

49
Chapter 6

Elastic scattering

50
6.1 Multiple Scattering
Elastic scattering of electrons and other charged particles is an important
component of any transport code. Elastic cross section is huge when particle
energy decreases, so multiple scattering (MSC) approach should be intro-
duced in order to have acceptable CPU performance of the simulation. A
universal interface G4VMultipleScattering is used by all Geant4 MSC pro-
cesses [1]:

G4eMultipleScattering;

G4hMultipleScattering;

G4MuMultipleScattering.

For concrete simulation the G4VMscModel interface is used, which is an

extension of the base G4VEmModel interface. The following models are
available:

G4UrbanMscModel - since Geant4 10.0 only one Urban model is avail-

able and it is applicable to all types of particles;

G4GoudsmitSaundersonModel - for electrons and positrons [2];

G4LowEWentzelVIModel - for all particles with low-energy limit 10 eV;

G4WentzelVIModel - for muons and hadrons, for muons should be in-

cluded in Physics List together with G4CoulombScattering process, for
hadrons large angle scattering is simulated by hadron elastic process.

The discussion on Geant4 MSC models is available in Ref.[3]. Below we will

describe models developed by L. Urban [4], because these models are used
in many Geant4 applications and have general components reused by other
models.

6.1.1 Introduction
MSC simulation algorithms can be classified as either detailed or condensed.
In the detailed algorithms, all the collisions/interactions experienced by the
particle are simulated. This simulation can be considered as exact, it gives
the same results as the solution of the transport equation. However, it can
be used only if the number of collisions is not too large, a condition fulfilled
only for special geometries (such as thin foils, or low density gas). In solid

51
or liquid media the average number of collisions is very large and the de-
tailed simulation becomes very inefficient. High energy simulation codes use
condensed simulation algorithms, in which the global effects of the collisions
are simulated at the end of a track segment. The global effects generally
computed in these codes are the net energy loss, displacement, and change
of direction of the charged particle. The last two quantities are computed
from MSC theories used in the codes and the accuracy of the condensed
simulations is limited by accuracy of MSC approximation.
Most particle physics simulation codes use the multiple scattering theo-
ries of Moli`ere [5], Goudsmit and Saunderson [6] and Lewis [7]. The theories
of Moli`ere and Goudsmit-Saunderson give only the angular distribution after
a step, while the Lewis theory computes the moments of the spatial distribu-
tion as well. None of these MSC theories gives the probability distribution
of the spatial displacement. Each of the MSC simulation codes incorporates
its own algorithm to determine the angular deflection, true path length cor-
rection, and spatial displacement of the charged particle after a given step.
These algorithms are not exact, of course, and are responsible for most of
the uncertainties of the transport codes. Also due to inaccuracy of MSC the
simulation results can depend on the value of the step length and generally
user has to select the value of the step length carefully.
A new class of MSC simulation, the mixed simulation algorithms (see
e.g.[8]), appeared in the literature recently. The mixed algorithm simulates
the hard collisions one by one and uses a MSC theory to treat the effects of
the soft collisions at the end of a given step. Such algorithms can prevent
the number of steps from becoming too large and also reduce the dependence
on the step length. Geant4 original implementation of a similar approach is
realized in G4WentzelVIModel [3].
The Urban MSC models used in Geant4 belongs to the class of condensed
simulations. Urban uses model functions to determine the angular and spatial
distributions after a step. The functions have been chosen in such a way as
to give the same moments of the (angular and spatial) distributions as are
given by the Lewis theory [7].

6.1.2 Definition of Terms

In simulation, a particle is transported by steps through the detector ge-
ometry. The shortest distance between the endpoints of a step is called
the geometrical path length, z. In the absence of a magnetic field, this is a
straight line. For non-zero fields, z is the length along a curved trajectory.
Constraints on z are imposed when particle tracks cross volume boundaries.
The path length of an actual particle, however, is usually longer than the ge-

52
ometrical path length, due to multiple scattering. This distance is called the
true path length, t. Constraints on t are imposed by the physical processes
acting on the particle.
The properties of the MSC process are determined by the transport mean
free paths, k , which are functions of the energy in a given material. The
k-th transport mean free path is defined as
Z 1
1 d()
= 2na [1 Pk (cos)] d(cos) (6.1)
k 1 d
where d()/d is the differential cross section of the scattering, Pk (cos)
is the k-th Legendre polynomial, and na is the number of atoms per volume.
Most of the mean properties of MSC computed in the simulation codes
depend only on the first and second transport mean free paths. The mean
value of the geometrical path length (first moment) corresponding to a given
true path length t is given by
  
t
hzi = 1 1 exp (6.2)
1

Eq. 6.2 is an exact result for the mean value of z if the differential cross
section has axial symmetry and the energy loss can be neglected. The trans-
formation between true and geometrical path lengths is called the path length
correction. This formula and other expressions for the first moments of the
spatial distribution were taken from either [8] or [9], but were originally cal-
culated by Goudsmit and Saunderson [6] and Lewis [7].
At the end of the true step length, t, the scattering angle is . The mean
value of cos is  
t
hcosi = exp (6.3)
1
The variance of cos can be written as
1 + 2e2
2 = hcos2 i hcosi2 = e2 (6.4)
3
where = t/1 and = 1 /2 . The mean lateral displacement is given
by a more complicated formula [8], but this quantity can also be calculated
relatively easily and accurately. The square of the mean lateral displacement
is
421
 
2 2 +1 1
hx + y i = + e e (6.5)
3 1 ( 1)

53
Here it is assumed that the initial particle direction is parallel to the the z
axis. The lateral correlation is determined by the equation
 
21 1
hxvx + yvy i = 1 e + e (6.6)
3 1 1

where vx and vy are the x and y components of the direction unit vector. This
equation gives the correlation strength between the final lateral position and
final direction.
The transport mean free path values have been calculated in Refs.[10],[11]
for electrons and positrons in the kinetic energy range 100 eV - 20 MeV in
15 materials. The Urban MSC model in Geant4 uses these values for kinetic
energies below 10 MeV. For high energy particles (above 10 MeV) the trans-
port mean free path values have been taken from a paper of R. Mayol and
F. Salvat [12]. When necessary, the model linearly interpolates or extrap-
olates the transport cross section, 1 = 1/1 , in atomic number Z and in
the square of the particle velocity, 2 . The ratio is a very slowly varying
function of the energy: > 2 for T > a few keV, and 3 for very high
energies (see [9]). Hence, a constant value of 2.5 is used in the model.
Nuclear size effects are negligible for low energy particles and they are
accounted for in the Born approximation in [12], so there is no need for extra
corrections of this kind in the Urban model.

6.1.3 Path Length Correction

As mentioned above, the path length correction refers to the transformation
t g and its inverse. The t g transformation is given by Eq. 6.2 if the
step is small and the energy loss can be neglected. If the step is not small
the energy dependence makes the transformation more complicated. For this
case Eqs. 6.3,6.2 should be modified as
 Z t 
du
hcosi = exp (6.7)
0 1 (u)

Z t
hzi = hcosiu du (6.8)
0
where is the scattering angle, t and z are the true and geometrical path
lengths, and 1 is the transport mean free path.
In order to compute Eqs. 6.7,6.8 the t dependence of the transport mean
free path must be known. 1 depends on the kinetic energy of the particle

54
which decreases along the step. All computations in the model use a linear
approximation for this t dependence:

1 (t) = 10 (1 t) (6.9)

Here 10 denotes the value of 1 at the start of the step, and is a constant.
It is worth noting that Eq. 6.9 is not a crude approximation. It is rather
good at low (< 1 MeV) energy. At higher energies the step is generally much
smaller than the range of the particle, so the change in energy is small and
so is the change in 1 . Using Eqs. 6.7 - 6.9 the explicit formula for hcosi
and hzi are:
1
hcosi = (1 t) 10 (6.10)

1 h
1+ 1
i
hzi = 1 (1 t) 10 (6.11)
(1 + 110 )
The value of the constant can be expressed using 10 and 11 where 11 is
the value of the transport mean free path at the end of the step
10 11
= (6.12)
t10
At low energies ( Tkin < M , M - particle mass) has a simpler form:
1
= (6.13)
r0
where r0 denotes the range of the particle at the start of the step. It can
easily be seen that for a small step (i.e. for a step with small relative energy
loss) the formula of hzi is
  
t
hzi = 10 1 exp (6.14)
10
Eq. 6.11 or 6.14 gives the mean value of the geometrical step length for a
given true step length. The actual geometrical path length is sampled in
the model according to the simple probability density function defined for
v = z/t [0, 1] :
f (v) = (k + 1)(k + 2)v k (1 v) (6.15)
The value of the exponent k is computed from the requirement that f (v)
must give the same mean value for z = vt as Eq. 6.11 or 6.14. Hence
3hzi t
k= (6.16)
t hzi

55
The value of z = vt is sampled using f (v) if k > 0, otherwise z = hzi is
used. The g t transformation is performed using the mean values. The
transformation can be written as
 
z
t(z) = hti = 1 log 1 (6.17)
1
if the geometrical step is small and
1h 1
i
t(z) = 1 (1 wz) w (6.18)

where
1
w =1+
10
if the step is not small, i.e. the energy loss should be taken into account.

6.1.4 Angular Distribution

The quantity u = cos is sampled according to a model function g(u). The
shape of this function has been chosen such that Eqs. 6.3 and 6.4 are satisfied.
The functional form of g is
g(u) = q[pg1 (u) + (1 p)g2 (u)] + (1 q)g3 (u) (6.19)
where 0 p, q 1, and the gi are simple functions of u = cos, normalized
over the range u [1, 1]. The functions gi have been chosen as
g1 (u) = C1 ea(1u) 1 u0 u 1 (6.20)
1
g2 (u) = C2 1 u u0 1 (6.21)
(b u)d
g3 (u) = C3 1u1 (6.22)
where a > 0, b > 0, d > 0 and u0 are model parameters, and the Ci are
normalization constants. It is worth noting that for small scattering angles,
, g1 (u) is nearly Gaussian (exp(2 /202 )) if 02 1/a, while g2 (u) has a
Rutherford-like tail for large , if b 1 and d is not far from 2 .

6.1.5 Determination of the Model Parameters

The parameters a, b, d, u0 and p, q are not independent. The requirement
that the angular distribution function g(u) and its first derivative be contin-
uous at u = u0 imposes two constraints on the parameters:
p g1 (u0 ) = (1 p) g2 (u0 ) (6.23)

56
 d
p a g1 (u0 ) = (1 p) g2 (u0 ) (6.24)
b u0
A third constraint comes from Eq. 6.7 : g(u) must give the same mean value
for u as the theory. It follows from Eqs. 6.10 and 6.19 that
1
q{phui1 + (1 p)hui2 } = [1 t] 10
(6.25)

where huii denotes the mean value of u computed from the distribution gi (u).
The parameter a was chosen according to a modified Highland-Lynch-Dahl
formula for the width of the angular distribution [13], [14].
0.5
a= (6.26)
1 cos(0 )

where 0 is r   
13.6M eV t t
0 = zch 1 + hc ln (6.27)
cp X0 X0
rms
when the original Highland-Lynch-Dahl formula is used. Here 0 = plane
is the width of the approximate Gaussian projected angle distribution, p,
c and zch are the momentum, velocity and charge number of the incident
particle, and t/X0 is the true path length in radiation length unit. The
correction term hc = 0.038 in the formula. This value of 0 is from a fit to
the Moli`ere distribution for singly charged particles with = 1 for all Z,
and is accurate to 11 % or better for 103 t/X0 100 (see e.g. Rev. of
Particle Properties, section 23.3).
The model uses a slightly modified Highland-Lynch-Dahl formula to com-
pute 0 . For electrons/positrons the modified 0 formula is

13.6M eV
0 = zch yc (6.28)
cp
where
 
t
y = ln (6.29)
X0
The correction term c and coeffitients ci are
1 + c2 y),
c = c0 (c (6.30)

c0 = 0.990395 0.168386Z 1/6 + 0.093286Z 1/3 , (6.31)

0.08778
c1 = 1 , (6.32)
Z

57
 c2 = 0.04078 + 0.00017315Z. (6.33)
This formula gives a much smaller step dependence in the angular dis-
tribution than the Highland form. The value of the parameter u0 has been
chosen as

u0 = 1 (6.34)
a
where
= d1 + d2 v + d3 v 2 + d4 v 3 (6.35)
with  
t
v = ln (6.36)
1
The parameters di -s have the form
1 2
di = di0 + di1 Z 3 + di2 Z 3 (6.37)

The numerical values of the dij constants can be found in the code.
The tail parameter d is the same as the parameter .
This (empirical) expression is obtained comparing the simulation results
to the data of the MuScat experiment [16]. The remaining three parame-
ters can be computed from Eqs. 6.23 - 6.25. The numerical value of the
parameters can be found in the code.
In the case of heavy charged particles (, , p, etc.) the mean transport
free path is calculated from the electron or positron 1 values with a scaling
applied. This is possible because the transport mean free path 1 depends
only on the variable P c, where P is the momentum, and c is the velocity
of the particle.
In its present form the model samples the path length correction and an-
gular distribution from model functions, while for the lateral displacement
and the lateral correlation only the mean values are used and all the other
correlations are neglected. However, the model is general enough to incorpo-
rate other random quantities and correlations in the future.

6.1.6 Step Limitation Algorithm

In Geant4 the boundary crossing is treated by the transportation process.
The transportation ensures that the particle does not penetrate in a new
volume without stopping at the boundary, it restricts the step size when the
particle leaves a volume. However, this step restriction can be rather weak
in big volumes and this fact can result a not very good angular distribution
after the volume. At the same time, there is no similar step limitation when

58
a particle enters a volume and this fact does not allow a good backscattering
simulation for low energy particles. Low energy particles penetrate too deeply
into the volume in the first step and then - because of energy loss - they are
not able to reach again the boundary in backward direction.
MSC step limitation algorithm has been developed [4] in order to achieve
optimal balance between simulation precision and CPU performance of sim-
ulation for different applications. At the start of a track or after entering in
a new volume, the algorithm restricts the step size to a value

fr max{r, 1 } (6.38)

where r is the range of the particle, fr is a parameter [0, 1], taking the max
of r and 1 is an empirical [Link] value of fr is constant for low energy
particles while for particles with 1 > lim an effective value is used given by
the scaling equation
 
1
fref f = fr 1 sc + sc (6.39)
lim
( The numerical values sc = 0.25 and lim = 1 mm are used in the equation.)
In order not to use very small - unphysical - step sizes a lower limit is given
for the step size as
 
1
tlimitmin = max , elastic (6.40)
nstepmax
with nstepmax = 25 and elastic is the elastic mean free path of the particle
(see later).

It can be easily seen that this kind of step limitation poses a real constraint
only for low energy particles. In order to prevent a particle from crossing a
volume in just one step, an additional limitation is imposed: after entering
a volume the step size cannot be bigger than
dgeom
(6.41)
fg
where dgeom is the distance to the next boundary (in the direction of the
particle) and fg is a constant parameter. A similar restriction at the start of
a track is
2dgeom
(6.42)
fg
At this point the program also checks whether the particle has entered a
new volume. If it has, the particle steps cannot be bigger than tlim =

59
fr max(r, ). This step limitation is governed by the physics, because tlim
depends on the particle energy and the material.
The choice of the parameters fr and fg is also related to performance.
By default fr = 0.02 and fg = 2.5 are used, but these may be set to any
other value in a simple way. One can get an approximate simulation of
the backscattering with the default value, while if a better backscattering
simulation is needed it is possible to get it using a smaller value for fr .
However, this model is very simple and it can only approximately reproduce
the backscattering data.

6.1.7 Boundary Crossing Algorithm

A special stepping algorithm has been implemented in order to improve the
simulation around interfaces. This algorithm does not allow big last steps
in a volume and big first steps in the next volume. The step length of these
steps around a boundary crossing can not be bigger than the mean free path
of the elastic scattering of the particle in the given volume (material). After
these small steps the particle scattered according to a single scattering law
(i.e. there is no multiple scattering very close to the boundary or at the
boundary).
The key parameter of the algorithm is the variable called skin. The
algorithm is not active for skin 0, while for skin > 0 it is active in
layers of thickness skin elastic before boundary crossing and of thickness
(skin1)elastic after boundary crossing (for skin = 1 there is only one small
step just before the boundary). In this active area the particle performs steps
of length elastic (or smaller if the particle reaches the boundary traversing a
smaller distance than this value).
The scattering at the end of a small step is single or plural and for these
small steps there are no path length correction and lateral displacement com-
putation. In other words the program works in this thin layer in microscopic
mode. The elastic mean free path can be estimated as
elastic = 1 rat (Tkin ) (6.43)
where rat(Tkin ) a simple empirical function computed from the elastic and
first transport cross section values of Mayol and Salvat [12]

0.001(M eV )2
rat (Tkin ) = (6.44)
Tkin (Tkin + 10M eV )
Tkin is the kinetic energy of the particle.
At the end of a small step the number of scatterings is sampled according
to the Poissons distribution with a mean value t/elastic and in the case of

60
plural scattering the final scattering angle is computed by summing the con-
tributions of the individual scatterings. The single scattering is determined
by the distribution
1
g(u) = C (6.45)
(2a + 1 u)2
where u = cos() , a is the screening parameter, C is a normalization con-
stant. The form of the screening parameter is the same as in the single
scattering (see there).

6.1.8 Implementation Details

The step length of a particles is determined by the physics processes or the
geometry of the detectors. The tracking/stepping algorithm checks all the
step lengths demanded by the (continuous or discrete) physics processes and
determines the minimum of these step lengths (see 3.2). The MSC model
should be called to compute step limit after all processes except the trans-
portation process. The following sequence of computations are performed to
make the step:
the minimum of all processes true step length limit t including one of
the MSC process is selected;
The conversion t g (geometrical step limit) is performed;
the minimum of obtained value g and the transportation step limit is
selected;
The final conversion g t is performed.
The reason for this ordering is that the physics processes feel the true path
length t traveled by the particle, while the transportation process (geometry)
uses the z step length.
A new optional mechanis was recently introduced allowing sample dis-
placemnt in vicinity of geometry boundary. If it is enabled and the trans-
portation limit the step due to geometry boundary then after initial sampling
of the displacenet an additional push of track is applied forcing end point be
at the boundary. Corresponding correction to the true step length is applied
according to the value of the push.
After the actual step of the particle is done, the MSC model is responsible
for sampling of scattering angle and relocation of the end-point of the step.
The scattering angle of the particle after the step of length t is sampled
according to the model function given in Eq. 6.19 . The azimuthal angle
is generated uniformly in the range [0, 2].

61
 After the simulation of the scattering angle, the lateral displacement is
computed using Eq. 6.5. Then the correlation given by Eq. 6.6 is used to
determine the direction of the lateral displacement. Before moving the
particle according to the displacement a check is performed to ensure that
the relocation of the particle with the lateral displacement does not take the
particle beyond the volume boundary.
Default MSC parameter values optimized per particle type are shown in
Table 6.1. Note, that there is three types of step limitation by multiple
scattering process:

Minimal - only fr parameter and range are used;

UseSafety - fr parameter, range and geometrical safety are used;

UseSafetyPlus - fr parameter, range and geometrical safety are used;

UseDistanceToBoundary - uses particle range, geometrical safety and

linear distance to geometrical boundary.

particle e+ , e muons, hadrons ions

StepLimitType fUseSafety fMinimal fMinimal
skin 0 0 0
fr 0.04 0.2 0.2
fg 2.5 0.1 0.1
LateralDisplacement true true false

Table 6.1: The default values of parameters for different particle type.

The parameters of the model can be changed via public functions of the base
class G4VMultipleSacttering. They can be changed for all multiple scattering
processes simultaneously via G4EmParameters class, G4EmProcessOptions
class, or via Geant4 UI commands. The following commands are available:

/process/msc/StepLimit UseDistanceToBoundary
/process/msc/LateralDisplacement false
/process/msc/MuHadLateralDisplacement false
/process/msc/DisplacementBeyondSafety true
/process/msc/RangeFactor 0.02
/process/msc/GeomFactor 2.5
/process/msc/Skin 2

62
6.1.9 Status of this document
09.10.98 created by L. Urban.
15.11.01 major revision by L. Urban.
18.04.02 updated by L. Urban.
25.04.02 re-worded by D.H. Wright
07.06.02 major revision by L. Urban.
18.11.02 updated by L. Urban, now it describes the new angle distribution.
05.12.02 grammar check and parts re-written by D.H. Wright
13.11.03 revision by L. Urban.
01.12.03 revision by V. Ivanchenko.
17.05.04 revision by L. Urban.
01.12.04 updated by L. Urban.
18.03.05 sampling z + mistyping corrections (M. Maire)
22.06.05 grammar, spelling check by D.H. Wright
12.12.05 revised by L. Urban, according to Geant4 8.0
14.12.05 updated implementation Details (M. Maire)
08.06.06 revised by L. Urban, according to Geant4 8.1
25.11.06 revised by L. Urban, according to Geant4 8.2
29.03.07 revised by L. Urban, for Geant4 8.3
13.06.07 modified introduction (M. Maire)
17.06.07 explain effective FR (L. Urban)
25.06.07 update description of options by V. Ivanchenko
05.12.07 revised by L. Urban, for Geant4 9.1
08.12.08 revised by L. Urban, for Geant4 9.2
11.12.08 minor revision by V. Ivanchenko
11.12.09 minor revision by V. Ivanchenko, for Geant4 9.3
09.12.09 revision by V. Ivanchenko, for Geant4 9.4
25.11.11 minor revision by V. Ivanchenko, for Geant4 9.5
03.12.13 minor revision by L. Urban, for Geant4 10.0

Bibliography
[1] J. Apostolakis et al., Geometry and physics of the Geant4 toolkit for high
and medium energy applications. Rad. Phys. Chem. 78 (2009) 859.

[2] O. Kadri, V. Ivanchenko, F. Gharbi, A. Trabelsi, Incorporation of the

Goudsmit-Saunderson electron transport theory in the Geant4 Monte
Carlo code, Nucl. Instrum. and Meth. B 267 (2009) 3624.

63
 [3] V.N. Ivanchenko et al., Geant4 models for simulation of multiple scat-
tering, J. Phys.: Conf. Ser. 219 (2010) 032045.

[4] L. Urban, A multiple scattering model, CERN-OPEN-2006-077, Dec

2006. 18 pp.

[5] G.Z. Moli`ere Z. Naturforsch. 3a (1948) 78.

[6] S. Goudsmit and J.L. Saunderson. Phys. Rev. 57 (1940) 24.

[7] H.W. Lewis. Phys. Rev. 78 (1950) 526.

[8] J.M. Fernandez-Varea et al. NIM B73 (1993) 447.

[9] I. Kawrakow and A.F. Bielajew NIM B 142 (1998) 253.

[10] D. Liljequist and M. Ismail. [Link]. 62 (1987) 342.

[11] D. Liljequist et al. [Link]. 68 (1990) 3061.

[12] R. Mayol and F. Salvat [Link] and [Link] Tables 65 (1997) 55..

[13] V.L. Highland NIM 129 (1975) 497.

[14] G.R. Lynch and O.I. Dahl NIM B58 (1991) 6.

[15] G. Shen et al. Phys. Rev. D 20 (1979) 1584.

[16] D. Attwood et al. NIM B 251 (2006) 41.

64
6.2 Discrete Processes for Charged Particles
Some processes for charged particles following the same interface G4V EmP rocess
as gamma processes described in section 5.1:
G4CoulombScattering;
G4eplusAnnihilation (with additional AtRest methods);
G4eplusPolarizedAnnihilation (with additional AtRest methods);
G4eeToHadrons;
G4NuclearStopping;
G4MicroElecElastic;
G4MicroElecInelastic.
Corresponding model classes follow the G4V EmM odel interface:
G4DummyModel (zero cross section, no secondaries);
G4eCoulombScatteringModel;
G4eSingleCoulombScatteringModel;
G4IonCoulombScatteringModel;
G4eeToHadronsModel;
G4PenelopeAnnihilationModel;
G4PolarizedAnnihilationModel;
G4ICRU49NuclearStoppingModel;
G4MicroElecElasticModel;
G4MicroElecInelasticModel.
Some processes from do not follow described EM interfaces but provide direct
implementations of the basic G4V DiscreteP rocess process:
G4AnnihiToMuPair;
G4ScreenedNuclearRecoil;
G4Cerenkov;
G4Scintillation;
G4SynchrotronRadiation;

65
6.2.1 Status of This Document
10.12.10 created by V. Ivanchenko
29.11.13 updated by V. Ivanchenko

66
6.3 Single Scattering
Single elastic scattering process is an alternative to the multiple scattering
process. The advantage of the single scattering process is in possibility of
usage of theory based cross sections, in contrary to the Geant4 multiple scat-
tering model [1], which uses a number of phenomenological approximations
on top of Lewis theory. The process G4CoulombScattering was created for
simulation of single scattering of muons, it also applicable with some physical
limitations to electrons, muons and ions. Because each of elastic collisions are
simulated the number of steps of charged particles significantly increasing in
comparison with the multiple scattering approach, correspondingly its CPU
performance is pure. However, in low-density media (vacuum, low-density
gas) multiple scattering may provide wrong results and single scattering pro-
cesses is more adequate.

6.3.1 Coulomb Scattering

The single scattering model of Wentzel [2] is used in many of multiple scat-
tering models including Penelope code [4]. The Wentzel for describing elastic
scattering of particles with charge ze (z = 1 for electron) by atomic nucleus
with atomic number Z based on simplified scattering potential
zZe2
V (r) = exp(r/R), (6.46)
r
where the exponential factor tries to reproduce the effect of screening. The
parameter R is a screening radius [3]
R = 0.885Z 1/3 rB , (6.47)
where rB is the Bohr radius. In the first Born approximation the elastic
scattering cross section ( W ) can be obtained as
d (W ) () (ze2 )2 Z(Z + 1)
= , (6.48)
d (pc)2 (2A + 1 cos)2
where p is the momentum and is the velocity of the projectile particle. The
screening parameter A according to Moliere and Bethe [3]
 2

A= (1.13 + 3.76(Z/)2 ), (6.49)
2pR
where is a fine structure constant and the factor in brackets is used to take
into account second order corrections to the first Born approximation.

67
 The total elastic cross section can be expressed via Wentzel cross section
(6.48) !
d() d (W ) () Z 1
= (qR ) 2 +1 , (6.50)
d d (1 + N )2 Z +1
12
where q is momentum transfer to the nucleus, RN is nuclear radius. This
term takes into account nuclear size effect [5], the second term takes into
account scattering off electrons. The results of simulation with the single
scattering model (Fig.6.1) are competitive with the results of the multiple
scattering.

Figure 6.1: Scattering of muons off 1.5 mm aluminum foil: data [6] - black
squares; simulation - colored markers corresponding different options of mul-
tiple scattering and single scattering model; in the bottom plot - relative
difference between the simulation and the data in percents; hashed area
demonstrates one standard deviation of the data.

6.3.2 Implementation Details

The total cross section of the process is obtained as a result of integration
of the differential cross section (6.50). The first term of this cross section
is integrated in the interval (0, ). The second term in the smaller interval
(0, m ), where m is the maximum scattering angle off electrons, which is
determined using the cut value for the delta electron production. Before

68
sampling of angular distribution the random choice is performed between
scattering off the nucleus and off electrons.

6.3.3 Status of This Document

06.12.07 created by V. Ivanchenko
08.12.10 added chapter on discrete processes by [Link]

Bibliography
[1] L. Urban, A multiple scattering model, CERN-OPEN-2006-077, Dec
2006. 18 pp.

[2] G. Wentzel, Z. Phys. 40 (1927) 590.

[3] H.A. Bethe, Phys. Rev. 89 (1953) 1256.

[4] J.M. Fernandez-Varea et al. NIM B 73 (1993) 447.

[5] A.V. Butkevich et al., NIM A 488 (2002) 282.

[6] D. Attwood et al. NIM B 251 (2006) 41.

69
6.4 Ion Scattering
The necessity of accurately computing the characteristics of interatomic scat-
tering arises in many disciplines in which energetic ions pass through mate-
rials. Traditionally, solutions to this problem not involving hadronic inter-
actions have been dominated by the multiple scattering, which is reasonably
successful, but not very flexible. In particular, it is relatively difficult to in-
troduce into such a system a particular screening function which has been
measured for a specific atomic pair, rather than the universal functions which
are applied. In many problems of current interest, such as the behavior of
semiconductor device physics in a space environment, nuclear reactions, par-
ticle showers, and other effects are critically important in modeling the full
details of ion transport. The process G4ScreenedNuclearRecoil provides sim-
ulation of ion elastic scattering [1]. This process is available with extended
electromagnetic example TestEm7.

6.4.1 Method
The method used in this computation is a variant of a subset of the method
described in Ref.[2]. A very short recap of the basic material is included here.
The scattering of two atoms from each other is assumed to be a completely
classical process, subject to an interatomic potential described by a potential
function
Z1 Z2 e2  r 
V (r) = (6.51)
r a
where Z1 and Z2 are the nuclear proton numbers, e2 is the electromagnetic
coupling constant (qe2 /40 in SI units), r is the inter-nuclear separation,
is the screening function describing the effect of electronic screening of the
bare nuclear charges, and a is a characteristic length scale for this screening.
In most cases, is a universal function used for all ion pairs, and the value of
a is an appropriately adjusted length to give reasonably accurate scattering
behavior. In the method described here, there is no particular need for
a universal function , since the method is capable of directly solving the
problem for most physically plausible screening functions. It is still useful
to define a typical screening length a in the calculation described below, to
keep the equations in a form directly comparable with our previous work even
though, in the end, the actual value is irrelevant as long as the final function
(r) is correct. From this potential V (r) one can then compute the classical
scattering angle from the reduced center-of-mass energy Ec a/Z1 Z2 e2
(where Ec is the kinetic energy in the center-of-mass frame) and reduced

70
impact parameter b/a
Z
c = 2 f (z) dz/z 2 (6.52)
x0

where 1/2
(z) 2

f (z) = 1 2 (6.53)
z z
and x0 is the reduced classical turning radius for the given and .
The problem, then, is reduced to the efficient computation of this scat-
tering integral. In our previous work, a great deal of analytical effort was
included to proceed from the scattering integral to a full differential cross
section calculation, but for application in a Monte-Carlo code, the scattering
integral c (Z1 , Z2 , Ec , b) and an estimated total cross section 0 (Z1 , Z2 , Ec )
are all that is needed. Thus, we can skip algorithmically forward in the orig-
inal paper to equations 15-18 and the surrounding discussion to compute the
reduced distance of closest approach x0 . This computation follows that in
the previous work exactly, and will not be reintroduced here.
For the sake of ultimate accuracy in this algorithm, and due to the rela-
tively low computational cost of so doing, we compute the actual scattering
integral (as described in equations 19-21 of [2]) using a Lobatto quadrature
of order 6, instead of the 4th order method previously described. This re-
sults in the integration accuracy exceeding that of any available interatomic
potentials in the range of energies above those at which molecular structure
effects dominate, and should allow for future improvements in that area. The
integral then becomes (following the notation of the previous paper)
4  
1 + 0 X x0
+ wi f (6.54)
30 i=1
q i

where 1/2
2 (x0 )

1
0 = + (6.55)
2 2 x20 2

wi [0.03472124, 0.1476903, 0.23485003, 0.1860249]

qi [0.9830235, 0.8465224, 0.5323531, 0.18347974]
Then

c = (6.56)
x0
The other quantity required to implement a scattering process is the total
scattering cross section 0 for a given incident ion and a material through

71
which the ion is propagating. This value requires special consideration for a
process such as screened scattering. In the limiting case that the screening
function is unity, which corresponds to Rutherford scattering, the total cross
section is infinite. For various screening functions, the total cross section
may or may not be finite. However, one must ask what the intent of defining
a total cross section is, and determine from that how to define it.
In Geant4, the total cross section is used to determine a mean-free-path
l which is used in turn to generate random transport distances between
discrete scattering events for a particle. In reality, where an ion is propagating
through, for example, a solid material, scattering is not a discrete process
but is continuous. However, it is a useful, and highly accurate, simplification
to reduce such scattering to a series of discrete events, by defining some
minimum energy transfer of interest, and setting the mean free path to be
the path over which statistically one such minimal transfer has occurred. This
approach is identical to the approach developed for the original TRIM code
[3]. As long as the minimal interesting energy transfer is set small enough
that the cumulative effect of all transfers smaller than that is negligible,
the approximation is valid. As long as the impact parameter selection is
adjusted to be consistent with the selected value of l , the physical result
isnt particularly sensitive to the value chosen.
Noting, then, that the actual physical result isnt very sensitive to the
selection of l , one can be relatively free about defining the cross section 0
from which l is computed. The choice used for this implementation is fairly
simple. Define a physical cutoff energy Emin which is the smallest energy
transfer to be included in the calculation. Then, for a given incident particle
with atomic number Z1 , mass m1 , and lab energy Einc , and a target atom
with atomic number Z2 and mass m2 , compute the scattering angle c which
will transfer this much energy to the target from the solution of
4 m1 m2 c
Emin = Einc 2
sin2 (6.57)
(m1 + m2 ) 2
. Then, noting that from eq. 6.54 is a number very close to unity, one
can solve for an approximate impact parameter b with a single root-finding
operation to find the classical turning point. Then, define the total cross
section to be 0 = b2 , the area of the disk inside of which the passage of an
ion will cause at least the minimum interesting energy transfer. Because this
process is relatively expensive, and the result is needed extremely frequently,
the values of 0 (Einc ) are precomputed for each pairing of incident ion and
target atom, and the results cached in a cubic-spline interpolation table.
However, since the actual result isnt very critical, the cached results can be
stored in a very coarsely sampled table without degrading the calculation at

72
all, as long as the values of the l used in the impact parameter selection are
rigorously consistent with this table.
The final necessary piece of the scattering integral calculation is the sta-
tistical selection of the impact parameter b to be used in each scattering
event. This selection is done following the original algorithm from TRIM,
where the cumulative probability distribution for impact parameters is

b2
 
P (b) = 1 exp (6.58)
0

where N 0 1/l where N is the total number density of scattering centers

in the target material and l is the mean free path computed in the conven-
tional way. To produce this distribution from a uniform random variate r on
(0,1], the necessary function is
s
log r
b= (6.59)
N l

This choice of sampling function does have the one peculiarity that it can
produce values of the impact parameter which are larger than the impact
parameter which results in the cutoff energy transfer, as discussed above
in the section on the total cross section, with probability 1/e. When this
occurs, the scattering event is not processed further, since the energy transfer
is below threshold. For this reason, impact parameter selection is carried out
very early in the algorithm, so the effort spent on uninteresting events is
minimized.
The above choice of impact

sampling is modified when the mean-free-path
l 2
is very short. If 0 > 2 where l is the approximate lattice constant of
the material, as defined by l = N 1/3 , the sampling is replaced by uniform
sampling on a disk of radius l/2, so that

l
b= r (6.60)
2
This takes into account that impact parameters larger than half the lattice
spacing do not occur, since then one is closer to the adjacent atom. This also
derives from TRIM.
One extra feature is included in our model, to accelerate the produc-
tion of relatively rare events such as high-angle scattering. This feature is a
cross-section scaling algorithm, which allows the user access to an unphys-
ical control of the algorithm which arbitrarily scales the cross-sections for
a selected fraction of interactions. This is implemented as a two-parameter

73
adjustment to the central algorithm. The first parameter is a selection fre-
quency fh which sets what fraction of the interactions will be modified. The
second parameter is the scaling factor for the cross-section. This is imple-
mentedby, for a fraction fh of interactions, scaling the impact parameter by

b = b/ scale. This feature, if used with care so that it does not provide
excess multiple-scattering, can provide between 10 and 100-fold improve-
ments to event rates. If used without checking the validity by comparing to
un-adjusted scattering computations, it can also provide utter nonsense.

6.4.2 Implementation Details

The coefficients for the summation to approximate the integral for in
eq.(6.54) are derived from the values in Abramowitz & Stegun [4], altered to
make the change-of-variable used for this integral. There are two basic steps
to the transformation. First, since the provided abscissas xi and weights wi
are for integration on [-1,1], with only one half of the values provided, and
in this work the integration is being carried out on [0,1], the abscissas are
transformed as:  
1 xi
yi (6.61)
2
Then, the primary change-of-variable is applied resulting in:
yi
qi = cos (6.62)
2
wi yi
wi = sin (6.63)
2 2
except for the first coefficient w1 where the sin() part of the weight is taken
into the limit of 0 as described in eq.(6.55). This value is just w1 = w1 /2.

6.4.3 Status of this document

06.12.07 created by V. Ivanchenko from paper of M.H. Mendenhall and
R.A. Weller
06.12.07 further edited by M. Mendenhall to bring contents of paper up-to-
date with current implementation.

Bibliography
[1] M.H. Mendenhall, R.A. Weller, An algorithm for computing screened
Coulomb scattering in Geant4, Nucl. Instr. Meth. B 227 (2005) 420.

74
[2] M.H. Mendenhall, R.A. Weller, Algorithms for the rapid computation
of classical cross sections for screened coulomb collisions, Nucl. Instr.
Meth. in Physics Res. B58 (1991) 11.

[3] J.P. Biersack, L.G. Haggmark, A Monte Carlo computer program for
the transport of energetic ions in amorphous targets, Nucl. Instr. Meth.
in Physics Res. 174 (1980) 257.

[4] M. Abramowitz, I. Stegun (Eds.), Handbook of Mathematical Functions,

Dover, New York, 1965, pp. 888, 920.

75
6.5 Single Scattering, Screened Coulomb Po-
tential and NIEL
Alternative model of Coulomb scattering of ions have been developed based
on [1] and references therein. The advantage of this model is the wide appli-
cability range in energy from 50 keV to 100 T eV per nucleon.

6.5.1 NucleusNucleus Interactions

As discussed in Ref. [1], at small distances from the nucleus, the potential
energy is a Coulomb potential, while - at distances larger than the Bohr
radius - the nuclear field is screened by the fields of atomic electrons. The
interaction between two nuclei is usually described in terms of an interatomic
Coulomb potential (e.g., see Section [Link] of Ref. [2] and Section 4.1 of
Ref. [3]), which is a function of the radial distance r between the two nuclei
zZe2
V (r) = I (rr ), (6.64)
r
where ez (projectile) and eZ (target) are the charges of the bare nuclei and
I is the interatomic screening function and rr is given by
r
rr = , (6.65)
aI
with aI the so-called screening length (also termed screening radius). In the
framework of the ThomasFermi model of the atom (e.g., see Ref. [1] and
references therein) - thus, following the approach of ICRU Report 49 (1993)
-, a commonly used screening length for z = 1 incoming particles is that from
ThomasFermi
CTF a0
aTF = , (6.66)
Z 1/3
and - for incoming particles with z 2 - that introduced by Ziegler, Biersack
and Littmark (1985) (and termed universal screening length):
CTF a0
aU = , (6.67)
z 0.23
+ Z 0.23
where
~2
a0 =
me2
is the Bohr radius, m is the electron rest mass and
 2/3
1 3
CTF = 0.88534
2 4

76
is a constant introduced in the ThomasFermi model.
The simple scattering model due to Wentzel [5] - with a single exponen-
tial screening-function I (rr ) {e.g., see Ref. [1] and references therein} - was
repeatedly employed in treating single and multiple Coulomb-scattering with
screened potentials. The resulting elastic differential cross section differs from
the Rutherford differential cross section by an additional term - the so-called
screening parameter - which prevents the divergence of the cross section when
the angle of scattered particles approaches 0 . The screening parameter As
[e.g., see Equation (21) of Bethe (1953)] - as derived by Moli`ere (1947, 1948)
for the single Coulomb scattering using a ThomasFermi potential - is ex-
pressed as
 2 "  2 #
~ zZ
As = 1.13 + 3.76 (6.68)
2 p aI
where aI is the screening length - from Eqs. (6.66, 6.67) for particles with
z = 1 and z 2, respectively; is the fine-structure constant; p (c) is
the momentum (velocity) of the incoming particle undergoing the scattering
onto a target supposed to be initially at rest; c and ~ are the speed of light
and the reduced Planck constant, respectively. When the (relativistic) mass
- with corresponding rest mass m - of the incoming particle is much lower
than the rest mass (M ) of the target nucleus, the differential cross section -
obtained from the WentzelMoli`ere treatment of the single scattering - is:
2
d WM () zZe2

1
= 2 . (6.69)
d 2 p c As + sin2 (/2)


Equation (6.69) differs from Rutherfords formula - as already mentioned -

for the additional term As to sin2 (/2). As discussed in Ref. [1], for 1
(i.e., at very large p) and with As 1, one finds that the cross section
approaches a constant:
2
2 zZe2 aI

WM
c . (6.70)
~c 1.13 + 3.76 (zZ)2
As discussed in Ref. [1] and references therein, for a scattering under the
action of a central potential (for instance that due to a screened Coulomb
field), when the rest mass of the target particle is no longer much larger than
the relativistic mass of the incoming particle, the expression of the differential
cross section must properly be re-written - in the center of mass system - in
terms of an effective particle with momentum equal to that of the incoming
particle (pin ) and rest mass equal to the relativistic reduced mass
mM
rel = , (6.71)
M1,2

77
where M1,2 is the invariant mass; m and M are the rest masses of the incoming
and target particles, respectively. The effective particle velocity is given by:
v"
u  2 #1
u rel c
r c = ct 1 + .
pin

Thus, one finds (e.g, see Ref. [1]):

2
d WM ( ) zZe2

1
= 2 , (6.72)
d 2 pin r c As + sin2 ( /2)


with  2 "  2 #
~ zZ
As = 1.13 + 3.76 (6.73)
2 pin aI r
and the scattering angle in the center of mass system.
The energy T transferred to the recoil target is related to the scattering
angle as T = Tmax sin2 ( /2) - where Tmax is the maximum energy which
can be transferred in the scattering (e.g., see Section 1.5 of Ref. [2]) -, thus,
assuming an isotropic azimuthal distribution one can re-write Eq. (6.72) in
terms of the kinetic recoil energy T of the target
2
d WM (T ) zZe2

Tmax
= . (6.74)
dT pin r c [Tmax As + T ]2

Furthermore, one can demonstrates that Eq. (6.74) can be re-written as

(e.g, see Ref. [1]);

d WM (T )
E2
2 2 1
= 2 zZe (6.75)
dT p M c [Tmax As + T ]2
2 4

with p and E the momentum and total energy of the incoming particle in the
laboratory. Equation (6.75) expresses - as already mentioned - the differential
cross section as a function of the (kinetic) energy T achieved by the recoil
target.

6.5.2 Nuclear Stopping Power

Using Eq. (6.75) the nuclear stopping power - in MeV cm1 - is obtained as

E2
    
dE 2 2 As As + 1
= 2 nA zZe 1 + ln (6.76)
dx nucl p2 M c4 As + 1 As

78
 4
10

3
10 in silicon

[MeV cm g ]
-1
2
10

2
1 208
10 Pb
115
In
0 56
10 Fe
Stopping power

28
Si
-1
10
12
C

-2 11
10 B

-3 alpha
10
proton

-4
10
-1 0 1 2 3 4 5 6 7 8
10 10 10 10 10 10 10 10 10 10
Kinetic Energy [MeV/nucleon]

Figure 6.2: Nuclear stopping power from Ref. [1] - in MeV cm2 g1 - calcu-
lated using Eq. (6.76) in silicon is shown as a function of the kinetic energy per
nucleon - from 50 keV/nucleon up 100 TeV/nucleon - for protons, -particle
and 11 B-, 12 C-, 28 Si-, 56 Fe-, 115 In-, 208 Pb-nuclei.

with nA the number of nuclei (atoms) per unit of volume and, finally, the
negative sign indicates that the energy is lost by the incoming particle (thus,
achieved by recoil targets). As discussed in Ref. [1], a slight increase of the
nuclear stopping power with energy is expected because of the decrease of
the screening parameter with energy.
For instance, in Fig. 6.2 the nuclear stopping power in silicon - in MeV cm2 g1
- is shown as a function of the kinetic energy per nucleon - from 50 keV/nucleon
up 100 TeV/nucleon - for protons, -particles and 11 B-, 12 C-, 28 Si-, 56 Fe-,
115
In-, 208 Pb-nuclei.
A comparison of the present treatment with that obtained from Ziegler,
Biersack and Littmark (1985) - available in SRIM (2008) [8] - using the so-
called universal screening potential (see also Ref. [9]) is discussed in Ref. [1]:
a good agreement is achieved down to about 150 keV/nucleon. At large en-
ergies, the non-relativistic approach due to Ziegler, Biersack and Littmark
(1985) becomes less appropriate and deviations from stopping powers cal-
culated by means of the universal screening potential are expected and ob-
served.
The non-relativistic approach - based on the universal screening potential
- of Ziegler, Biersack and Littmark (1985) was also used by ICRU (1993) to
calculate nuclear stopping powers due to protons and -particles in materi-
als. ICRU (1993) used as screening lengths those from Eqs. (6.66, 6.67) for
protons and -particles, respectively. As discussed in Ref. [1], the stopping
powers for protons (-particles) from Eq. (6.76) are less than 5% larger

79
 3
10
in silicon

2
10

NIEL [MeV cm g ]
-1
1
10

2
208
0 Pb
10 115
In
56
-1 Fe
10 28
Si

-2 12
10 C

11

-3 B
10 alpha

proton
-4
10

-1 0 1 2 3 4 5 6 7 8
10 10 10 10 10 10 10 10 10 10
Kinetic Energy [MeV/nucleon]

Figure 6.3: Non-ionizing stopping power from Ref. [1] - in MeV cm2 g1 -
calculated using Eq. (6.79) in silicon is shown as a function of the kinetic
energy per nucleon - from 50 keV/nucleon up 100 TeV/nucleon - for protons,
-particles and 11 B-, 12 C-, 28 Si-, 56 Fe-, 115 In-, 208 Pb-nuclei. The threshold
energy for displacement is 21 eV in silicon.

than those reported by ICRU (1993) from 50 keV/nucleon up to 8 MeV

(19 MeV/nucleon). At larger energies the stopping powers from Eq. (6.76)
differ from those from ICRU - as expected - due to the complete relativistic
treatment of the present approach (see Ref. [1]).
The simple screening parameter used so far [Eq. (6.73)] - derived by
Moli`ere (1947) - can be modified by means of a practical correction, i.e.,
 2 "  2 #
~ zZ
As = 1.13 + 3.76 C , (6.77)
2 pin aI r

to achieve a better agreement with low energy calculations of Ziegler, Biersack

and Littmark (1985). For instance - as discussed in Ref. [1] -, for -particles
and heavier ions, with
C = (10zZ)0.12 (6.78)
the stopping powers obtained from Eq. (6.76) - in which As replaces As - differ
from the values of SRIM (2008) by less than 4.7 (3.6) % for -particles (lead
ions) in silicon down to about 50 keV/nucleon. With respect to the tabulated
values of ICRU (1993), the agreement for -particles is usually better than
4% at low energy down to 50 keV/nucleon - a 5% agreement is achieved at
about 50 keV/nucleon in case of a lead medium. At very high energy, the
stopping power is slightly affected when As replaces As (a further disvussion
is found in Ref. [1]).

80
6.5.3 Non-Ionizing Energy Loss due to Coulomb Scat-
tering
A relevant process - which causes permanent damage to the silicon bulk struc-
ture - is the so-called displacement damage (e.g., see Chapter 4 of Ref. [2],
Ref. [10] and references therein). Displacement damage may be inflicted when
a primary knocked-on atom (PKA) is generated. The interstitial atom and
relative vacancy are termed Frenkel-pair (FP). In turn, the displaced atom
may have sufficient energy to migrate inside the lattice and - by further col-
lisions - can displace other atoms as in a collision cascade. This displacement
process modifies the bulk characteristics of the device and causes its degra-
dation. The total number of FPs can be estimated calculating the energy
density deposited from displacement processes. In turn, this energy density
is related to the Non-Ionizing Energy Loss (NIEL), i.e., the energy per unit
path lost by the incident particle due to displacement processes.
In case of Coulomb scattering on nuclei, the non-ionizing energy-loss can
be calculated using the WentzelMoli`ere differential cross section [Eq. (6.75)]
discussed in Sect. 6.5.1, i.e.,
NIEL Tmax
d WM (T )

dE
Z
= nA T L(T ) dT , (6.79)
dx nucl Td dT

where E is the kinetic energy of the incoming particle, T is the kinetic energy
transferred to the target atom, L(T ) is the fraction of T deposited by means
of displacement processes. The expression of L(T ) - the so-called Lindhard
partition function - can be found, for instance, in Equations (4.94, 4.96) of
Section [Link] in Ref. [2] (see also references therein). Tde = T L(T ) is the
so-called damage energy, i.e., the energy deposited by a recoil nucleus with
kinetic energy T via displacement damages inside the medium. The integral in
Eq. (6.79) is computed from the minimum energy Td - the so-called threshold
energy for displacement, i.e., that energy necessary to displace the atom from
its lattice position - up to the maximum energy Tmax that can be transferred
during a single collision process. Td is about 21 eV in silicon. For instance, in
Fig. 6.3 the non-ionizing energy loss - in MeV cm2 g1 - in silicon is shown
as a function of the kinetic energy per nucleon - from 50 keV/nucleon up
100 TeV/nucleon - for protons, -particles and 11 B-, 12 C-, 28 Si-, 56 Fe-, 115 In-,
208
Pb-nuclei.
A further discussion on the agreement with the results obtained by Jun
and collaborators (2003) - using a relativistic treatment of Coulomb scat-
tering of protons with kinetic energies above 50 MeV and up to 1 GeV upon
silicon - can be found in Ref. [1].

81
6.5.4 G4IonCoulombScatteringModel
As discussed sofar, high energetic particles may inflict permanent damage to
the electronic devices employed in a radiation environment. In particular the
nuclear energy loss is important for the formation of defects in semiconductor
devices. Nuclear energy loss is also responsible for the displacement damage
which is the typical cause of degradation for silicon devices. The electromag-
netic model G4IonCoulombScatteringModel was created in order to simulate
the single scattering of protons, alpha particles and all heavier nuclei inci-
dent on all target materials in the energy range from 50100 keV/nucleon to
10 TeV.

6.5.5 The Method

The differential cross section previously described is calculated by means
of the class G4IonCoulombCrossSection where a modified version of the
Wentzels cross section is used. To solve the scattering problem of heavy
ions it is necessary to introduce an effective particle whose mass is equal to
the relativistic reduced mass of the system defined as
m1 m2 c2
r . (6.80)
Ecm
where m1 and m2 are incident and target rest masses respectively and Ecm
(in Eq. (6.71) M1,2 = Ecm /c2 ) is the total center of mass energy of the
two particles system. The effective particle interacts with a fixed scattering
center with interacting potential expressed by Eq. (6.64) . The momentum
of the effective particle is equal to the momentum of the incoming particle
calculated in the center of mass system (pr p1cm ). Since the target particle
is inside the material it can be considered at rest in the laboratory as a
consequence the magnitude of pr is calculated as
m2 c2
pr p1cm = p1lab , (6.81)
Ecm
with Ecm given by
p
Ecm = (m1 c2 )2 + (m2 c2 )2 + 2E1lab m2 c2 , (6.82)
where p1lab and E1lab are the momentum and the total energy of the incom-
ing particle in the laboratory system respectively. The velocity (r ) of the
effective particle is obtained by the relation
!2
1 r c2
=1+ . (6.83)
r2 pr c

82
The modified Wentzels cross section is then equal to:
2
Z1 Z2 e2

d(r ) 1
= (6.84)
d pr c r (2As + 1 cos r )2
(in Eq. (6.72) pin pr ) where Z1 and Z2 are the nuclear proton numbers
of projectile and of target respectively; As is the screening coefficient [see
Eq. (6.73)] and r is the scattering angle of the effective particle which is
equal the one in the center of mass system (r 1cm ). Knowing the scat-
tering angle the recoil kinetic energy of the target particle after scattering is
calculated by ! 2
2 p1lab c
T = m2 c (1 cos r ). (6.85)
Ecm
The momentum and the total energy of the incident particle after scattering
in the laboratory system are obtained by the usual Lorentzs transformations.

6.5.6 Implementation Details

In the G4IonCoulombScatteringModel the scattering off electrons is not con-
sidered: only scattering off nuclei is simulated. Secondary particles are gen-
erated when T [Eq. (6.85)] is greater then a given threshold for displacement
Td ; it is not cut in range. The user can set this energy threshold Td by the
method SetRecoilThreshold(G4double Td ). The default screening coefficient
As is given by Eq. (6.73). If the user wants to use the one given by Eq. (6.77)
the condition SetHeavyIonCorr(1) must be set. When Z1 = 1 the Thomas-
Fermi screening length [aT F see Eq. (6.66)] is used in the calculation of As .
For Z1 2 the screening length is the universal one [aU see Eq. (6.67)].
In the G4IonCoulombCrossSection the total differential cross section is ob-
tained by the method NuclearCrossSection() where the Eq. (6.84) is inte-
grated in the interval (0, ):
2
Z1 Z2 e2

1
= (6.86)
pr c r As (As + 1)
The cosine of the scattering angle is chosen randomly in the interval (-1, 1)
according to the distribution of the total cross section and it is given by the
method SampleCosineTheta() which returns (1 cos r ).

6.5.7 Status of This Document

02.12.10 created by C. Consolandi and P.G. Rancoita
10.12.10 minor edition by V. Ivanchenko

83
Bibliography
[1] M. Boschini et al., Nuclear and Non-Ionizing Energy-Loss for Coulomb
Scattered Particles from Low Energy up to Relativistic Regime in
Space Radiation Environment, Proc. of the ICATPP Conference on
Cosmic Rays for Particle and Astroparticle Physics, October 78
2010, Villa Olmo, Como, Italy, World Scientific, Singapore (2011);
arXiv:1011.4822v3 [[Link]-ph], available at the web site:
[Link]

[2] C. Leroy and P.G. Rancoita, Principles of Radiation Interaction in Mat-

ter and Detection, 2nd Edition, World Scientific (Singapore) 2009.

[3] ICRU, Stopping Powers and Ranges for Protons and Alpha Particles.
ICRU Report 49, 1993.

[4] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping Range of Ions in
Solids, Vol. 1, Pergamon Press (New York) 1985.

[5] G. Wentzel, Z. Phys. 40 (1926), 590593.

[6] von G. Moli`ere, Z. Naturforsh. A2 (1947), 133145; A3 (1948), 78.

[7] H.A. Bethe, Phys. Rev. 98 (1953), 12561266.

[8] J.F. Ziegler, M.D. Ziegler, J.P. Biersack, The Stopping and Range
of Ions in Matter, SRIM version 2008.03 (2008), available at:
[Link]

[9] J.F. Ziegler, M.D. Ziegler, J.P. Biersack, The Stopping and Range of
Ions in Matter, SRIM Co. (Chester.) 2008.

[10] C. Leroy and P.G. Rancoita, Reports on Progress in Physics 70, 4 (2007)
493625.

[11] S.R. Messenger et al., IEEE Trans. on Nucl. Sci. 50 (2003), 19191923.

[12] I. Jun et al., IEEE Trans. on Nucl. Sci. 50 (2003) 19241928.

84
6.6 Electron Screened Single Scattering and
NIEL
The present treatment[1] of electronnucleus interaction is based on numer-
ical and analytical approximations of the Mott differential cross section. It
accounts for effects due to screened Coulomb potentials, finite sizes and finite
rest masses of nuclei for electron with kinetic energies above 200 keV and up
to ultra high. This treatment allows one to determine both the total and
differential cross sections, thus, to calculate the resulting nuclear and non-
ionizing stopping powers (NIEL). Above a few hundreds of MeV, neglecting
the effects of finite sizes and rest masses of recoil nuclei the stopping power
and NIEL result to be largely underestimated, while, above a few tens of
MeV prevents a further large increase, thus, resulting in approaching almost
constant values at high energies.
The non-ionizing energy-loss (NIEL) is the energy lost from a particle
traversing a unit length of a medium through physical process resulting in
permanent displacement damages (e.g. see Ref.[2]). The nuclear stopping
power and NIEL deposition - due to elastic Coulomb scatterings - from pro-
tons, light- and heavy-ions traversing an absorber were previously dealt[3]
and is available in Geant4 (6.5) (see also Sections 1.6, 1.6.1, [Link].4.2,
[Link] of Ref.[4]). In the present model included in GEANT4, the nuclear
stopping power and NIEL deposition due to elastic Coulomb scatterings of
electrons are treated up to ultra relativistic energies.

6.6.1 Scattering Cross Section of Electrons on Nuclei

The scattering of electrons by unscreened atomic nuclei was treated by Mott
extending a method - dealing with incident and scattered waves on point-like
nuclei - of Wentzel and including effects related to the spin of electrons. The
differential cross section (DCS) - the so-called Mott differential cross section
(MDCS) - was expressed by Mott as two conditionally convergent infinite
series in terms of Legendre expansions. In MottWentzel treatment, the
scattering occurs on a field of force generating a radially dependent Coulomb
- unscreened (screened) in Mott (Wentzel) - potential. Furthermore, the
MDCS was derived in the laboratory reference system for infinitely heavy
nuclei initially at rest with negligible spin effects and must be numerically
evaluated for any specific nuclear target. Effects related to the recoil and
finite rest mass of the target nucleus (M ) were neglected. Thus, in this
framework the total energy of electrons has to be smaller or much smaller
than M c2 .

85
 The MDCS is usually expressed as:

d Mott () d Rut Mott

= R , (6.87)
d d
where RMott is the ratio between the MDCS and Rutherfords formula [RDCS,
see Equation (1) of Ref.[1]]. For electrons with kinetic energies from several
keV up to 900 MeV and target nuclei with 1 6 Z 6 90, Lijian, Quing and
Zhengming[5] provided a practical interpolated expression [Eq. (6.99)] for
RMott with an average error less than 1%; in the present treatment, that ex-
pression - discussed in Sect. 6.6.1 - is the one assumed for RMott in Eq. (6.87)
hereafter. The analytical expression derived by McKinley and Feshbach[6]
for the ratio with respect to Rutherfords formula [Equation (7) of Ref.[6]] is
given by:

RMcF = 1 2 sin2 (/2) + Z sin(/2) [1 sin(/2)] (6.88)

with the corresponding differential cross section (McFDCS)

d McF d Rut McF

= R . (6.89)
d d
Furthermore, for M c2 much larger than the total energy of incoming electron
energies the distinction between laboratory (i.e., the system in which the tar-
get particle is initially at rest) and center-of-mass (CoM) systems disappears
(e.g., see discussion in Section 1.6.1 of Ref.[4]). Furthermore, in the CoM of
the reaction the energy transferred from an electron to a nucleus initially at
rest in the laboratory system (i.e., its recoil kinetic energy T ) is related with
the maximum energy transferable Tmax as

T = Tmax sin2 ( /2) (6.90)

[e.g., see Equations (1.27, 1.95) at page 11 and 31, respectively, of Ref.[4]],
where is the scattering angle in the CoM system. In addition, one obtains
Tmax
dT = d . (6.91)
4
Since for M c2 much larger than the electron energy is , one finds that
Eq. (6.90) can be approximated as

T Tmax sin2 (/2) , (6.92)

T
= sin2 (/2) = (6.93)
Tmax

86
and
Tmax
dT d. (6.94)
4
Using Eqs. (6.88, 6.93, 6.94), Rutherfords formula and Eq. (6.89) can be
respectively rewritten as:
 22
d Rut Ze Tmax
= = , (6.95)
dT pc T2
 22 " r #
d McF Ze Tmax T T
= = 1 ( +Z)+Z (6.96)
T pc T2 Tmax Tmax
 22
Ze Tmax McF
= R (T )
pc T2

with " r #
T T
RMcF (T ) = 1 ( +Z)+Z . (6.97)
Tmax Tmax
Finally, in a similar way the MDCS [Eq. (6.87)] is

d Mott (T ) d Rut Mott

= R (T )
dT dT
 22
Ze Tmax Mott
= R (T ) (6.98)
pc T2

with RMott (T ) from Eq. (6.101).

Interpolated Expression for RMott

Recently, Lijian, Quing and Zhengming[5] provided a practical interpolated
expression [Eq. (6.99)] which is a function of both and for electron energies
from several keV up to 900 MeV, i.e.,
4
X
Mott
R = aj (Z, )(1 cos )j/2 , (6.99)
j=0

where
6
X
aj (Z, ) = bk,j (Z)( )k1 , (6.100)
k=1

and c = 0.7181287 c is the mean velocity of electrons within the above

mentioned energy range. The coefficients bk,j (Z) are listed in Table 1 of

87
 2.2

2.0

1.8

1.6 Pb

1.4

RMott 1.2
Fe

1.0
Si

0.8
Li

0.6

0.4

0.2

0.0
0 20 40 60 80 100 120 140 160 180

Scattering Angle [degree]

Figure 6.4: RMott obtained from Eq. (6.99) at 100 MeV for Li, Si, Fe and Pb
nuclei as a function of scattering angle.

Ref.[5] for 1 6 Z 6 90. RMott obtained from Eq. (6.99) at 100 MeV is shown
in Fig. 6.4 for Li, Si, Fe and Pb nuclei as a function of scattering angle.
Furthermore, it has to be remarked that the energy dependence of RMott
from Eq. (6.99) was studied and observed to be negligible above 10 MeV
[for instance, see Eq. (6.100)].
Finally, from Eqs.(6.90, 6.99) [e.g., see also Equation (1.93) at page 31
of Ref.[4]], one finds that RMott can be expressed in terms of the transferred
energy T as
4  j/2
Mott
X 2T
R (T ) = aj (Z, ) . (6.101)
j=0
T max

Screened Coulomb Potentials

The simple scattering model due to Wentzel - with a single exponential
screening function [e.g., see Equation (2.71) at page 95 of Ref.[4]] - was re-
peatedly employed in treating single and multiple Coulomb scattering with
screened potentials. Neglecting effects like those related to spin and finite
size of nuclei, for proton and nucleus interactions on nuclei it was shown
that the resulting elastic differential cross section of a projectile with bare
nuclear-charge ez on a target with bare nuclear-charge eZ differs from the
Rutherford differential cross section (RDCS) by an additional term - the
so-called screening parameter - which prevents the divergence of the cross
section when the angle of scattered particles approaches 0 [e.g., see Sec-
tion 1.6.1 of Ref.[4]]. For z = 1 particles the screening parameter As,M is

88
expressed as
 2 "  2 #
~ Z
As,M = 1.13 + 3.76 (6.102)
2 p aTF

where , c and ~ are the fine-structure constant, speed of light and reduced
Planck constant, respectively; p (c) is the momentum (velocity) of the in-
coming particle undergoing the scattering onto a target supposed to be ini-
tially at rest - i.e., in the laboratory system -; aTF is the screening length
suggested by ThomasFermi
CTF a0
aTF = (6.103)
Z 1/3
with
~2
a0 =
me2
the Bohr radius, m the electron rest mass and
 2/3
1 3
CTF = 0.88534
2 4

a constant introduced in the ThomasFermi model [e.g., see Ref.[3] , Equa-

tions (2.73, 2,82) - at page 95 and 99, respectively - of Ref.[4], see also
references therein]. The modified Rutherfords formula [d WM ()/d], i.e.,
the differential cross section - obtained from the WentzelMoli`ere treatment
of the single scattering on screened nuclear potential - is given by [e.g., see
Equation (2.84) of Ref.[4] and Ref.[3], see also references therein]:
2
d WM () zZe2

1
= 2 (6.104)
d 2 p c As,M + sin2 (/2)


d Rut 2
= F (). (6.105)
d
with
sin2 (/2)
F() = . (6.106)
As,M + sin2 (/2)
F() - the so-called screening factor - depends on the scattering angle and
the screening parameter As,M . As discussed in Sect. 6.6.1, the term As,M
(the screening parameter) cannot be neglected in the DCS [Eq. (6.105)] for
scattering angles () within a forward (with respect to the electron direction)
angular region narrowing with increasing energy from several degrees (for

89
high-Z material) at 200 keV down to less than or much less than a mrad
above 200 MeV.
An approximated description of elastic interactions of electrons with screened
Coulomb fields of nuclei can be obtained by the factorization of the MDCS,
i.e., involving Rutherfords formula [d Rut /d] for particle with z = 1, the
screening factor [F()] and the ratio RMott between the RDCS and MDCS:
Mott
dsc () d Rut 2
F () RMott . (6.107)
d d
Thus, the corresponding screened differential cross section derived using the
analytical expression from McKinley and Feshbach[6] can be approximated
with
McF
dsc () d Rut 2
F () RMcF . (6.108)
d d
Zeitler and Olsen[7] suggested that for electron energies above 200 keV the
overlap of spin and screening effects is small for all elements and for all
energies; for lower energies the overlapping of the spin and screening effects
may be appreciable for heavy elements and large angles.

Finite Nuclear Size

The ratio between the actual measured and that expected from the point-
like differential cross section expresses the square of nuclear form factor (|F |)
which, in turn, depends on the momentum transfer q, i.e., that acquired by
the target initially at rest:
p
T (T + 2M c2 )
q= , (6.109)
c
with T from Eq. (6.90) or for M c2 larger or much larger than the electron
energy from its approximate expression Eq. (6.92).
The approximated (factorized) differential cross section for elastic inter-
actions of electrons with screened Coulomb fields of nuclei [Eq. (6.107)] ac-
counting for the effects due to the finite nuclear size is given by:
Mott
dsc,F () d Rut 2
F () RMott |F (q)|2 . (6.110)
d d
Thus, using the analytical expression derived by McKinley and Feshbach[6]
[Eq. (6.88)] one obtains that the corresponding screened differential cross

90
section [Eq. (6.108)] accounting for the finite nuclear size effects
McF
dsc,F () d Rut 2
F () RMcF |F (q)|2 (6.111)
d d
d Rut 2
= F () |F (q)|2
d
1 2 sin2 (/2) + Z sin(/2) [1 sin(/2)](6.112)

.

In terms of kinetic energy, one can respectively rewrite Eqs. (6.110, 6.111) as
Mott
dsc,F (T ) d Rut 2
= F (T ) RMott (T ) |F (q)|2 (6.113)
dT dT
McF
dsc,F (T ) d Rut (T ) 2
F (T ) RMcF (T ) |F (q)|2 (6.114)
dT dT
with d Rut /dT from Eq. (6.95), RMott (T ) from Eq. (6.101), RMcF (T ) from
Eq. (6.97) and, using Eqs. (6.90, 6.92, 6.106),
T
F(T ) = .
Tmax As,M + T

For instance, the form factor Fexp is

 2
1  qrn 2
Fexp (q) = 1 + , (6.115)
12 ~

where rn is the nuclear radius, rn can be parameterized by

rn = 1.27A0.27 fm (6.116)

with A the atomic weight. Equation (6.116) provides values of rn in agreement

up to heavy nuclei (like Pb and U) with those available, for instance, in
Table 1 of Ref.[8] .

Finite Rest Mass of Target Nucleus

The DCS treated in Sects. [Link].1 is based on the extension of MDCS to
include effects due to interactions on screened Coulomb potentials of nuclei
and their finite size. However, the electron energies were considered small (or
much smaller) with respect to that (M c2 ) corresponding to rest mass (M )
target nuclei.
The Rutherford scattering on screened Coulomb fields - i.e., under the
action of a central forces - by massive charged particles at energies large or

91
much larger than M c2 was treated by Boschini et al.[3] in the CoM system
(e.g., see also Sections 1.6, 1.6.1, [Link] of Ref.[4] and references therein).
It was shown that the differential cross section [d WM ( )/d with the
scattering angle in the CoM system] is that one derived for describing the in-
teraction on a fixed scattering center of a particle with i) momentum pr equal
to the momentum of the incoming particle (i.e., the electron in the present
treatment) in the CoM system and ii) rest mass equal to the relativistic re-
duced mass rel [e.g., see Equations (1.80, 1.81) at page 28 of Ref.[4]]. rel is
given by
mM
rel = (6.117)
M1,2
mM c
= q p , (6.118)
m2 c2 + M 2 c2 + 2 M m2 c4 + p2 c2

where p is the momentum of the incoming particle (the electron in the present
treatment) in the laboratory system: m is the rest mass of the incoming
particle (i.e., the electron rest mass); finally, M1,2 is the invariant mass - e.g.,
Section 1.3.2 of Ref.[4] - of the two-particle system. Thus, the velocity of the
interacting particle is [e.g., see Equation (1.82) at page 29 of Ref.[4]]
v"
u  2 #1
u rel c
r c = ct 1 + . (6.119)
pr

For an incoming particle with z = 1, d WM ( )/d is given by

2
d WM ( ) Ze2

1
= 2 , (6.120)
d 2 pr r c As + sin2 ( /2)


with  2 "  2 #
~ Z
As =
1.13 + 3.76 (6.121)
2 pr aTF r
the screening factor [e.g., see Equations (2.87, 2.88) at page 103 of Ref.[4]].
Equation (6.120) can be rewritten as

d WM ( ) d Rut ( ) 2
= FCoM ( ) (6.122)
d d
with 2
d Rut ( ) Ze2

1

= 4 (6.123)
d 2pr r c sin ( /2)

92
the corresponding RDCS for the reaction in the CoM system [e.g., see Equa-
tion (1.79) at page 28 of Ref.[4]] and

sin2 ( /2)
FCoM ( ) = (6.124)
As + sin2 ( /2)

the screening factor. Using, Eqs. (6.90, 6.91), one can respectively rewrite
Eqs. (6.123, 6.124, 6.122, 6.120) as
2
d Rut Ze2

Tmax
= (6.125)
dT pr r c T2
T
FCoM (T ) = (6.126)
Tmax As + T

d WM (T ) d Rut
= FCoM (T ) (6.127)
dT dT
 2 2
d WM (T ) Ze Tmax
= . (6.128)
dT pr r c (Tmax As + T )2

[e.g., see Equation (2.90) at page 103 of Ref.[4] or Equation (13) of Ref.[3]].
To account for the finite rest mass of target nucleus the factorized MDCS
[Eq. (6.110)] has to be re-expressed in the CoM system using as:

( ) d Rut ( ) 2
Mott
dsc,F,CoM 2


FCoM ( ) RMott
CoM ( ) |F (q)| , (6.129)
d d
where F (q) is the nuclear form factor (Sect. 6.6.1) with q the momentum
transfer to the recoil nucleus [Eq. (6.109)]; finally, as discussed in Sect. 6.6.1,
RMott exhibits almost no dependence on electron energy above 10 MeV,
thus, since at low energies and r , RMott
CoM ( ) is obtained replacing
and r with and r , respectively, in Eq. (6.99).
Using the analytical expression derived by McKinley and Feshbach[6] , one
finds that the corresponding screened differential cross section accounting for
the finite nuclear size effects [Eqs. (6.111, 6.112)] can be re-expressed as

( )
McF
dsc,F,CoM d Rut ( ) 2 2
FCoM ( ) RMcF
CoM ( ) |F (q)| (6.130)
d d
with
2
RMcF 2

CoM ( ) = 1r sin ( /2)+Z r sin( /2) [1sin( /2)] . (6.131)

93
 -4
5.0x10

g ]
-1
2
Nuclear stopping power [MeV cm
-4
4.0x10 12
C

28
Si

56
-4 Fe
3.0x10
7
Li

-4
2.0x10

0 1 2 3 4 5 6
10 10 10 10 10 10 10
Kinetic Energy [MeV]

Figure 6.5: In MeV cm2 /g, nuclear stopping powers in 7 Li, 12 C, 28 Si and 56 Fe
- calculated from Eq. (6.136) - and divided by the density of the material as
a function of the kinetic energy of electrons from 200 keV up to 1 TeV.

In terms of kinetic energy T , from Eqs. (6.90, 6.91) one can respectively
rewrite Eqs. (6.129, 6.130) as
Mott
dsc,F,CoM (T ) d Rut 2 2
= FCoM (T ) RMott
CoM (T ) |F (q)| (6.132)
dT dT
McF
dsc,F,CoM (T ) d Rut (T ) 2 2
FCoM (T ) RMcF
CoM (T ) |F (q)| (6.133)
dT dT

with d Rut /dT from Eq. (6.125), FCoM (T ) from Eq. (6.126) and RMcF
CoM (T )

replacing with r in Eq. (6.97), i.e.,
" r #
T T
RMcF
CoM (T ) = 1r ( +Z)+Zr . (6.134)
Tmax r Tmax

Finally, as discussed in Sect. 6.6.1, RMott (T ) exhibits almost no dependence

on electron energy above 10 MeV, thus, since at low energies and
r , RMott
CoM (T ) is obtained replacing with r in Eq. (6.101).

6.6.2 Nuclear Stopping Power of Electrons

Using Eq. (6.132), the nuclear stopping power - in MeV cm1 - of Coulomb
electronnucleus interaction can be obtained as
 Mott Z Tmax Mott
dE dsc,F,CoM (T )
= nA T dT (6.135)
dx nucl 0 dT
with nA the number of nuclei (atoms) per unit of volume [e.g., see Equa-
tion (1.71) of Ref.[4]] and, finally, the negative sign indicates that the energy

94
is lost by the electron (thus, achieved by recoil targets). Using the analytical
approximation derived by McKinley and Feshbach[6], i.e., Eq. (6.133), for
the nuclear stopping power one finds
McF Tmax McF
dsc,F,CoM (T )

dE
Z
= nA T dT. (6.136)
dx nucl 0 dT

As already mentioned in Sect. 6.6.1, the large momentum transfers -

corresponding to large scattering angles - are disfavored by effects due to
the finite nuclear size accounted for by means of the nuclear form factor
(Sect.6.6.1). For instance, the ratios of nuclear stopping powers of electrons
in silicon are shown in Ref.[1] as a function of the kinetic energies of electrons
from 200 keV up to 1 TeV. These ratios are the nuclear stopping powers
calculated neglecting i) nuclear size effects (i.e., for |Fexp |2 = 1) and ii) effects
due to the finite rest mass of the target nucleus [i.e., in Eq. (6.136) replacing
McF McF
dsc,F,CoM (T )/dT with dsc,F (T )/dT from Eq. (6.114)] both divided by that
one obtained using Eq. (6.136). Above a few tens of MeV, a larger stopping
power is found assuming |Fexp |2 = 1 and, in addition, above a few hundreds
of MeV the stopping power largely decreases when the effects of nuclear rest
mass are not accounted for.
In Fig. 6.5 , the nuclear stopping powers in 7 Li, 12 C, 28 Si and 56 Fe are
shown as a function of the kinetic energy of electrons from 200 keV up to
1 TeV. These nuclear stopping powers in MeV cm2 /g are calculated from
Eq. (6.136) and divided by the density of the medium.

6.6.3 Non-Ionizing Energy-Loss of Electrons

In case of Coulomb scattering of electrons on nuclei, the non-ionizing energy-
loss can be calculated using (as discussed in Sect. [Link].2) the MDCRS or
its approximate expression McFDCS [e.g., Eqs. (6.132, 6.133), respectively],
once the screened Coulomb fields, finite sizes and rest masses of nuclei are
accounted for, i.e., in Mev/cm
NIEL Tmax Mott
dsc,F,CoM (T )

dE
Z
= nA T L(T ) dT (6.137)
dx n,Mott Td dT
or NIEL Tmax McF
dsc,F,CoM (T )

dE
Z
= nA T L(T ) dT (6.138)
dx n,McF Td dT
[e.g., see Equation (4.113) at page 402 and, in addition, Sections [Link].1.2
of Ref.[4]], where T is the kinetic energy transferred to the target nucleus,

95
L(T ) is the fraction of T deposited by means of displacement processes. The
Lindhard partition function, L(T ), can be approximated using the so-called
NorgettRobintsonTorrens expression [e.g., see Equations (4.121, 4.123) at
pages 404 and 405, respectively, of Ref.[4] (see also references therein)]. Tde =
T L(T ) is the so-called damage energy, i.e., the energy deposited by a recoil
nucleus with kinetic energy T via displacement damages inside the medium. In
Eqs. (6.137, 6.138) the integral is computed from the minimum energy Td -
the so-called threshold energy for displacement, i.e., that energy necessary
to displace the atom from its lattice position - up to the maximum energy
Tmax that can be transferred during a single collision process. For instance,
Td is about 21 eV in silicon requiring electrons with kinetic energies above
220 kev.
As already discussed with respect to nuclear stopping powers in Sect. 6.6.2,
the large momentum transfers (corresponding to large scattering angles) are
disfavored by effects due to the finite nuclear size accounted for by the nu-
clear form factor. For instance, the ratios of NIELs for electrons in silicon are
shown in Ref.[1] as a function of the kinetic energy of electrons from 220 keV
up to 1 TeV. These ratios are the NIELs calculated neglecting i) nuclear size
effects (i.e., for |Fexp |2 = 1) and ii) effects due to the finite rest mass of the tar-
McF McF
get nucleus [i.e., in Eq. (6.138) replacing dsc,F,CoM (T )/dT with dsc,F (T )/dT
from Eq. (6.114)] both divided by that one obtained using Eq. (6.138). Above
10 MeV, the NIEL is 20% larger assuming |Fexp |2 = 1 and, in addition,
above (100200) MeV the calculated NIEL largely decreases when the effects
of nuclear rest mass are not accounted for.

6.7 G4eSingleScatteringModel
The G4eSingleScatteringModel performs the single scattering interaction of
electrons on nuclei. The differential cross section (DCS) for the energy trans-
ferred is define in the G4ScreeningMottCrossSection class. In this class the
analytical McKinley and Feshbach [6] Mott differential cross Section approx-
imation is implemented. This CDS is modified by the introduction of the
Molieres [9] screening coefficient. In addition the exponential charge distri-
bution Nuclear Form Factor is applied [10]. This treatment is fully performed
in the center of mass system and the usual Lorentz transformations are ap-
plied to obtained the energy and momentum quantities in the laboratory
system after scattering. This model well simulates the interacting process
for low scattering angles and it is suitable for high energy electrons (from
200 keV) incident on medium light target nuclei. The nuclear energy loss
(i.e. nuclear stopping power) is calculated for every single interaction. In

96
addition the production of secondary scattered nuclei is simulated from a
threshold kinetic energy which can be decided by the user (threshold energy
for displacement).

6.7.1 The method

In the G4eSingleScatteringModel the method ComputeCrossSectionPerAtom()
performs the total cross section computation. The SetupParticle() and the
DefineMaterial() methods are called to defined the incident and target par-
ticles. Before the total cross section computation, the SetupKinematic()
method of the G4ScreeningMottCrossSection class calculates all the physi-
cal quantities in the center of mass system (CM). The scattering in the CM
system is equivalent to the one of an effective particle which interacts with
a fixed scattering center. The effective particle rest mass is equal to the
relativistic reduced mass of the system whose expression is calculated by:
M c2
=m (6.139)
Ecm
where m and M are rest masses of the electron and of the target nuclei
respectively. Ecm is the total center of mass energy and, since the target is
at rest before scattering, its expression is calculated by:
p
Ecm = (mc2 )2 + (M c2 )2 + 2E M c2 (6.140)
where E = mc2 is the total energy of the electron before scattering in the
laboratory system. The momentum and the scattering angle of the effective
particle are equal to the corresponding quantities calculated in the center of
mass system (p pcm , cm ) of the incident electron:
M c2
pc = p c (6.141)
Ecm
where p is the momentum of the incident electron calculated in the labora-
tory system. The velocity of the effective particle is related with its momen-
tum by the following expression:
1  c2 2
= 1 + (6.142)
2 pc
The integration of the DCS is performed by the NuclearCrossSection() method
of the G4ScreeningMottCrossSection:
Z max
d()
tot = 2 sin d (6.143)
min d

97
The integration is performed in the scattering range [0 ;] but the user can
decide to vary the minimum (min ) and the maximum (max ) scattering an-
gles. The DCS is then given by:
!2
d() Ze2 RM cF |FN (q)|2
=
2 (6.144)
d c2 2 2As + 2 sin2 (/2)
where Z is the atomic number of the nucleus, As is the screening coefficient
whose expression has been given by Moliere[9] :
 2   2 
~ Z
As = 1.13 + 3.76 (6.145)
2p aT F
where aT F is the Thomas-Fermi screening length given by:
0.88534 a0
aT F = (6.146)
Z 1/3
and a0 is the Bhor radius. RM cF is the ratio of the Mott to the Rutherfor
DCS given by McKinley and Feshbach approximation [6]:
 
2 2


RM cF = 1 sin (/2) + Z sin(/2) 1 sin(/2) (6.147)

The nuclear form factor for the exponential charge distribution is given by
[10]:
" #2
(qRN )2
FN (q) = 1 + (6.148)
12~2
where RN is the nuclear radius that is parameterized by:
RN = 1.27A0.27 fm. (6.149)
q is the momentum transferred to the nucleus and it is calculated as:
p
qc = T (T + 2M c2 ) (6.150)
where T is the kinetic energy transferred to the nucleus. This kinetic energy
is calculated in the GetNewDirection() method as:
2M c2 (p c)2
T = 2
sin2 /2. (6.151)
Ecm
The scattering angle calculation is performed in the GetScatteringAngle()
method of G4ScreeningMottCrossSection class. By means of AngleDistribu-
tion() function the scattering angle is chosen randomly according to the total

98
cross section distribution (p.d.f. probability density function) by means of
the inverse transform method.
In the SampleSecondary() method of G4eSingleScatteringModel the ki-
netic energy of the incident particle after scattering is then calculated as

Enew = E T where E is the electron incident kinetic energy (in lab.); in
addition the new particle direction and momentum are obtained from the
scattering angle information.

6.7.2 Implementation Details

The scattering angle probability density function f () (p.d.f.) is performed
by the AngleDistribution() of G4ScreeningMottCrossSection class where the
inverse transform method is applied. The normalized cumulative function of
the cross section is calculated as a function of the scattering angle in this
way:
2 d(t)
Z Z
n () f ()d = sin tdt (6.152)
tot 0 d
The normalized cumulative function n () depends on the DCS and its val-
ues range in the interval [0;1]. After this calculation a random number r,
uniformly distributed in the same interval [0;1], is chosen in order to fix the
cumulative function value (i.e. r n ()). This number is the probability
to find the scattering angle in the interval [; + d]. The scattering angle
is then given by the inverse function of n ().

The threshold energy for displacement Th can by set by the user in her/his
own Physic class by adding the electromagnetic model:
G4eSingleCoulombScatteringModel* mod=
new G4eSingleCoulombScatteringModel();
mod->SetRecoilThreshold(Th);

If the energy lost by the incident particle is grater then this threshold value
a new secondary particle is created for transportation processes. The energy
lost is added to ProposeNonIonizingEnergyDeposit().

In next patches the NIEL and the Lijians et al. Motts approximations[5]
calculations will be included in the G4eSingleCoulombScatteringModel.

6.8 Status of this Document

17.11.11 created by C. Consolandi and P.G. Rancoita

99
Bibliography
[1] M. Boschini et al., Nuclear and Non-Ionizing Energy-Loss of Electrons
with Low and Relativistic Energies in Materials and Space Environ-
ment,Proc. of the ICATPP Conference on Cosmic Rays for Particle and
Astroparticle Physics, October 37 (2011), Villa Olmo, Como, Italy, S.
Giani, C. Leroy, L. Price, P.G. Rancoita and R. Ruchri, Editors, World
Scientific, Singapore (2012); arXiv 1111.4042.

[2] C. Leroy and P.G. Rancoita, Particle Interaction and Displacement

Damage in Silicon Devices operated in Radiation Environments,
Rep. Prog. in Phys. 70 (no. 4)(2007), 403625, doi: 10.1088/0034-
4885/70/4/R01.

[3] M. Boschini et al., Nuclear and Non-Ionizing Energy-Loss for Coulomb

Scattered Particles from Low Energy up to Relativistic Regime in Space
Radiation Environment, Proc. of the ICATPP Conference on Cos-
mic Rays for Particle and Astroparticle Physics, October 78 (2010),
Villa Olmo, Como, Italy, S. Giani, C. Leroy and P.G. Rancoita, Edi-
tors, World Scientific, Singapore (2011), 923, IBSN: 978-981-4329-02-6;
arXiv 1011.4822.

[4] C. Leroy and P.G. Rancoita, Principles of Radiation Interaction in Mat-

ter and Detection, 3rd Edition, World Scientific (Singapore) 2011.

[5] T. Lijian, H. Quing and L. Zhengming, Radiat. Phys. Chem. 45 (1995),

235245.

[6] [Link]. McKinley and H. Feshbach, Phys. Rev. 74 (1948), 17591763.

[7] E. Zeitler and A. Olsen, Phys. Rev. 136 (1956), A1546-A1552.

[8] H. De Vries, C.W. De Jager, and C. De Vries, Atomic Data and Nuclear
Data Tables 36 (1987), 495.

[9] von G. Moliere, Z. Naturforsh A2 (1947), 133-145; A3 (1948), 78.

[10] A.V. Butkevich et al. Nucl. Instr. and Meth. in Phys. Res. A 488 (2002),
282-294.

100
Chapter 7

Energy loss of Charged

Particles

101
7.1 Mean Energy Loss
Energy loss processes are very similar for e + /e , + / and charged
hadrons, so a common description for them was a natural choice in Geant4
[1], [2]. Any energy loss process must calculate the continuous and discrete
energy loss in a material. Below a given energy threshold the energy loss is
continuous and above it the energy loss is simulated by the explicit produc-
tion of secondary particles - gammas, electrons, and positrons.

7.1.1 Method
Let
d(Z, E, T )
dT
be the differential cross-section per atom (atomic number Z) for the ejection
of a secondary particle with kinetic energy T by an incident particle of total
energy E moving in a material of density . The value of the kinetic energy
cut-off or production threshold is denoted by Tcut . Below this threshold the
soft secondaries ejected are simulated as continuous energy loss by the inci-
dent particle, and above it they are explicitly generated. The mean rate of
energy loss is given by:
Z Tcut
dEsof t (E, Tcut ) d(Z, E, T )
= nat T dT (7.1)
dx 0 dT
where nat is the number of atoms per volume in the material. The total cross
section per atom for the ejection of a secondary of energy
T > Tcut is Z Tmax
d(Z, E, T )
(Z, E, Tcut ) = dT (7.2)
Tcut dT
where Tmax is the maximum energy transferable to the secondary particle.
If there are several processes providing energy loss for a given particle, then
the total continuous part of the energy loss is the sum:
tot
dEsof t (E, Tcut )
X dEsof t,i (E, Tcut )
= . (7.3)
dx i
dx

These values are pre-calculated during the initialization phase of Geant4

and stored in the dE/dx table. Using this table the ranges of the particle
in given materials are calculated and stored in the Range table. The Range
table is then inverted to provide the InverseRange table. At run time, values
of the particles continuous energy loss and range are obtained using these

102
tables. Concrete processes contributing to the energy loss are not involved in
the calculation at that moment. In contrast, the production of secondaries
with kinetic energies above the production threshold is sampled by each
concrete energy loss process.
The default energy interval for these tables extends from 100eV to 10T eV
and the default number of bins is 77. For muons and for heavy particles en-
ergy loss processes models are valid for higher energies and can be extended.
For muons uppper limit may be set to 1000P eV .

7.1.2 General Interfaces

There are a number of similar functions for discrete electromagnetic pro-
cesses and for electromagnetic (EM) packages an additional base classes were
designed to provide common computations [2]. Common calculations for dis-
crete EM processes are performed in the class G4V EnergyLossP rocess. De-
rived classes (7.1) are concrete processes providing initialisation. The physics
models are implemented using the G4V EmM odel interface. Each process
may have one or many models defined to be active over a given energy range
and set of G4Regions. Models are implementing computation of energy loss,
cross section and sampling of final state. The list of EM processes and models
for gamma incident is shown in Table 7.1.

7.1.3 Step-size Limit

Continuous energy loss imposes a limit on the step-size because of the energy
dependence of the cross sections. It is generally assumed in MC programs
(for example, Geant3) that the cross sections are approximately constant
along a step, i.e. the step size should be small enough, so that the change
in cross section along the step is also small. In principle one must use very
small steps in order to insure an accurate simulation, however the computing
time increases as the step-size decreases.
For EM processes the exact solution is available (see 7.3) but is is not
implemented yet for all physics processes including hadronics. A good com-
promise is to limit the step-size by not allowing the stopping range of the
particle to decrease by more than 20 % during the step. This condition
works well for particles with kinetic energies > 1 MeV, but for lower energies
it gives too short step-sizes, so must be relaxed. To solve this problem a
lower limit on the step-size was introduced. A smooth StepFunction, with
2 parameters, controls the step size. At high energy the maximum step
size is defined by Step/Range R (parameter dRoverRange). By default
R = 0.2. As the particle travels the maximum step size decreases gradually

103
 Table 7.1: List of process and model classes for charged particles.
EM process EM model Ref.
G4eIonisation G4MollerBhabhaModel 8.1
G4LivermoreIonisationModel 9.8
G4PenelopeIonisationModel 10.1.7
G4PAIModel 7.5
G4PAIPhotModel 7.5
G4ePolarizedIonisation G4PolarizedMollerBhabhaModel 17.1
G4MuIonisation G4MuBetheBlochModel 13.1
G4PAIModel 7.5
G4PAIPhotModel 7.5
G4hIonisation G4BetheBlochModel 12.1
G4BraggModel 12.1
G4ICRU73QOModel 12.2.1
G4PAIModel 7.5
G4PAIPhotModel 7.5
G4ionIonisation G4BetheBlochModel 12.1
G4BetheBlochIonGasModel 12.1
G4BraggIonModel 12.1
G4BraggIonGasModel 12.1
G4IonParametrisedLossModel 12.2.4
G4NuclearStopping G4ICRU49NuclearStoppingModel 12.1.3
G4mplIonisation G4mplIonisationWithDeltaModel
G4eBremsstrahlung G4SeltzerBergerModel 8.2.1
G4eBremsstrahlungRelModel 8.2.2
G4LivermoreBremsstrahlungModel 9.9
G4PenelopeBremsstrahlungModel 10.1.6
G4ePolarizedBremsstrahlung G4PolarizedBremsstrahlungModel 17.1
G4MuBremsstrahlung G4MuBremsstrahlungModel 13.2
G4hBremsstrahlung G4hBremsstrahlungModel
G4MuPairProduction G4MuPairProductionModel 13.3
G4hPairProduction G4hPairProductionModel

104
until the range becomes lower than R (parameter finalRange). Default final-
Range R = 1mm. For the case of a particle range R > R the StepFunction
provides limit for the step size Slim by the following formula:
 R 
Slim = R R + R (1 R ) 2 . (7.4)
R
In the opposite case of a small range Slim = R. The figure below shows the
ratio step/range as a function of range if step limitation is determined only
by the expression (7.4).

step

range

dRoverRange

finalRange

range

The parameters of StepFunction can be overwritten using an UI command:

/process/eLoss/StepFunction 0.2 1 mm
To provide more accurate simulation of particle ranges in physics constructors
G4EmStandardPhysics option3 and G4EmStandardPhysics option4 more strict
step limitation is chosen for different particle types.

7.1.4 Run Time Energy Loss Computation

The computation of the mean energy loss after a given step is done by using
the dE/dx, Range, and InverseRange tables. The dE/dx table is used if
the energy deposition (T ) is less than allowed limit T < T0 , where is

105
linearLossLimit parameter (by default = 0.01), T0 is the kinetic energy
of the particle. In that case
dE
T = s, (7.5)
dx
where T is the energy loss, s is the true step length. When a larger
percentage of energy is lost, the mean loss can be written as
T = T0 fT (r0 s) (7.6)
where r0 the range at the beginning of the step, the function fT (r) is the
inverse of the Range table (i.e. it gives the kinetic energy of the particle for
a range value of r). By default spline approximation is used to retrieve a
value from dE/dx, Range, and InverseRange tables. The spline flag can be
changed using an UI command:

/process/em/spline false
After the mean energy loss has been calculated, the process computes the
actual energy loss, i.e. the loss with fluctuations. The fluctuation models are
described in Section 7.2.
If deexcitation module (see 14.1) is enabled then simulation of atomic de-
excitation is performed using information on step length and ionisation cross
section. Fluorescence gamma and Auger electrons are produced above the
same threshold energy as -electrons and bremsstrahlung gammas. Following
UI commands can be used to enable atomic relaxation:
/process/em/deexcitation myregion true true true
/process/em/fluo true
/process/em/auger true
/process/em/pixe true
/process/em/deexcitationIgnoreCut true

The last command means that production threshold for electrons and gam-
mas are not checked, so full atomic de-excitation decay chain is simulated.

7.1.5 Energy Loss by Heavy Charged Particles

To save memory in the case of positively charged hadrons and ions energy
loss, dE/dx, Range and InverseRange tables are constructed only for pro-
ton, antiproton, muons, pions, kaons, and Generic Ion. The energy loss for
other particles is computed from these tables at the scaled kinetic energy
Tscaled :
Mbase
Tscaled = T , (7.7)
Mparticle

106
where T is the kinetic energy of the particle, Mbase and Mparticle are the masses
of the base particle (proton or kaon) and particle. For positively changed
hadrons with non-zero spin proton is used as a based particle, for negatively
charged hadrons with non-zero spin - antiproton, for charged particles with
zero spin - K + or K correspondingly. The virtual particle Generic Ion is
used as a base particle for for all ions with Z > 2. It has mass, change and
other quantum numbers of the proton. The energy loss can be defined via
scaling relation:
dE 2 dE
(T ) = qef f (F1 (T ) (Tscaled ) + F2 (T, qef f )), (7.8)
dx dx base
where qef f is particle effective change in units of positron charge, F1 and
F2 are correction function taking into account Birks effect, Block correction,
low-energy corrections based on data from evaluated data bases [5]. For a
hadron qef f is equal to the hadron charge, for a slow ion effective charge is
different from the charge of the ions nucleus, because of electron exchange
between transporting ion and the media. The effective charge approach is
used to describe this effect [3]. The scaling relation (7.7) is valid for any
combination of two heavy charged particles with accuracy corresponding to
high order mass, charge and spin corrections [4].

7.1.6 Status of This Document

09.10.98 created by L. Urban
01.12.03 revised by V. Ivanchenko
02.12.03 spelling and grammar check by D.H. Wright
09.12.05 minor update by V. Ivanchenko
14.06.07 added formula of StepFunction (M. Maire)
15.06.07 updated last sub-charter, list of processes and models by V. Ivanchenko
11.12.08 revised by V. Ivanchenko
09.12.10 revised by V. Ivanchenko

Bibliography
[1] S. Agostinelli et al., Geant4 a simulation toolkit Nucl. Instr. Meth.
A506 (2003) 250.

[2] J. Apostolakis et al., Geometry and physics of the Geant4 toolkit for high
and medium energy applications. Rad. Phys. Chem. 78 (2009) 859.

107
[3] J.F. Ziegler and J.M. Manoyan, Nucl. Instr. and Meth. B35 (1988) 215.

[4] ICRU (A. Allisy et al), Stopping Powers and Ranges for Protons and
Alpha Particles, ICRU Report 49, 1993.

[5] ICRU (R. Bimbot et al), Stopping of Ions Heavier than Helium, Journal
of the ICRU Vol5 No1 (2005) Report 73.

108
7.2 Energy Loss Fluctuations
The total continuous energy loss of charged particles is a stochastic quan-
tity with a distribution described in terms of a straggling function. The
straggling is partially taken into account in the simulation of energy loss
by the production of -electrons with energy T > Tcut (Eq.7.2). However,
continuous energy loss (Eq.7.1) also has fluctuations. Hence in the current
GEANT4 implementation different models of fluctuations implementing the
G4V EmF luctuationM odel interface:

G4BohrFluctuations;

G4IonFluctuations;

G4PAIModel;

G4PAIPhotModel;

G4UniversalFluctuation.

The last model is the default one used in main Physics List and will be de-
scribed below. Other models have limited applicability and will be described
in chapters for ion ionisation and PAI models.

7.2.1 Fluctuations in Thick Absorbers

The total continuous energy loss of charged particles is a stochastic quantity
with a distribution described in terms of a straggling function. The strag-
gling is partially taken into account in the simulation of energy loss by the
production of -electrons with energy T > Tc . However, continuous energy
loss also has fluctuations. Hence in the current GEANT4 implementation
two different models of fluctuations are applied depending on the value of
the parameter which is the lower limit of the number of interactions of the
particle in a step. The default value chosen is = 10. In the case of a high
range cut (i.e. energy loss without delta ray production) for thick absorbers
the following condition should be fulfilled:

E > Tmax (7.9)

where E is the mean continuous energy loss in a track segment of length s,

and Tmax is the maximum kinetic energy that can be transferred to the atomic
electron. If this condition holds the fluctuation of the total (unrestricted)
energy loss follows a Gaussian distribution. It is worth noting that this

109
condition can be true only for heavy particles, because for electrons, Tmax =
T /2, and for positrons, Tmax = T , where T is the kinetic energy of the
particle. In order to simulate the fluctuation of the continuous (restricted)
energy loss, the condition should be modified. After a study, the following
conditions have been chosen:
E > Tc (7.10)
Tmax <= 2 Tc (7.11)
where Tc is the cut kinetic energy of -electrons. For thick absorbers the
straggling function approaches the Gaussian distribution with Bohrs vari-
ance [4]:
Zh2 2
 
2 2 2
= 2re me c Nel 2 Tc s 1 , (7.12)
2
where re is the classical electron radius, Nel is the electron density of the
medium, Zh is the charge of the incident particle in units of positron charge,
and is the relativistic velocity.

7.2.2 Fluctuations in Thin Absorbers

If the conditions 7.10 and 7.11 are not satisfied the model of energy fluctua-
tions in thin absorbers is applied. The formulas used to compute the energy
loss fluctuation (straggling) are based on a very simple physics model of the
atom. It is assumed that the atoms have only two energy levels with binding
energies E1 and E2 . The particle-atom interaction can be an excitation with
energy loss E1 or E2 , or ionisation with energy loss distributed according to
a function g(E) 1/E 2 :
Z Tup
E0 Tup 1
g(E) dE = 1 = g(E) = . (7.13)
E0 Tup E0 E 2
The macroscopic cross section for excitation (i = 1, 2) is
fi ln[2mc2 ()2 /Ei ] 2
i = C (1 r) (7.14)
Ei ln[2mc2 ()2 /I] 2
and the ionisation cross section is
Tup E0
3 = C r (7.15)
E0 Tup ln( TEup0 )
where E0 denotes the ionisation energy of the atom, I is the mean ionisation
energy, Tup is the production threshold for delta ray production (or the max-
imum energy transfer if this value smaller than the production threshold),

110
Ei and fi are the energy levels and corresponding oscillator strengths of the
atom, and C and r are model parameters.
The oscillator strengths fi and energy levels Ei should satisfy the con-
straints
f1 + f2 = 1 (7.16)
f1 lnE1 + f2 lnE2 = lnI. (7.17)
The cross section formulas 7.14,7.15 and the sum rule equations 7.16,7.17
can be found e.g. in Ref. [1]. The model parameter C can be defined in the
following way. The numbers of collisions (ni , i = 1, 2 for excitation and 3 for
ionisation) follow the Poisson distribution with a mean value hni i. In a step
of length x the mean number of collisions is given by

hni i = x i (7.18)

The mean energy loss in a step is the sum of the excitation and ionisation
contributions and can be written as
 Z Tup 
dE
x = 1 E1 + 2 E2 + Eg(E)dE x. (7.19)
dx E0

From this, using Eq. 7.14 - 7.17, one can see that

C = dE/dx. (7.20)

The other parameters in the fluctuation model have been chosen in the follow-
ing way. Z f1 and Z f2 represent in the model the number of loosely/tightly
bound electrons
f2 = 0 f or Z = 1 (7.21)
f2 = 2/Z f or Z2 (7.22)
E2 = 10 eV Z 2 (7.23)
E0 = 10 eV . (7.24)
Using these parameter values, E2 corresponds approximately to the K-shell
energy of the atoms ( and Zf2 = 2 is the number of K-shell electrons).
The parameters f1 and E1 can be obtained from Eqs. 7.16 and 7.17. The
parameter r is the only variable in the model which can be tuned. This
parameter determines the relative contribution of ionisation and excitation to
the energy loss. Based on comparisons of simulated energy loss distributions
to experimental data, its value has been fixed as

r = 0.55 (7.25)

111
7.2.3 Width Correction Algorithm
This simple parametrization and sampling in the model give good values for
the most probable energy loss in thin layers. The width of the energy loss
distribution (Full Width at Half Maximum, FWHM) in most of the cases
is too small. In order to get good FWHM values a relatively simple width
correction algorithm has been applied. This algorithm rescales the energy
levels E1 , E2 and the number of excitations n1 , n2 in such a way that the
mean energy loss remains the same. Using this width correction scheme the
model gives not only good most probable energy loss, but good FWHM value
too.
Width correction algorithm is in the model since version 9.2. The updated
version in the model (in version 9.4) causes an important change in the
behaviour of the model: the results become much more stable, i.e. the results
do not change practically when the cuts and/or the stepsizes are changing.
Another important change: the (unphysical) second peak or shoulder in the
energy loss distribution which can be seen in some cases (energy loss in thin
gas layers) in older versions of the model disappeared. Limit of validity
of the model for thin targets: the model gives good (reliable) energy loss
distribution if the mean energy loss in the target is (f ew times) Iexc ,
where Iexc is the mean excitation energy of the target material.
This simple model of energy loss fluctuations is rather fast and can be
used for any thickness of material. This has been verified by performing
many simulations and comparing the results with experimental data, such as
those in Ref.[2]. As the limit of validity of Landaus theory is approached,
the loss distribution approaches the Landau form smoothly.

7.2.4 Sampling of Energy Loss

If the mean energy loss and step are in the range of validity of the Gaussian
approximation of the fluctuation (7.10 and 7.11), the Gaussian sampling is
used to compute the actual energy loss (7.12). For smaller steps the energy
loss is computed in the model under the assumption that the step length (or
relative energy loss) is small and, in consequence, the cross section can be
considered constant along the step. The loss due to the excitation is

Eexc = n1 E1 + n2 E2 (7.26)

where n1 and n2 are sampled from a Poisson distribution. The energy loss
due to ionisation can be generated from the distribution g(E) by the inverse

112
transformation method :
Z E
u = F (E) = g(x)dx (7.27)
E0

E0
E = F 1 (u) = E0
(7.28)
1 u TupTup
where u is a uniformly distributed random number [0, 1]. The contribution
coming from the ionisation will then be
n3
X E0
Eion = E0
(7.29)
j=1 1 uj TupTup
where n3 is the number of ionisations sampled from the Poisson distribution.
The total energy loss in a step will be E = Eexc + Eion and the energy
loss fluctuation comes from fluctuations in the number of collisions ni and
from the sampling of the ionisation loss.

7.2.5 Status of This Document

30.01.02 created by L. Urban
28.08.02 updated by V. Ivanchenko
17.08.04 moved to common to all charged particles (M. Maire)
04.12.04 spelling and grammar check by D.H. Wright
04.05.05 updated by L. Urban
09.12.05 updated by V. Ivanchenko
29.03.07 updated by L. Urban
11.12.08 updated by V. Ivanchenko
03.12.10 updated by L. Urban
29.11.13 updated by V. Ivanchenko

Bibliography
[1] H. Bichsel [Link]. 60 (1988) 663.
[2] K. Lassila-Perini, L. Urban [Link]. A362(1995) 416.
[3] geant3 manual Cern Program Library Long Writeup W5013 (1994)
[4] ICRU (A. Allisy et al), Stopping Powers and Ranges for Protons and
Alpha Particles, ICRU Report 49 (1993).

113
7.3 Correcting the Cross Section for Energy
Variation
As described in Sections 7.1 and 3.2.2 the step size limitation is provided
by energy loss processes in order to insure the precise calculation of the
probability of particle interaction. It is generally assumed in Monte Carlo
programs that the particle cross sections are approximately constant during a
step, hence the reaction probability p at the end of the step can be expressed
as
p = 1 exp (ns(Ei )) , (7.30)
where n is the density of atoms in the medium, s is the step length, Ei is the
energy of the incident particle at the beginning of the step, and (Ei ) is the
reaction cross section at the beginning of the step.
However, it is possible to sample the reaction probability from the exact
expression  Z Ef 
p = 1 exp n(E)ds , (7.31)
Ei
where Ef is the energy of the incident particle at the end of the step, by using
the integral approach to particle transport. This approach is available for pro-
cesses implemented via the G4V EnergyLossP rocess and G4V EmP rocess
interfaces.
The Monte Carlo method of integration is used for sampling the reaction
probability [1]. It is assumed that during the step the reaction cross section
smaller, than some value (E) < m . The mean free path for the given step
is computed using m . If the process is chosen as the process happens at the
step, the sampling of the final state is performed only with the probability p =
(Ef )/m , alternatively no interaction happen and tracking of the particle
is continued. To estimate the maximum value m for the given tracking step
at Geant4 initialisation the energy Em of absoluted maximum max of the
cross section for given material is determined and stored. If at the tracking
time particle energy E < Em , then m = (E). For higher initial energies
if E > Em then m = min((E), (E)). In the opposit case m = max .
Here is a parameter of the algorithm. Its optimal value is connected with
the value of the dRoverRange parameter (see sub-chapter 7.1), by default
= 1 R = 0.8. Note, that described method is precise if the cross section
has only one maximum, which is a typical case for electromagnetic processes.
The integral variant of step limitation is the default for the G4eIonisation,
G4eBremsstrahlung and some otehr process but is not automatically acti-
vated for others. To do so the boolean UI command can be used:

114
/process/eLoss/integral true

The integral variant of the energy loss sampling process is less dependent
on values of the production cuts [2] and allows to have less step limitation,
however it should be applied on a case-by-case basis because may require
extra CPU.

7.3.1 Status of This Document

01.12.03 integral method subsection added by V. Ivanchenko
17.08.04 moved to common to all charged particles by M. Maire
25.11.06 revised by V. Ivanchenko
09.12.10 revised by V. Ivanchenko

Bibliography
[1] [Link] et al., Proc. of Int. Conf. MC91: Detector and event
simulation in high energy physics, Amsterdam 1991, pp. 79-85. (HEP
INDEX 30 (1992) No. 3237).

[2] J. Apostolakis et al., Geometry and physics of the Geant4 toolkit for high
and medium energy applications. Rad. Phys. Chem. 78 (2009) 859.

115
7.4 Conversion from Cut in Range to Energy
Threshold
In Geant4 charged particles are tracked to the end of their range. The dif-
ferential cross section of -electron productions and bremsstrahlung grow
rapidly when secondary energy decrease. If all secondary particles will be
tracked the CPU performance of any Monte Carlo code will be pure. The
traditional solution is to use cuts. The specific of Geant4 [1] is that user
provides value of cut in term of cut in range, which is unique for defined
G4Region or for the complete geometry [2].
Range is used, rather than energy, as a more natural concept for designing
a coherent policy for different particles and materials. Definition of the cer-
tain value of the cut in range means the requirement for precision of spatial
radioactive dose deposition. This conception is more strict for a simulation
code and provides less handles for user to modify final results. At the same
time, it ensures that simulation validated in one geometry is valid also for
the other geometries.
The value of cut is defined for electrons, positrons, gamma and protons.
At the beginning of initialization of Geant4 physics the conversion is per-
formed from unique cut in range to cuts (production thresholds) in kinetic
energy for each G4MaterialCutsCouple [2]. At that moment no energy loss
or range table is created, so computation should be performed using original
formulas. For electrons and positrons ionization above 10keV a simplified
Berger-Seltzer energy loss formula (8.2) is used, in which the density correc-
tion term is omitted. The contribution of the bremsstrahlung is added using
empirical parameterized formula. For T < 10keV the linear dependence of
ionization losses on electron velocity is assumed, bremsstrahlung contribution
is neglected. The stopping range is defined as
Z T
1
R(T ) = dE. (7.32)
0 (dE/dx)

The integration has been done analytically for the low energy part and numer-
ically above an energy limit 1 keV . For each cut in range the corresponding
kinetic energy can be found out. If obtained production threshold in kinetic
energy cannot be below the parameter lowlimit (default 1 keV ) and above
highlimit (default 10 GeV ). If in specific application lower threshold is re-
quired, then the allowed energy cut needs to be extended:

G4ProductionCutsTable::GetProductionCutsTable()SetEnergyRange(lowlimit,highlimit);
or via UI commands

116
/cuts/setM inCutEnergy 100 eV
/cuts/setM axCutEnergy 100 T eV

In contrary to electrons, gammas has no range, so some approximation should

be used for range to energy conversion. An approximate empirical formula is
used to compute the absorption cross section of a photon in an element abs .
Here, the absorption cross section means the sum of the cross sections of
the gamma conversion, Compton scattering and photoelectric effect. These
processes are the destructive processes for photons: they destroy the pho-
ton or decrease its energy. The coherent or Rayleigh scattering changes the
direction of the gamma only; its cross section is not included in the absorp-
tion cross section. The AbsorptionLength Labs vector is calculated for every
material as
Labs = 5/abs . (7.33)
The factor 5 comes from the requirement that the probability of having no
destructive interaction should be small, hence

exp(Labs abs ) = exp(5) = 6.7 103 . (7.34)

The photon cross section for a material has a minimum at a certain energy
Emin . Correspondingly Labs has a maximum at E = Emin , the value of the
maximal Labs is the biggest meaningful cut in absorption length. If the cut
given by the user is bigger than this maximum, a warning is printed and the
cut in kinetic energy is set to the highlimit.
The cut for proton is introduced with Geant4 v9.3. The main goal of
this cut is to limit production of all recoil ions including protons in elastic
scattering processes. A simple linear conversion formula is used to compute
energy threshold from the value of cut in range, in particular, the cut in
range 1 mm corresponds to the production threshold 100keV .
The conversion from range to energy can be studied using G4EmCalculator
class. This class allows access or recalculation of energy loss, ranges and other
values. It can be instantiated and at any place of user code and can be used
after initialisation of Physics Lists:

G4EmCalculator calc;
[Link](range, particle, material);

here particle and material may be string names or corresponding const point-
ers to G4ParticleDefinition and G4Material.

117
7.4.1 Status of This Document
09.10.98 created by L. Urban.
27.07.01 minor revision [Link]
17.08.04 moved to common to all charged particles (mma)
04.12.04 minor re-wording by D.H. Wright
18.05.07 rewritten by V. Ivanchenko
11.12.08 minor revision by V. Ivanchenko, Geant4 v9.2
11.12.09 minor revision by V. Ivanchenko, Geant4 v9.3
09.12.10 minor revision by V. Ivanchenko, Geant4 v9.4

Bibliography
[1] Geant4 Collaboration (S. Agostinelli et al.), Nucl. Instr. Meth. A506
(2003) 250.

[2] J. Allison et al., IEEE Trans. Nucl. Sci., 53 (2006) 270.

118
7.5 Photoabsorption Ionization Model
7.5.1 Cross Section for Ionizing Collisions
The Photoabsorption Ionization (PAI) model describes the ionization energy
loss of a relativistic charged particle in matter. For such a particle, the
differential cross section di /d for ionizing collisions with energy transfer
can be expressed most generally by the following equations [1]:

2Ze4 2mv 2
 
di f ()
= ln
d mv 2 |()|2 |1 2 |
# )
1 2 ||2 F ()
arg(1 2 ) + , (7.35)
2 2

f ( )
Z
F () =
2 d ,

|( )|
0

m2 ()
f () = .
2 2 ZN ~2
Here m and e are the electron mass and charge, ~ is Plancks constant,
= v/c is the ratio of the particles velocity v to the speed of light c, Z
is the effective atomic number, N is the number of atoms (or molecules)
per unit volume, and = 1 + i2 is the complex dielectric constant of the
medium. In an isotropic non-magnetic medium the dielectric constant can
be expressed in terms of a complex index of refraction, n() = n1 + in2 ,
() = n2 (). In the energy range above the first ionization potential I1
for all cases of practical interest, and in particular for all gases, n1 1.
Therefore the imaginary part of the dielectric constant can be expressed in
terms of the photoabsorption cross section ():
N ~c
2 () = 2n1 n2 2n2 = ().

The real part of the dielectric constant is calculated in turn from the disper-
sion relation
Z
2N ~c ( )
1 () 1 = V.p. 2 2
d ,
0
where the integral of the pole expression is considered in terms of the princi-
pal value. In practice it is convenient to calculate the contribution from the
continuous part of the spectrum only. In this case the normalized photoab-
sorption cross section

119
 max 1
2 2 ~e2 Z
Z


() = () ( )d , max 100 keV
mc I1

is used, which satisfies the quantum mechanical sum rule [2]:

Z max
2 2 ~e2 Z
( )d =
.
I1 mc
The differential cross section for ionizing collisions is expressed by the pho-
toabsorption cross section in the continuous spectrum region:

2mv 2
 
di ()
= ln
d 2 |()|2 |1 2 |
# )
1 2 ||2
Z
1
( )
arg(1 2 ) + 2 d , (7.36)
2 I1 |( )|2

N ~c
2 () =
(),

Z max
2N ~c ( )
1 () 1 = V.p. 2 2
d .
I1

For practical calculations using Eq. 7.35 it is convenient to represent the

photoabsorption cross section as a polynomial in 1 as was proposed in [3]:
4
(i)
X
() = ak k ,
k=1
(i)
where the coefficients, ak result from a separate least-squares fit to experi-
mental data in each energy interval i. As a rule the interval borders are equal
to the corresponding photoabsorption edges. The dielectric constant can now
be calculated analytically with elementary functions for all , except near
the photoabsorption edges where there are breaks in the photoabsorption
cross section and the integral for the real part is not defined in the sense of
the principal value.

The third term in Eq. (7.35), which can only be integrated numerically,
results in a complex calculation of di /d. However, this term is dominant

120
for energy transfers > 10 keV , where the function |()|2 1. This is clear
from physical reasons, because the third term represents the Rutherford cross
section on atomic electrons which can be considered as quasifree for a given
energy transfer [4]. In addition, for high energy transfers, () = 1p2 / 2
1, where p is the plasma energy of the material. Therefore the factor |()|2
can be removed from under the integral and the differential cross section of
ionizing collisions can be expressed as:

2mv 2
 
di
()
= ln
d 2 |()|2 |1 2 |
# )
2 2 Z
1 || 1
arg(1 2 ) + 2 ( )d . (7.37)

2 I1

This is especially simple in gases when |()|2 1 for all > I1 [4].

7.5.2 Energy Loss Simulation

For a given track length the number of ionizing collisions is simulated by a
Poisson distribution whose mean is proportional to the total cross section of
ionizing collisions:
Z max
d( )
i = d .
I1 d
The energy transfer in each collision is simulated according to a distribution
proportional to
Z max
d( )
i (> ) = d .
d
The sum of the energy transfers is equal to the energy loss. PAI ionisation is
implemented according to the model approach (class G4PAIModel) allowing
a user to select specific models in different regions. Here is an example physics
list:
const G4RegionStore* theRegionStore = G4RegionStore::GetInstance();
G4Region* gas = theRegionStore->GetRegion("VertexDetector");
...
if (particleName == "e-")
{
...
G4eIonisation* eion = new G4eIonisation();

121
 G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");

// here 0 is the highest priority in region gas

eion->AddEmModel(0,pai,pai,gas);
...
}
...

It shows how to select the G4PAIModel to be the preferred ionisation model

for electrons in a G4Region named VertexDetector. The first argument in
AddEmModel is 0 which means highest priority.
The class G4PAIPhotonModel generates both -electrons and photons as
secondaries and can be used for more detailed descriptions of ionisation space
distribution around the particle trajectory.

7.5.3 Photoabsorption Cross Section at Low Energies

The photoabsorption cross section, (), where is the photon energy, is
used in Geant4 for the description of the photo-electric effect, X-ray trans-
portation and ionization effects in very thin absorbers. As mentioned in the
discussion of photoabsorption ionization (see section 7.5), it is convenient to
represent the cross section as a polynomial in 1 [5] :
4
(i)
X
() = ak k . (7.38)
k=1

Using cross sections from the original Sandia data tables, calculations of pri-
mary ionization and energy loss distributions produced by relativistic charged
particles in gaseous detectors show clear disagreement with experimental
data, especially for gas mixtures which include xenon.

Therefore a special investigation was performed [6] by fitting the coefficients

(i)
ak to modern data from synchrotron radiation experiments in the energy
range of 10 50 eV . The fits were performed for elements typically used
in detector gas mixtures: hydrogen, fluorine, carbon, nitrogen and oxygen.
Parameters for these elements were extracted from data on molecular gases
such as N2 , O2 , CO2 , CH4 , and CF4 [7, 8]. Parameters for the noble gases
were found using data given in the tables [9, 10].

122
7.5.4 Status of this document
01.12.05 expanded discussion by V. Grichine
08.05.02 re-written by D.H. Wright
16.11.98 created by V. Grichine
20.11.12 updated by V. Ivanchenko

Bibliography
[1] Asoskov V.S., Chechin V.A., Grichine V.M. at el, Lebedev Institute
annual report, v. 140, p. 3 (1982)

[2] Fano U., and Cooper J.W. [Link]., v. 40, p. 441 (1968)

[3] Biggs F., and Lighthill R., Preprint Sandia Laboratory, SAND 87-0070
(1990)

[4] Allison W.W.M., and Cobb J. [Link]., v.30,p.253 (1980)

[5] Biggs F., and Lighthill R., Preprint Sandia Laboratory, SAND 87-0070
(1990)

[6] Grichine V.M., Kostin A.P., Kotelnikov S.K. et al., Bulletin of the Lebe-
dev Institute no. 2-3, 34 (1994).

[7] Lee L.C. et al., J.Q.S.R.T., v. 13, p. 1023 (1973).

[8] Lee L.C. et al., Journ. of Chem. Phys., v. 67, p. 1237 (1977).

[9] G.V. Marr and J.B. West, Atom. Data Nucl. Data Tabl., v. 18, p. 497
(1976).

[10] J.B. West and J. Morton, Atom. Data Nucl. Data Tabl., v. 30, p. 253
(1980).

123
Chapter 8

Electron and Positron Incident

124
8.1 Ionization
8.1.1 Method
The G4eIonisation class provides the continuous and discrete energy losses
of electrons and positrons due to ionization in a material according to the
approach described in Section 7.1. The value of the maximum energy trans-
ferable to a free electron Tmax is given by the following relation:

E mc2 f or e+

Tmax = (8.1)
(E mc2 )/2 f or e

where mc2 is the electron mass. Above a given threshold energy the energy
loss is simulated by the explicit production of delta rays by Moller scattering
(e e ), or Bhabha scattering (e+ e ). Below the threshold the soft electrons
ejected are simulated as continuous energy loss by the incident e .

8.1.2 Continuous Energy Loss

The integration of 7.1 leads to the Berger-Seltzer formula [1]:
  
dE 2 2 1 2( + 1)
= 2re mc nel 2 ln + F (, up ) (8.2)
dx T <Tcut (I/mc2 )2

with
re classical electron radius: e2 /(40 mc2 )
mc2 mass energy of the electron
nel electron density in the material
I mean excitation energy in the material
E/mc2
2 1 (1/ 2 )
1
Tcut minimum energy cut for -ray production
c Tcut /mc2
max maximum energy transfer: for e+ , /2 for e
up min(c , max )
density effect function.
In an elemental material the electron density is
Nav
nel = Z nat = Z .
A

125
Nav is Avogadros number, is the material density, and A is the mass of a
mole. In a compound material
X X Nav wi
nel = Zi nati = Zi ,
i i
Ai

where wi is the proportion by mass of the ith element, with molar mass Ai .
The mean excitation energies I for all elements are taken from [2].
The functions F are given by :

F + (, up ) = ln( up ) (8.3)
2  2 3   2 3 4 
up 3up y up up up up
 
2 3
+ 2up up y + y
2 3 2 3 4

F (, up ) = 1 2 (8.4)
2
up
 
1  up 
+ ln [( up )up ] + + 2 + (2 + 1) ln 1
up 2
where y = 1/( + 1).
The density effect correction is calculated according to the formalism of
Sternheimer [3]:
x is a kinetic variable of the particle : x = log10 () = ln( 2 2 )/4.606,
and (x) is defined by

for x < x0 : (x) = 0

for x [x0 , x1 ] : (x) = 4.606x C + a(x1 x)m (8.5)
for x > x1 : (x) = 4.606x C

where the matter-dependent constants are calculated as follows:

p
hp = plasma energy of the medium = 4nel re3 mc2 / = 4nel re ~c
C = 1 + 2 ln(I/hp )
xa = C/4.606
a = 4.606(xa x0 )/(x1 x0 )m
m = 3.
(8.6)
For condensed media

for C 3.681 x0 = 0.2 x1 = 2
I < 100 eV
 for C > 3.681 x0 = 0.326C 1.0 x1 = 2
for C 5.215 x0 = 0.2 x1 = 3
I 100 eV
for C > 5.215 x0 = 0.326C 1.5 x1 = 3

126
and for gaseous media
for C < 10. x0 = 1.6 x1 =4
for C [10.0, 10.5[ x0 = 1.7 x1 =4
for C [10.5, 11.0[ x0 = 1.8 x1 =4
for C [11.0, 11.5[ x0 = 1.9 x1 =4
for C [11.5, 12.25[ x0 = 2. x1 =4
for C [12.25, 13.804[ x0 = 2. x1 =5
for C 13.804 x0 = 0.326C 2.5 x1 = 5.

8.1.3 Total Cross Section per Atom and Mean Free

Path
The total cross section per atom for Moller scattering (e e ) and Bhabha
scattering (e+ e ) is obtained by integrating Eq. 7.2. In Geant4 Tcut is
always 1 keV or larger. For delta ray energies much larger than the excitation
energy of the material (T I), the total cross section becomes [1] for Moller
scattering,
2re2 Z
(Z, E, Tcut ) = (8.7)
2 ( 1)
( 1)2 1
   
1 1 2 1 1 x
x + ln ,
2 2 x 1x 2 x
and for Bhabha scattering (e+ e ),
2re2 Z
(Z, E, Tcut ) = (8.8)
( 1)
   
1 1 B3 2 B4 3
1 + B1 ln x + B2 (1 x) (1 x ) + (1 x ) .
2 x 2 3
Here
= E/mc2 B1 = 2 y2
2 = 1 (1/ 2 ) B2 = (1 2y)(3 + y 2 )
x = Tcut /(E mc2 ) B3 = (1 2y)2 + (1 2y)3
y = 1/( + 1) B4 = (1 2y)3 .
The above formulas give the total cross section for scattering above the
threshold energies

thr thr
TMoller = 2Tcut and TBhabha = Tcut . (8.9)
In a given material the mean free path is then
= (nat )1 or = ( i nati i )1 .
P
(8.10)

127
8.1.4 Simulation of Delta-ray Production
Differential Cross Section
For T I the differential cross section per atom becomes [1] for Moller
scattering,
d 2re2 Z
= (8.11)
d 2 ( 1)
( 1)2 1 1 2 1
    
1 1 2 1
+ +
2 2 1 1 2
and for Bhabha scattering,
2re2 Z
 
d 1 B1 2
= + B2 B3 + B4 . (8.12)
d ( 1) 2 2

Here = T /(E mc2 ). The kinematical limits of are

Tcut 1 Tcut
0 = 2
for e e 0 = 1 for e+ e .
E mc 2 E mc2

Sampling
The delta ray energy is sampled according to methods discussed in Chapter
2. Apart from normalization, the cross section can be factorized as
d
= f ()g(). (8.13)
d
For e e scattering
1 0
f () = (8.14)
2 1 20
2
 
4 2 2 2
g() = ( 1) (2 + 2 1) + (8.15)
9 2 10 + 5 1 (1 )2
and for e+ e scattering
1 0
f () = (8.16)
2 1 0
B 0 B 1 + B 2 2 B 3 3 + B 4 4
g() = . (8.17)
B0 B1 0 + B2 20 B3 30 + B4 40
Here B0 = 2 /( 2 1) and all other quantities have been defined above.
To choose , and hence the delta ray energy,

128
 1. is sampled from f ()

2. the rejection function g() is calculated using the sampled value of

3. is accepted with probability g().

After the successful sampling of , the direction of the ejected electron is

generated with respect to the direction of the incident particle. The az-
imuthal angle is generated isotropically and the polar angle is calculated
from energy-momentum conservation. This information is used to calculate
the energy and momentum of both the scattered incident particle and the
ejected electron, and to transform them to the global coordinate system.

8.1.5 Status of this document

9.10.98 crea