Difference between revisions of "Simulations of Particle Interactions with Matter"

From New IAC Wiki
Jump to navigation Jump to search
 
(37 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
==Class Admin==
 
==Class Admin==
===[[Forest_SimPart_Syllabus]]===
 
  
The ability to simulate nature using a computer has become a useful skill for physicists who work in disciplines ranging from basic research to video games.  This class will teach you these fundamentals using a practical approach, based on a specific UNIX based programming environment, to simulate fundamental physics processes ranging from ionization and the photoelectric effect to producing anti-matter.
 
 
===Homework===
 
Homework is due at the beginning of class on the assigned day. If you have a documented excuse for your absence, then you will have 24 hours to hand in the homework after released by your doctor.
 
 
=== Class Policies===
 
http://wiki.iac.isu.edu/index.php/Forest_Class_Policies
 
 
=== Instructional Objectives ===
 
 
;Course Catalog Description
 
: Simulations of Particle Interactions with Matter 3 credits. Lecture course with monte-carlo computation requirements. Topics include: Stopping power, interactions of electrons and photons with matter, hadronic interactions, and radiation detection devices.
 
Prequisites:Math 3360. Phys 3301 or 5561.
 
;Course Description
 
: A practical course applying theoretical descriptions of fundamental particle interactions such as the photoelectric effect, compton scattering and pair production to describe multiple interactions of particles with matter using the montecarlo method. A software package known as [http://geant4.web.cern.ch/geant4/ GEANT from CERN] will be used that is freely available under the UNIX environment. The course assumes that the student has very limited experience with the UNIX environment and no experience with GEANT. Homework problems involve modifying and compiling example programs written in C++. A final project is required in which the student chooses a process to compare the predictions of GEANT with experimental data. A report is written in a format that would be publishable in a scientific journal. Publishing the report is not required but left as an option for the student.
 
 
===Objectives and Outcomes===
 
 
[[SPIM_ObjectiveNoutcomes]]
 
  
 +
[[TF_SPIM_ClassAdmin]]
  
 
== Homework Problems==
 
== Homework Problems==
Line 37: Line 18:
  
 
=Interactions of Electrons and Photons with Matter=
 
=Interactions of Electrons and Photons with Matter=
=== Bremsstrahlung===
 
;Definition:  Radiation produced when a charged particle is deflected by the electric field of nuclei in a material.
 
:Note: There is also electron-electron brehmstrahlung but the interaction is with the electric field of the materials atomic electrons.
 
  
The Cross section formula is given in Formula 3Cs, pg 928 of reference [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day8/BremXsectFormula_Rev.Mod.Phys_vol31_pg920_1959.pdf H.W. Koch & J.W Motz, Rev. Mod. Phys., vol 31 (1959) pg 920] as
+
[[TF_SPIM_e-gamma]]
  
;Note: Bethe & Heitler first calculated this radiation in 1934 which is why you will sometimes hear Bremsstrahlung radiation refererd to as Bethe-Heitler.
+
= Hadronic Interactions =
  
:<math>d \sigma = 4 Z^2r_e^2 \alpha \frac{d \nu}{\nu} \left \{ \left (1 + \left( \frac{E}{E_0} \right )^2 \right ) \left [ \frac{\phi_1(\gamma)}{4} - \frac{1}{3} \ln Z -f(Z)\right ] - \frac{2E}{3E_0} \left [  \frac{\phi_2(\gamma)}{4} - \frac{1}{3} \ln Z -f(Z)\right ] \right \} </math>
+
[[TF_SPIM_HadronicInteractions]]
  
where
+
= Final Project=
  
: <math>E_0</math> = initial total energy of the electron
+
A final project will be submitted that will be graded with the following metrics:
:<math>E</math> = final total energy of the electron
 
: <math>\nu = \frac{E_0-E}{h}</math> = energy of the emitted photon
 
:<math>Z</math> = Atomic number = number of protons in target material
 
: <math>\gamma = \frac{100 m_ec^2 h \nu}{E_0 E Z^{1/3}}</math> = charge screening parameter
 
Coulomb correction to using the Born approximation (approximation assumes the incident particle is a plane wave interacting with a static E-field the correction accounts for changes iin the plane wave due to the presence of the field)  Charge screening and the coulomb correction are different effects that have been shown to be additive/independent. [[File:Haug_2008.pdf]]
 
:<math>f(Z) = (Z \alpha)^2 \sum_1^{\infty} \frac{1}{ n [ n^2 + (Z \alpha)^2]}</math>
 
: <math>\sim (Z \alpha)^2 \left \{ \frac{1}{1+(Z \alpha)^2} +0.20206 - 0.0369(Z \alpha)^2 + 0.0083 (Z \alpha)^4 - 0.002 (Z \alpha)^6\right \}</math>
 
:<math>\alpha = \frac{1}{137}</math>
 
: <math>\phi_1</math> and <math>\phi_2</math>  = screening functions that depend on Z
 
  
if <math>Z \ge 5</math>
+
1.) The document must be less than 15 pages.
  
:<math>\phi_1(\gamma) = 20.863 - 2 \ln[1+(0.55 \gamma)^2] - 4[1-0.6e^{-0.98} - 0.4e^{-3 \gamma/2}]</math>
+
2.) The document must contain references in a bibliography (5 points) .
:<math>\phi_2(\gamma) = \phi_1(\gamma) - \frac{2}{3}(1+6.5 \gamma + 6 \gamma^2)</math>
 
  
 +
3.) A comparison must be made between GEANT4's prediction and either the prediction of someone else or an experimental result(30 points).
  
For Z<5 see [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day8/Tsai_ScreeningFunctions_Rev.Mod.Phys._vol46_pg815_1974.pdf Tsai, Rev.Mod. Phys., vol 46 (1974) pg 815]
+
4.) The graphs must be of publication quality with font sizes similar or larger than the 12 point font (10 points).
  
:if <math>3 \ge Z < 5</math> use Equation 3.46 and 3.47
+
5.) The document must be grammatically correct (5 points).
  
:if <math> Z < 2</math> use Equation 3.25 and 3.26
+
6.) The document format must contain the following sections: An abstract of 5 sentences (5 points) , an Introduction(10 points), a Theory section (20 points) , if applicable a section describing the experiment that was simulated, a section delineating the comparisons that were made, and a conclusion( 15 points).
  
;Note: Energy loss via Bethe-Bloch is due to coulomb deflection and is a '''continuous''' process while Bremsstrahlung is a '''discrete''' process (emission of photons)
+
=Resources=
  
;We now know 2 ways charged particles can loose energy when passsing through matter.
+
[http://geant4.web.cern.ch/geant4/  GEANT4 Home Page]
  
;Energy loss: <math>\left ( \frac{dE}{dx} \right )_{tot} = \left ( \frac{dE}{dx} \right )_{rad} + \left ( \frac{dE}{dx} \right )_{col}</math>
+
[http://root.cern.ch ROOT Home page]
:<math>{rad} \Rightarrow</math> : Bremsstrahlung
 
:<math>{col} \Rightarrow</math> : Bethe-Bloch (collision)
 
  
;<math>-\left ( \frac{dE}{dx} \right )_{rad} = N \int_o^{\nu_0} \left ( h \nu \right ) \left ( \frac{d \sigma}{d \nu} \right ) d \nu</math>
+
[http://conferences.fnal.gov/g4tutorial/g4cd/Documentation/WorkshopExercises/  Fermi Lab Example]
  
where
 
  
: <math>N = \frac{ number\; atoms}{cm^3} = \frac{\rho N_a}{A} = \frac{density \;\times\; Avagadros\;\;Number}{Atomic number}</math>
+
[http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html NIST Range Tables]
:<math>\left ( h \nu \right )</math> = Energy of emmitted photon
 
:<math>\left ( \frac{d \sigma}{d \nu} \right )</math> = Probabitlity of Energy loss
 
  
 +
[http://ie.lbl.gov/xray/  X-ray specturm]
  
The quantity<math> \Phi_{rad} </math> is defined such that
+
[[Installing_GEANT4.9.3_Fsim]]
  
:<math>\Phi_{rad} \equiv \frac{1}{E_0} \int_o^{\nu_0} \left ( h \nu \right ) \left ( \frac{d \sigma (E_0,\nu)}{d \nu} \right ) d \nu</math>
+
== Saving/restoring Random number seed==
  
<math>\Phi_{rad}</math> is a macroscopic function of a given material rather than just the energy <math>\nu</math> which we will use to define a common property of materials known as the radiation length <math>\left ( R_0=\frac{1}{N\Phi_{rad}} \right )</math>
+
You save the current state of the random number generator with the command
  
:<math>\Phi_{rad} = 4Z^2r_e^2 \alpha \left \{ { \left [  \ln(\frac{2E}{m_ec^2}) - \frac{1}{3} - f(z) \right ] \gamma >> 1 \atop  \left [  \ln(183E^{-1/3}) + \frac{1}{18} - f(z) \right ] \gamma \sim 0} \right .</math>
+
/random/setSavingFlag 1
  
where
+
/run/beamOn 100
  
: <math>\gamma >> 1</math> case is no screening and <math>1 \ll \frac{E_0}{m_e c^2}  < \frac{1}{\alpha Z^{1/3}}</math>
+
/random/saveThisRun
: <math>\gamma \; \sim \; 0</math> case has <math> \frac{E_0}{m_e c^2}  \gg \frac{1}{\alpha Z^{1/3}}</math>
 
  
The energy loss equation becomes
+
A file is created called
  
: <math>- \frac{dE}{dx} = N E_0 \Phi_{rad}</math>
+
currentEvent.rndm
  
;Note: for intermediate value of <math>\gamma</math> you need to integrate numerically
+
/control/shell mv currentEvent.rndm currentEvent10.rndm
:<math>\left ( \frac{dE}{dx} \right )_{rad}\propto Z^2E</math> : Bremsstrahlung
 
:<math>\left ( \frac{dE}{dx} \right )_{col}\propto Z \ln (E)</math> : Bethe-Bloch
 
  
The illustration below shows the relative contributions of Bethe-Bloch and Bremsstrahlung to the total energy loss according to the above functional dependence.  At low energies the physics of collisions dominates the loss (Bethe-Bloch) and as energy increases the discrete loss by radiation begins to dominate. 
 
  
[[Image:SPIM_Bethe-Brem_Eloss-vs-Energy.jpg | 600 px]]
+
You can restore the random number generator and begin generating random number from the last save time
  
 +
/random/resetEngineFrom currentEvent.rndm
  
;Critical Energy<math> E_C</math>:
+
==Building GEANT4.11==
At the critical energy <math>E_C</math> the two energy loss processes contribute equally to the total energy lost by a charged particle interacting with matter.
 
:<math>E_C \equiv</math> energy at which <math>\left ( \frac{dE}{dx} \right )_{rad} = \left ( \frac{dE}{dx} \right )_{col}</math>
 
  
In the PDG
+
===4.11.2===
 +
[[TF_GEANT4.11]]
  
:<math>E_C \sim \frac{800 MeV}{Z+1.2}</math>
+
==Building GEANT4.10==
  
;Examples:
 
{| border="3"  cellpadding="20" cellspacing="0"
 
|colspan="2"|
 
====Critical Energy E_C====
 
|-
 
|Material
 
|<math>E_C</math> (MeV)
 
|-
 
|Pb
 
|9.51
 
|-
 
|Fe
 
|27.4
 
|-
 
|Cu
 
|24.8
 
|-
 
|Al
 
|51
 
|}
 
  
==== Electron-Electron Bremsstrahlung ====
+
===4.10.02===
; Electron electron bremsstrahlung: The radiation produced as 2 electrons pass near eachother
 
: <math>d \sigma</math> is essentially the same except you have <math>z=1</math> thereby adding a <math>Z</math> term and not a <math>Z^2</math> term
 
  
reference:pg 947 from Koch and Motz, Rev. Mod. Phys, vol 31 (1959) pg 920 [[Image:SPIM_Koch andMotz_RevModPhysv31_1959pg920.pdf]]
+
[[TF_GEANT4.10.2]]
 +
===4.10.01===
  
as a result
+
[[TF_GEANT4.10.1]]
  
<math>d \sigma_{tot} = \frac{Z(Z+1)}{Z^2} d \sigma_{Brem}</math>
+
==Building GEANT4.9.6==
: = <math>4 Z(Z+1)r_e^2 \alpha \frac{d \nu}{\nu} \left \{ \left (1 + \left( \frac{E}{E_0} \right )^2 \right ) \left [ \frac{\phi_1(\gamma)}{4} - \frac{1}{3} \ln Z -f(Z)\right ] - \frac{2E}{3E_0} \left [  \frac{\phi_2(\gamma)}{4} - \frac{1}{3} \ln Z -f(Z)\right ] \right \} </math>
 
  
Most calculations ignore electron-electron Brehmsstrahung because its linear in Z and doesn;t become important until low Z where measured atomic form factors are actually used and not Form factors calulated by the Thomas-Fermi-Moliere Model (Z>4).
+
[[TF_GEANT4.9.6]]
  
 +
==Building GEANT4.9.5==
  
==== Coherent Bremsstrahlung ====
+
[[TF_GEANT4.9.5]]
 
 
 
 
The above equations represent the physics of incoherent bremsstrahlung production.
 
 
 
 
 
http://137.99.79.133/halld/tagger/references/inelasticBrems-94.pdf
 
 
 
==== Radiation Length (Xo)====
 
 
 
;Radiation Length<math> (X_0)</math>: The distance an electron travels through matter until loosing <math>\frac{1}{e}</math> of its energy due to radiation <math>\left ( \frac{dE}{dx} \right )_{rad}</math>.
 
 
 
in the high energy limit where <math>\left ( \frac{dE}{dx} \right )_{col}</math> can be ignored <math>( E > E_C )</math>
 
 
 
:<math>\frac{dE}{dx} = NE_0 \Phi_{rad} \Rightarrow \int_{E_0}^{E} \frac{dE}{E} = -\int_0^X N \Phi_{rad} dx</math>
 
: <math>\Rightarrow \ln(\frac{E}{E_0}) =-N \Phi_{rad}X</math>
 
 
 
or
 
 
 
: <math>E = E_0 e^{-N \Phi_{rad}X} = E_0 e^{-\frac{X}{X_0}}</math>
 
 
 
where
 
 
 
: <math>X_0 \equiv \frac{1}{N \Phi_{rad}}</math> = Radiation Length of a given material
 
 
 
ie:
 
:if <math> X = X_0</math>  Then <math>E=\frac{1}{e} E_0</math> = Energy of electron after traveling a distance of <math>X_0</math> through the material
 
 
 
  
{| border="3"  cellpadding="20" cellspacing="0"
+
An old version of Installation notes for versions prior to 9.5
|colspan="2"|
 
  
====Table of Radiation Lengths for several materials====
+
[http://brems.iac.isu.edu/~tforest/NucSim/Day3/ Old Install Notes]
|-
 
|-
 
|Material
 
|<math>X_0</math> (cm)
 
|-
 
|Air
 
|30,050
 
|-
 
|Al
 
|8.9
 
|-
 
|Cu
 
|1.43
 
|-
 
|Fe
 
|1.76
 
|-
 
|H2O
 
|36.1
 
|-
 
|NaI
 
|2.59
 
|-
 
|Pb
 
|0.56
 
|-
 
|Polystyrene
 
|42.9
 
|-
 
|Scintillators
 
|42.2
 
|-
 
|}
 
  
  
;If we have complete screening <math> (\gamma=0)</math>
+
Visualization Libraries:
:Then <math>\frac{1}{X_0} = N \Phi_{rad} = 4 \alpha r_e^2 \frac{N_A}{A} \left \{ Z^2 \left [ L_{rad} - f(Z)\right ] + ZL_{rad}^{\prime}\right \}</math>
 
: = <math>\frac{Z^2 \left [ L_{rad} - f(Z)\right ] + ZL_{rad}^{\prime}}{716.408 \frac{g}{cm^2}A}</math>
 
  
where
+
[http://www.opengl.org/ OpenGL]
  
: <math>L_{rad} \equiv  \frac{1}{4} [\phi_1(\gamma=0) - \frac{4}{3} \ln(Z)] = \left \{ {1 + \int_0^{m_e} [ 1-\frac{F(q)}{Z}]^2 \frac{dq}{q} \;\;\; Z\le4\atop \ln(184.15 Z^{-1/3} \;\;\; Z>4} \right . </math>= radiation logarithm for elastic Atomic scattering
+
[http://geant4.kek.jp/~tanaka/DAWN/About_DAWN.html  DAWN]
:<math>L_{rad}^{\prime} \equiv \frac{1}{4} [ \phi_2(\gamma=0) - \frac{8}{3}\ln Z]= \left \{ {1 + \frac{1}{2}\int_0^{m_e} \frac{G_2^{inel}(t)}{Z}\frac{dt}{t} \;\;\; Z\le4\atop \ln(1194 Z^{-2/3} \;\;\; Z>4} \right .</math> = radiation logarithm for inelastic Atomic scattering
 
:<math>f(Z) = \alpha^2 Z^2 \left [ \frac{1}{1+\alpha^2 Z^2} + 0.20206 - 0.0369 \alpha^2 Z^2 + 0.0083 \alpha^4 Z^4 -0.002 \alpha^6 Z^6 \right ]</math>  :Z < 92
 
  
  
;Quick <math>X_0</math> Estimates
+
[http://doc.coin3d.org/Coin/ Coin3D]
: <math> X_0 = \frac{716.4 \left ( \frac{g}{cm^2} \right ) A }{Z(Z+1) \ln \left ( \frac{287}{\sqrt{Z}} \right )}</math>
 
  
;Examples of Radiation length
+
==Compiling G4 with ROOT==
  
;<math>\frac{1}{e} = \frac{1}{2.72} \sim \frac{1}{3}</math> : an electron has lost 1/3 of its original energy after traveling 1 radiation length (1 <math>X_0</math>) through the material
+
These instruction describe how you can create a tree within ExN02SteppingVerbose to store tracking info in an array (max number of steps in a track is set to 100 for the desired particle)
  
;<math>\frac{1}{e^2}  \sim \frac{1}{7}</math> : an electron has lost 1/7 of its original energy after traveling 2 radiation lengths (2 <math>X_0</math>) through the material
+
[[G4CompileWRootforTracks]]
  
;<math>\frac{1}{e^3}  \sim \frac{1}{20}</math> : an electron has lost 1/20 of its original energy after traveling 3 radiation lengths (3 <math>X_0</math>) through the material
+
==Using SLURM==
  
; After 2.3 radiation lengths the electron energy is down by a factor of 10 from its original value.
+
http://slurm.schedmd.com/quickstart.html
  
====Bremsstrahlung in GEANT 4 ====
 
  
GEANT4 uses an energy cut off  <math>(T_c, k_c)</math> to decide whether to use a continuous energy loss algorithm (msc, Bethe-Bloch, soft photons) or to generate a secondary particle (photon) and use Bremsstrahliung.
+
https://rc.fas.harvard.edu/resources/documentation/convenient-slurm-commands/
  
:<math>T_C</math> = incident particle K.E. cutof = secondary particle production threshold
+
===simple batch script for one process job===
:<math>k_c</math> = photon energy cutoff below which photons are treated as continuous energy loss.
 
  
;if <math>E_{secondary}<T_C</math> then no photon is created and the effect of the soft photon reaction is treated as a continuous energy loss via
+
create the file submit.sbatch below
:<math>E_{loss}^{Brem} (Z,T,k) = \int_0^{k_c} \left ( \frac{ d \sigma(z,T,k)}{dk}\right ) k dk</math> = continuous energy loss via "soft" photon emission
 
  
:<math>\frac{ d \sigma(Z,T,k)}{dk}</math> = cross sections parametrerized by the Evaluated Electrons Data Library (EEDL)
+
<pre>
: reference: J. Tuli, "Evaluated Nuclear Structure Data File", BNL-NCS - 51655 -Rev 87, 1987 from Brookhaven Nat. Lab
+
#!/bin/sh
: see  [http://www.nndc.bnl.gov National Nuclear Data Center]:
+
#SBATCH --time=1
 +
cd src/PI
 +
./PI_MC 100000000000000
 +
</pre>
  
; Note: Soft photons are photons created in the scattering process which have less energy than the energy of the particles participating in the interaction.  Soft photon are not energtic enough to be detected.
+
the execute
  
To improve simulation speed though, GEANT 4 actually uses a fit to the above cross sections such that
+
:sbatch submit.sbatch
  
:<math>E_{loss}^{Brem} (Z,T,k) = (2 - C_{th}Z^{1/4}) \frac{Z(Z+\epsilon_{\ell})(T+m)^2}{T+2m} \left [ \frac{k_C}{T}\right ]^{\beta} \frac{a+b\frac{T}{T_{\ell m}}}{1+c \frac{T}{T_{\ell m}}} \frac{f_{\ell}}{N_a}</math>
+
check if its running with
  
where
+
:squeue
  
: <math>m</math> = electron mass
+
to kill a batch job
:<math>T</math> = kinetic energy of incident particle
 
:<math>N_a</math> = Avagadros number
 
:<math>\epsilon , \beta, C_{th}, a, b, c</math> = constants
 
:<math>f_{\ell}</math> = polynomial (in log(T) ) chosen to fit the data
 
  
;if <math>E_{secondary}>T_C</math> then a photon is created and tracked
+
:scancel JOBID
  
: The energy of the emitted photon is determine by sampling a probability distribution from
+
===On minerve===
:S. Seltzer and M. Berger, Atomic Data & Nucle. Data tables, vol 35 (1986) pg 345-418
 
  
and
+
Sample script to submit 10 batch jobs.
  
: the angular distribution (<math>\cos(\theta)</math>)  is sampled according to
+
the filename is minervesubmit and you run like
:E. Acosta, Appl. Phys. Letter, vol 80 # 17 (2002) pg 3228-3230
+
source minervesubmit
  
 +
<pre>
 +
cd /home/foretony/src/GEANT4/geant4.9.5/Simulations/N02wROOT/batch
 +
qsub submit10mil
 +
qsub submit20mil
 +
qsub submit30mil
 +
qsub submit40mil
 +
qsub submit50mil
 +
qsub submit60mil
 +
qsub submit70mil
 +
qsub submit80mil
 +
qsub submit90mil
 +
qsub submit100mil
 +
</pre>
  
Note
+
The file submit10mil looks like this
:* The MC program PENELOPE was used to generate the energy distributions that are sampled
+
<pre>
:* GEANT4 uses a modified version of  base equations for <math>e^+e^-</math> bremsstrahlung with model corrections for <math>e^+</math>
+
#!/bin/sh
 +
#PBS -l nodes=1
 +
#PBS -A FIAC
 +
#PBS -M foretony@isu.edu
 +
#PBS -m abe
 +
#
 +
source /home/foretony/src/GEANT4/geant4.9.5/geant4.9.6-install/bin/geant4.sh
 +
cd /home/foretony/src/GEANT4/geant4.9.5/Simulations/N02wROOT/batch/10mil
 +
../../exampleN02 run1.mac > /dev/null
 +
</pre>
  
; LPM effect:  There is also a correction kown as the Landau Pomeranchuk Migdal (LPM) effect which corrects for multiple scattering experiences by the electron during the scattering which causes the emission of a photon.
 
  
====Bragg's Rule for compound materials ====
+
use
  
The radiation length for compounds and mixtures is determined by parallel weighting (resistors in parallel)
+
qstat
  
:<math>\frac{1}{L_{rad}} = \omega_1 \left ( \frac{1}{L_{rad}} \right )_1 + \omega_2 \left ( \frac{1}{L_{rad}} \right )_2</math>
+
to check that the process is still running
  
where
+
use
  
: <math>\omega_1</math> = fraction , by wieght, of each element in the mixture/compound.
+
qdel jobID#
: <math>= \frac{a_i A_i}{A_m}</math>
 
: <math>a_1</math> = # of atoms of element "i"
 
:<math>A_i</math> = atomic # of element "i"
 
:<math>A_m = \sum a_i A_i</math> = effective atomic mass of the compound/mixture
 
  
=== Photo-electric effect===
+
if you want to kill the batch job, the jobID number shows up when you do stat.
  
The photo-electric effect identifies the physics process by which bound electrons in an atom are liberated by an interaction with an incident photon.
+
for example
 +
<pre>
 +
[foretony@minerve HW10]$ qstat
 +
Job id                    Name            User            Time Use S Queue
 +
------------------------- ---------------- --------------- -------- - -----
 +
27033.minerve            submit          foretony        00:41:55 R default       
 +
[foretony@minerve HW10]$ qdel 27033
 +
[foretony@minerve HW10]$ qstat
 +
</pre>
  
[[Image:SPIM_PhotoElectricEffectSchematic.jpg | 400 px]]
+
==Definitions of Materials==
  
:<math>E_f = E-E_B = h \nu - E_B</math>
+
[[File:MCNP_Compendium_of_Material_Composition.pdf]]
  
where
+
==Minerve2 GEANT 4.10.1 Xterm error==
  
:<math>h\nu</math> = incident photon energy
 
:<math>E_B</math> = electron binding energy
 
  
==== Moseley's Law ====
+
On OS X El Capitan V 10.11.4 using XQuartz
  
Moseley's law approximates the binding energies of electrons in atoms as
+
~/src/GEANT4/geant4.10.1/Simulations/B2/B2a/exsmpleB2a
  
: <math>E_B = 13.605 \frac{\left ( Z-k_s \right )^2 }{n^2}</math> (eV)
+
<pre>
  
X-ray electron shells are labeled K,L,M
+
# Use this open statement to create an OpenGL view:
 +
/vis/open OGL 600x600-0+0
 +
/vis/sceneHandler/create OGL
 +
/vis/viewer/create ! ! 600x600-0+0
 +
libGL error: failed to load driver: swrast
 +
X Error of failed request:  BadValue (integer parameter out of range for operation)
 +
  Major opcode of failed request:  150 (GLX)
 +
  Minor opcode of failed request:  3 (X_GLXCreateContext)
 +
  Value in failed request:  0x0
 +
  Serial number of failed request:  25
 +
  Current serial number in output stream:  26
 +
</pre>
  
  
{| border="3"  cellpadding="20" cellspacing="0"
 
| Shell || n || Spect. Notation (low E)|| Spect. Notation (High E) || k_s
 
|-
 
|K || 1 || <math>1S_0</math> ||  || 3
 
|-
 
|L || 2 || <math>2P_{3/2} , 2P_{1/2}</math> ||  3P{3/2}|| 5
 
|-
 
|}
 
 
;Example: Binding Energies for Argon (A=18) "Chemical Rubber Company Handbook of Chemistry and Physics", CRC press. Boca Raton FL, 81st ed, 2000.
 
 
{| border="3"  cellpadding="20" cellspacing="0"
 
| Shell || n || Spect. Notation || colspan="3"| Binding Energy (eV)
 
|-
 
| || || || Measured ||GEANT4 || Moseley
 
|-
 
|<math>K</math> || 1 || <math>1S</math> || 3218  || 3178 || 3061
 
|-
 
|<math>L_{I}</math> || 2 || <math>2S</math> ||  328 ||313.5 || 575
 
|-
 
|<math>L_{II}</math> || 2 || <math>2P_{1/2}</math> ||  251 || 247  || 575
 
|-
 
|<math>L_{III}</math> || 2 || <math>2P_{3/2}</math> ||  248.4|| 247 || 575
 
|-
 
|}
 
 
 
Binding energies for a few common elements
 
 
{| border="3"  cellpadding="20" cellspacing="0"
 
| Element || colspan="5"| Binding Energy (eV)
 
|-
 
|  || n=1 || colspan="2"| n=2|| colspan="2"| n=3
 
|-
 
|B || 201 || 14.2 || 8.3  ||  ||
 
|-
 
|C || 298 || 17.9 || 11.4  ||  ||
 
|-
 
|N || 450 || 26 || 15  ||  ||
 
|-
 
|O || 548 || 33 || 13.6  ||  ||
 
|-
 
|Na || 1083 || 71|| 38  || 5.2 ||
 
|-
 
|Mg || 1313 || 94 || 55  || 7.7 ||
 
|-
 
|Al || 1573 || 126 || 81  || 11 || 6
 
|-
 
|Si || 1854|| 157 || 107  || 15 || 8.2
 
|-
 
|P || 2167 || 195 || 141  || 20 || 10.5
 
|-
 
|S || 2490 || 236 || 172  || 21.3 || 10.4
 
|-
 
|Cl || 2844 || 279 || 210  || 25|| 13
 
|-
 
|K || 3615 || 386 || 303  || 41 || 25
 
|-
 
|}
 
 
==== Photo-electric cross section====
 
 
;the most general expression
 
 
:<math>\frac{d \sigma}{d \Omega} = 32 \frac{e^2}{4 \pi \hbar c} | \vec{k}^{\prime} | \frac{\hbar}{mc} \frac{c}{\omega}  \left ( \frac{z}{a_o} \right )^5 \frac{ \left ( \hat{\epsilon} \cdot \vec{k}^{\prime} \right ) ^2}{ \left [ \frac{z^2}{a_0^2} + q^2 \right ]^4 } </math>
 
 
where
 
 
: <math>\vec{k}^{\prime} = \frac{\vec{p}}{h}</math> = scattered electron wave number [<math>\frac{1}{m}</math> ]
 
: <math>\omega = c | \vec{k} |</math> = incident photon wave frequency
 
: <math>a_o = \frac{\hbar^2}{(mc)^2}</math>
 
: <math>\hat{\epsilon}</math> = incident photon polarization
 
:<math>\vec{q} = \vec{k} - \vec{k}^{\prime}</math> = momentum given to the atom divided by Plank's constant (h)
 
 
if the electron's K.E. after emission is larger than its binding energy
 
 
then
 
 
: <math>k^{\prime} \approx \left ( \frac{2m \omega}{\hbar} \right )^{1/2}</math>
 
:<math>a_0 = \frac{r_0}{\alpha^2}</math>
 
: <math>\hat{\epsilon} \cdot {\vec{k}}^{\prime} = k^{\prime} \sin(\theta) \cos(\phi)</math>
 
 
<math>\Rightarrow</math>
 
 
:<math>\frac{d \sigma}{d \Omega} = \alpha^4 r_0^2 Z^5 \left ( \frac{mc^2}{\hbar \omega}\right)^{7/2} \frac{4 \sqrt{2} \sin^2(\theta) \cos^2(\phi)}{\left [ 1 - \frac{v}{c}\cos(\theta)\right ]^4} </math>
 
 
For K shell emmission
 
 
: <math>\sigma = 4 \sqrt{2} \frac{8 \pi r_o^2}{3} \alpha^4 Z^5 \left( \frac{mc^2}{\hbar \omega}\right)^{7/2}</math>
 
 
at higher energies <math>(\hbar \omega \gg mc^2 )</math> (the ultra-relativistic limit
 
 
: <math>\sigma = \frac{3}{2} \frac{8 \pi r_o^2}{3} \alpha^4 Z^5 \left( \frac{mc^2}{\hbar \omega}\right)</math>
 
 
References
 
 
http://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html
 
 
http://www-amdis.iaea.org/LANL/
 
 
http://www.nist.gov/pml/data/ionization/index.cfm
 
 
==== Mass Attenuation Coefficient====
 
 
The mass attenuation coefficient <math>(\frac{\mu}{\rho})</math> is used to describe the attenutation of a photon interacting with matter via the photo-electric effect ( absorption), coherent scattering (Rayleigh), incoherent scattering (Compton), or pair production.
 
 
:<math>\mu</math> = linear photo-electric attenutaion
 
:<math>\rho</math> = density of the material
 
: <math>\frac{\mu}{\rho} = \frac{\sigma_{PE}}{amu}</math>
 
 
 
;Eample: Below is an example of the mass atenuation coefficeint as a function of the incident photon energy
 
 
[[Image:SPIM_MassAttenCoef_H2O.jpg | 400 px]]
 
 
a 10 keV photon (0.01 MeV) will have <math>\frac{\mu}{\rho} = 1 \frac{cm^2}{g}</math> when traveling through water <math>(\rho = 1 \frac{g}{cm^3})</math>
 
 
:<math>\Rightarrow \mu = 1 \frac{cm^2}{g} \rho = 1 \frac{1}{cm}</math> = attenuation coefficient
 
 
:<math> I = I_0 e^{-\mu x}</math> = intensity of light
 
:if <math>I= \frac{I_0}{2}</math>
 
:<math>\Rightarrow x = \lambda_{1/2}</math> = half length = <math>- \frac{1}{\mu} \ln (\frac{1}{2})</math> = 0.69 cm
 
 
This means that 1/2 of the photons impinging on water get absorbed by the water atoms after a depth of 0.69 cm.
 
 
;Scaling
 
: Sometimes when <math>\frac{\mu}{\rho}</math> is not available for your material you can scale a <math>\frac{\mu}{\rho}</math> of a material with similar atomic number using the equation
 
 
: <math>\left ( \frac{\mu}{\rho} \right)_{unknown} = \left ( \frac{\mu}{\rho} \right)_{known} \left [ \frac{Z^n}{A \rho}\right]</math>
 
 
where the coefficient<math> n</math> varies with the photon energy from 4 <math>\rightarrow</math> 5 according to:
 
 
[[Image:SPIM_ScalingMassAttCoeff.jpg | 400 px]]
 
 
==== GEANT4 ====
 
GEANT4 uses a parameterization of photon absorption cross sections to determine the mean free path, atomic shell data to determine the ejected electron energy, and the k-shell angular distribution to determine the direction of the ejected electron.
 
 
<u>The fit to the photoabsorption cross sections</u>
 
 
: The photoabsorption cross section is parametrized according to
 
:<math>\sigma(Z,E_{\gamma}) = \frac{a(Z,E_{\gamma})}{E_{\gamma}} +  \frac{b(Z,E_{\gamma})}{E^2_{\gamma}} + \frac{c(Z,E_{\gamma})}{E^3_{\gamma}} + \frac{d(Z,E_{\gamma})}{E^4_{\gamma}}</math>
 
 
where
 
 
:<math>a(Z,E_{\gamma}) ,\; b(Z,E_{\gamma}) ,\; c(Z,E_{\gamma}), \;d(Z,E_{\gamma})</math> are determined by a least squares fit to the data as outlined in
 
 
F. Biggs & R. Lighthill, Sandia Lab Preprint, SAND 87-0070 (1990)
 
 
[[File:Biggs_Lighthill_SandiaPreprint_87-0070.pdf]]
 
 
You select <math>E_{\gamma}</math> by sampling from a distribution generated by the above cross section.
 
 
<u>The mean free path (<math>\lambda</math>) </u>of the photon through the material is given by
 
 
:<math>\lambda(E_{\gamma}) = \frac{1}{\sum_i N_i \sigma(Z_i,E_{\gamma})}</math>
 
 
where
 
 
:<math>N_i = \frac{\mbox{number of Atoms of} \quad i^{th} \quad\mbox{element } \quad Z_i \quad \mbox{in the material}}{\mbox{Volume}}</math>
 
 
<u> K.E. of ejected electron </u>
 
 
Given that a photo-electric event happens then the energy of the ejected electron is given by
 
 
:<math>T_{p.e.} = E_{\gamma}-B_{shell}(Z_i)</math>
 
 
where
 
 
:<math>E_{\gamma}</math> = energy of the incident photon
 
:<math>B_{shell}(Z_i)</math> = electron shell energy from the closest available atomic shell as tabulated in data/G4EMLOW/fluor/binding.dat
 
 
The shell is selected according to the shell cross sections
 
 
 
<u> Electron direction</u>
 
 
The ejected electron is chosen by an  angle according to the Souter-Gavrila distribution ([http://www.physics.isu.edu/~tforest/Classes/NucSim/Day9/PhotoElec_X-sect_Gavrila_Phys.Rev_vol113_pg514_1959.pdf Gavrila_M._Phys.Rev._vol113_1959_pg514])  in the  "standard" package such that
 
 
:<math>\cos(\theta) = \frac{rnd + \beta}{rnd\times \beta +1}</math>
 
 
where rnd is a random number chosen such that
 
 
:<math>\frac{1-\cos^2(\theta)}{(1-\beta \cos(\theta))^2} \times \left [ 1 + b(1-\cos(\theta) \right ] < rnd \times  \left \{ { \gamma^2(1+b(1-\beta)) \quad \gamma < 2 \atop \gamma^2(1+b(1+\beta) ) \quad \mbox{otherwise} } \right .</math>
 
:<math>b = \frac{\gamma (\gamma -1) ( \gamma -2)}{2}</math>
 
:<math>\gamma = 1+ \frac{\mbox{K.E.}_{e^-}}{m_ec^2}</math>
 
 
===== Physics Models=====
 
======G4PhotoElectricEffect======
 
 
This model will generate an ionized electron which is about the same as the incident photon energy.  Don't use this one if your simulation is sensitive to atomic energy levels (ie; looking at keV energy effects).  This should be O.K. if you are just interested in attenuating photons.
 
 
======G4LowEnergyPhotoElectric======
 
 
This process will generate ionized electrons for each possible electron binding energy which is less than the incident photon energy.  It should be cross section weighted.
 
 
This model seems to break when <math>E_{\gamma} \le</math> 100 keV  (GEANT4 version 4.8)
 
 
======PAI Model ======
 
 
PhotoAbsorption Ionizaton (PAI) Model
 
 
The PAI model uses a least squares fit  of a 4th order polynomial in <math>\frac{1}{\omega}</math> to the experimental photoabsorption data for the cross section such than
 
 
:<math>\sigma(\omega) = \sum_i^4 \frac{a_k(E)}{\omega^k}</math>
 
 
where
 
:<math>a_k(E)</math> = fit coefficent for energy bin <math>E</math>
 
:<math>\omega</math> = energy transfered in the ionization collision
 
 
=== Compton Scattering ===
 
 
Compton scattering is like the photo-electric effect except the photon isn't absorbed but scattered by atomic electrons. 
 
 
"Ideal" compton scattering is described in terms of free electrons.
 
 
[[Image:SPIM_IdealComptonScattering.jpg | 400 px]]
 
 
The collision is elastic
 
 
:<math>\lambda^{\prime} = \lambda + \lambda_C (1-\cos(\theta)) = \frac{2 \pi}{\omega^{\prime}} = \frac{ch}{E_{\gamma}^{\prime}} = \frac{12,400 \mbox{Angstroms}}{E_{\gamma}^{\prime}}</math>
 
 
: <math>\lambda_C</math> = electron compton wavelength = <math>\frac{h}{m_ec} = 2.43 \times 10^{-12} m</math>
 
: <math>E_k = \hbar \omega \frac{\lambda_C}{\lambda} \frac{1- \cos(\theta)}{1 + \frac{\lambda_C}{\lambda} \left (1 - \cos(\theta) \right )}</math> = electron final kinetic energy
 
: <math>\phi = \cot \left [ \left ( 1+ \frac{\lambda_C}{\lambda}\right ) tan(\frac{\theta}{2}\right ]</math> = ejected electron angle w.r.t  original photon direction
 
 
;Note: <math>\phi_{max} = \frac{\pi}{2}</math> : No electrons can be backscattered in the compton process.
 
 
; The photon can backscattered:<math>\theta = \pi</math> = Max energy transfered to the <math>e^-</math>
 
:<math>E_k(max) = \frac{2 \hbar \omega \lambda_C}{\lambda + 2 \lambda_C}</math> = The max energy transfer point corresponds to the "compton" edge
 
 
;Example: Find <math> E_k(max)</math> of the compton edge for a given <math>E_{\gamma}</math>.
 
 
:<math>E_k(max) = \frac{\hbar \omega }{\frac{\lambda}{2 \lambda_C} + 1} = \frac{E_{\gamma} }{\frac{\lambda}{2 \lambda_C} + 1}</math>
 
 
 
: <math>\lambda = \frac{c}{\nu} = \frac{c h}{E_{\gamma}}</math>
 
 
: <math>E_k(max) =\frac{E_{\gamma} }{\frac{c h}{2 \lambda_C E_{\gamma}} + 1} =\frac{E_{\gamma} }{\frac{3 \times 10^8 \frac{m}{s} 4.14 \times 10^{-15} \frac{eV}{s}}{2 \times 2.43 \times 10^{-12} m E_{\gamma}} + 1}</math>
 
: <math>\approx 4 \times 10^{-6} \frac{E_{\gamma}^2}{1 + 4 \times 10^{-6} E_{\gamma}}</math>
 
 
: If <math>E_{\gamma}</math> = 8 keV
 
:Then <math>E_k(max) \approx 256</math> eV = max energy lost by photon and given to electron
 
 
 
 
[[Image:SPIM_EnergyDistributionComptonElectrons.jpg | 400 px]]
 
 
==== Cross Section ====
 
 
The Klein-Nishina formula (Oskar Klein & Yoshio Nashina, Z. fur Phys., vol 52 (1929), pg 853 ) is given as
 
 
:<math>\frac{d \sigma}{d \Omega} = \frac{r_e^2}{2} \frac{1 + \cos^2(\theta) + \frac{\xi^2 \left [ 1+ \cos(\theta) \right ]^2}{1 + \xi \left( 1+ \cos(\theta) \right)}}{\left[ 1+ \xi (1-\cos(\theta) ) \right ]^2}</math>
 
 
where
 
 
: <math>\xi = \frac{h \nu}{m_e c^2} = \frac{E_{\gamma}}{E_0^{e^-}} = \frac{E_{\gamma}}{0.511 MeV} \approx 2\frac{E_{\gamma}}{MeV}</math>
 
 
;Note: The above cross section is for a free electron.  Multiple by <math>Z</math> (the number of electrons in the target) to get the atomic cross section.
 
 
After integrating over <math>d \Omega</math>
 
 
:<math>\sigma_{compt} = 2 \pi r_e^2 \left \{ \frac{1+\xi}{\xi^2} \left [ \frac{2(1+\xi)}{1 + 2\xi} - \frac{1}{\xi} \ln(1+2 \xi) \right ]\right \}</math>
 
 
[[Image:SPIM_ComptonScatt_KleinNishiwaXsect.jpg | 400 px]]
 
 
==== Energy Distribution ====
 
 
The compton electron energy distribution can be evaluated from the differential cross section below
 
 
:<math>\frac{d \sigma}{d E^{e^-}} = \frac{\pi r_e^2}{m_e c^2 \xi^2} \left [ 2 + \frac{s^2}{\xi^2 (1-s)^2} + \frac{s}{1-s} \left ( s - \frac{2}{\xi}\right )\right ]</math>
 
 
where
 
 
:<math>s = \frac{E^{e^-}}{h \nu} =  \frac{E^{e^-}}{E_{\gamma}}</math>
 
:<math>\xi \approx 2 E_{\gamma}</math>
 
:<math>r_e^2 = 0.794</math> barns
 
:<math>m_e c^2 = 0.511</math> MeV
 
 
====GEANT 4 ====
 
 
GEANT 4 parametrized the Compton cross section to reproduce the data down to 10 keV using the expression
 
 
:<math>\sigma(Z,E_{\gamma}) = \left [ P_1(Z) \frac{\log(1+2 \xi)}{\xi} + \frac{P_2(Z) + P_3(Z) \xi + P_4(Z) \xi^2}{1 + a \xi + b \xi^2 + c \xi^3} \right ]</math>
 
 
 
: <math>P_i(Z) = d_i Z + e_i Z^2 + f_i Z^3</math>
 
: <math>1 \le Z \le 100</math>
 
:<math> a,b,c, d_i, e_i, f_i</math> are determined from fit
 
 
The data used in the fit may be found in
 
:Hubbell, Grimm, & Overbo, J. Phys. Chem. Ref. Data 9, (1980) pg 1023
 
:H. Storm, Nucl. Data Tables, A7 (1970) pg 565
 
 
In addition to the default and low energy models which come with GEANT4 (as was available with the Photo Electric effect), there is also a model called "G4LECS" which may be installed.
 
 
; Models
 
 
:a.) G4ComptonScattering: listed as "compt" in the process name.  No Rayleigh scattering in the model.
 
 
:b.) G4LowEnergyCompton:  Process name is "LowEnCompton" in the tracking code.  It has errors in the treatment of Rayleigh scattering and does not account for doppler broadening (the effect of bound electron momentum on the scattered particle energies).
 
 
:c.) G4LECS:  Bound electron effects are corrected for on an :shell-by-shell" basis.  Rayleigh scattering is modeled using the coherent scattering cross section and form factor data.  The Doppler broadening effect is included ( a result of the compton telescope simulation work).
 
 
;Note: Thomson & Rayleigh scattering are classical processes related to Compton scattering.  Klein-Nishina formula reduces to the Thomson cross section at low energies such that <math>\sigma_{Thompson} = \frac{8 \pi}{3} r_e^2</math>.  Thomson scattering produces polarized light because at these low non-relativistic energies the particle that absorbs the photon emits it in a direction perpendicular to its motion, that motion is the results of sees the oscillating E & M wave from the incident photon.  Rayleigh scattering (why our sky is blue) is photon scattering from an atom as a whole, coherent scattering.  No energy is transfered to the Medium in either case, <math>\gamma</math> only changes direction.  Rayleigh scattering is the elastic scattering of the photon from a particle that is '''smaller''' than the wavelength of the light.  Mei theory is used to describe elastic scattering from particles that are '''larger''' than the wavelength of the incident photon.
 
 
=== Pair Production ===
 
 
Pair production is similar to the Bremsstrahlung process.
 
 
 
 
Remember, in Bremsstrahlung the incident charged particle interacts with the <math>\vec{E}</math> of the Nucleus (or shell electron)
 
 
 
[[Image:SPIM_BremProcessDiagram.jpg | 400 px]]
 
 
 
In pair production a photon interacts with the <math>\vec{E}</math> of the Nucleus.
 
 
[[Image:SPIM_PairProductionProcessDiagram.jpg | 400 px]]
 
 
when the recoil of the atom is taken into account
 
 
:<math>E_{threshold} = 2 m_e(1+\frac{m_e}{M_A})</math> = Threshold energy for pair production from an atom of mass <math>M_A</math>
 
 
;Note: You can also have photon-electron pair production analogous to electron-electron bremsstrahlung production.
 
 
 
==== Pair Production Cross Section====
 
 
The pair production cross section is given by Equation 6.35 in  ([http://www.physics.isu.edu/~tforest/Classes/NucSim/Day11/PairProd_HighEnergy_Bethe_Phys.Rev._vol93_pg768_1954.pdf  Bethe, Phys. Rev., vol. 93 (1964) pg 768])
 
 
===== at small angles =====
 
 
A version which assumes small angles is given in Eq 7.35 of the same reference as the triple differential cross section:
 
 
:<math>\frac{d \sigma}{d \epsilon_1 d \theta_1 d \theta_2} = 8 \left ( \frac{\pi a}{\sinh (\pi a)} \right )^2 \frac{a^2}{2 \pi} \frac{e^2}{\hbar c} \left ( \frac{\hbar}{m_e c }\right )^2 \frac{\epsilon_1 \epsilon_2}{k^3}  \theta_1 \theta_2 </math>
 
: <math>\times \left \{ \frac{V^2(x)}{q^4} \left [ k^2 (u^2 + v^2) \xi \eta - 2 \epsilon_1 \epsilon_2 (u^2 \xi^2 + v^2 \eta^2 ) + 2 (\epsilon_1^2 + \epsilon_2^2)uv \xi \eta cos(\phi) \right ]  \right . </math>
 
:<math>\left . + a^2W^2(x) \xi^2 \eta^2 \left [ k^2(1 - (u^2+v^2)\xi \eta - 2 \epsilon_1 \epsilon_2 (u^2 \xi^2 + v^2 \eta^2) -2 (\epsilon_1^2 + \epsilon_2^2) u v \xi \eta \cos(\phi)\right ]\right \}</math>
 
 
where
 
 
:<math>k =</math> photon momentum/energy
 
:<math>\theta_1</math> = scattering angle of <math>e^+</math>
 
:<math>\theta_2</math> = scattering angle of <math>e^-</math>
 
: <math>\phi = \phi_1 - \phi_2 = \phi</math> angle between the <math>e^+</math> and <math>e^-</math> pair
 
:<math>\epsilon_1 = \sqrt{p_1^2 + m_e^2}</math> = Energy of the positron
 
:<math>\epsilon_2 = \sqrt{p_2^2 + m_e^2}</math> = Energy of the electron
 
:<math>u = \epsilon_1 \theta_1</math>
 
:<math>v=\epsilon_2 \theta_2</math>
 
:<math>\xi = \frac{1}{1+u^2}</math>
 
:<math>\eta= \frac{1}{1+v^2}</math>
 
:<math>q^2 = u^2 + v^2 + 2 u v \cos(\phi)</math>
 
: <math>x= 1-q^2 \xi \eta</math>
 
:<math>V(x) = 1 + \frac{a^2}{(1!)^2} +  \frac{a^2 (1+a^2) x^2}{(2!)^2} + \frac{a^2 (1+a^2)(2^2+a^2)x^4 x^2}{(3!)^2} + \cdots</math>
 
:<math>W(x) = \frac{1}{a^2} \frac{d V(x)}{d x}</math>
 
: <math>a = \frac{Ze^2}{\hbar c}</math>
 
 
;Note: The above equations for the differential cross section are using "natural" units where <math>c  \equiv 1</math>
 
 
 
===== Davies' version integrates over all angles =====
 
[http://www.physics.isu.edu/~tforest/Classes/NucSim/Day8/Davies_Phys.Rev._vol93_pg788_1954.pdf Davies] published a version which has been integrated over angles and includes some screening effects ( see Eq 35):
 
 
:<math>\frac{d \sigma}{d \epsilon_1} = 2 a^2 \frac{e^2}{\hbar c} \left ( \frac{\hbar}{m_e c }\right )^2  \frac{\epsilon_1^2 +\epsilon_2^2 + + \frac{2}{3} \epsilon_1 \epsilon_2 }{k^2} \left [ 2 \log \left (\frac{2\epsilon_1 \epsilon_2}{k^2} \right ) - 1 - 2 f(Z)\right ]</math>
 
: <math>= 4 Z^2 \alpha r_2^2 \frac{\epsilon_1^2 +\epsilon_2^2 + + \frac{2}{3} \epsilon_1 +\epsilon_2 }{( h \nu)^3} \left [ 2 \log \left (\frac{2\epsilon_1 \epsilon_2}{h \nu m_e c^2} \right ) - 1 - 2 f(Z)\right ]</math>
 
 
where
 
 
:<math>f(Z) = a^2 \sum_1^{\infty} \frac{1}{\nu (Z^2 +a^2)} \sim \frac{1}{1+a^2} + 0.20206 - 0.0369 a^2 + 0.0083 a^4 - 0.002 a^6</math>
 
 
If you integrate over all positron (<math>\epsilon_1</math> ) energies you get Eq. 44 (no screening)
 
 
:<math>\sigma_{e^+e^-} = 4 Z^2 \alpha r_e^2 \left [ \frac{7}{9} \right ( \ln(\frac{2 h \nu}{m_e c^2}) -f(Z) \left ) - \frac{763}{378} \right ]</math>
 
 
and Eq. 45 (complete screening)
 
 
:<math>\sigma_{e^+e^-} = 4 Z^2 \alpha r_e^2 \left [ \frac{7}{9} \right ( \ln(\frac{183}{Z^{\frac{1}{3}}}) -f(Z) \left ) - \frac{7}{378} \right ]</math>
 
 
Davies expressions were shown to work well at high energies (<math>E_{\gamma} > 88</math> MeV)
 
 
===== Overbo's low energy Cross sections =====
 
 
At low energies ( <math>E_\gamma < 5</math> MeV), [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day11/Overbo_Egamma.lt.5MeV_Phys.Rev.A8_pg668_1973.pdf Overbo] published an exact calculation in the case that the Atomic field is unscreened.  Overbo then provided a fit to the results of his calculation which he claims is valid to within 0.1% for the energy range <math>(3 MeV < E_{\gamma} < 5 MeV)</math>.  The fit is given in Eq. 7.1 of his paper as
 
 
:<math>\sigma_{e^+e^-} = \sigma_B \left ( 1 + a + \frac{b}{k-2}\right )</math>
 
 
where
 
 
 
: <math>a= -0.44 (\alpha Z)^2 - 0.07 (\alpha Z)^4</math>
 
: <math>b= 5.06 (\alpha Z)^2 - 2.1 (\alpha Z)^4</math>
 
:<math>\sigma_B = \alpha Z^2 r_e^2 \frac{2 \pi}{3} \left ( \frac{k-2}{k} \right ) ^3  \left [ 1 + \frac{\epsilon}{2} + \frac{23 \epsilon^2}{40} + \frac{11 \epsilon^3}{60} + \frac{29\epsilon^4}{960} \right ]=</math> Low Energy unscreened Born approximation total cross section for pair production
 
:<math> k = \frac{h \nu}{m_e c^2} = \frac{E_{\gamma}}{0.511 MeV}</math> = incident photon energy in units of the electron rest mass energy
 
:<math>\epsilon \equiv \frac{2k -4}{2+k+2\sqrt{2k}}</math>
 
 
===== Intermediate Energy Cross sections =====
 
 
For <math>5 MeV < E_{\gamma} < 80 MeV</math> the Gradstein semi-ephirical formula is used from G. White Gradstein, Natl. Bur. Standard., Circ 583 (1957) pg 1.
 
 
:<math>\sigma = \sigma_{BH} - \Delta_e + \frac{b^2}{k} \ln(k -0.57)</math>
 
 
where
 
 
:<math>\Delta_e</math> = empirical constant = 4.02 barns for Pb
 
:<math>b^2</math> = empirical constant = 16.8 barns for Pb
 
:<math>\sigma_{BH}</math> = [[Simulations_of_Particle_Interactions_with_Matter#Bremsstrahlung Bethe-Heitler cross section]]
 
 
===== Triplet production =====
 
; Triplet Production: identifies  photon-electron pair production.  The recoiling electron track adds to the two <math>e^+e^-</math> tracks making three total particle tracks  (kind of like pair production ionization).
 
 
A description for how to calculate the cross section for this process is illustrated for 1<math>0 < E_{\gamma} < 20 MeV</math>  in [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day11/Wright_5MeV.gt.Egamma.lt.20MeV_Phys.Rev.C36_pg562_1987.pdf L.E. Wright, Phys. Rev. C 36 (1987) pg 582].  Unfortunately, analytic expression is not given but one could construct tables of <math>Tq(92)</math> and <math>Tq(1)</math> in equation 65 of L.E. Wright's paper.
 
 
An example for  <math>E_{\gamma} =6 MeV</math> may be found in [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day11/Sud_Phys.Rev.A49_pg4624_1994.pdf Sud & Vargus, Phys. Rev. A49 (1994) pg 4624].
 
 
===== GEANT4 Pair production =====
 
 
GEANT4 uses the pair production cross section given in Tsai, Rev. Mod. Phys, vol 46 (1974) pg 815.
 
 
:<math>\frac{d \sigma}{d \epsilon} = Z(Z+ \eta) \alpha r_e^2 \left \{ \left [ \epsilon^2 +(1- \epsilon)^2 (\phi_1 - 4f(Z) ) \right ] + \frac{2}{3} \epsilon(1-\epsilon) (\phi_2 - 4f(Z) )\right \}</math>
 
 
where
 
 
:<math>E_{\gamma} =</math> energy of incident photon
 
:<math>E_{e^-} =</math> KE of create electron
 
:<math>\epsilon = \frac{E_{e^-} +m_ec^2}{ E_{\gamma}} =</math> fraction of <math>E_{\gamma}</math> taken away by <math>e^-</math>
 
:<math>\eta=</math> triplet production correction
 
:<math>f(Z) =</math> high energy coulomb correction from Davies above
 
:<math>\phi_1</math> & <math>\phi_2</math> = electron screening functions
 
 
The formula GEANT4 uses may also be found on pg 541 in [http://www.physics.isu.edu/~tforest/Classes/NucSim/Day11/G4_PairProd_Baro_Rad.Phys.Chem._vol44_pg531_1994.pdf J. Bono , Radiation Physics Chemistry, vol 44 (1994) pg 531]
 
 
:<math>\frac{d \sigma}{d\epsilon} = \frac{2}{3} Z(Z+ \eta) \alpha r_e^2 C_r  \left [ 2(\frac{1}{2} - \epsilon)^2 \phi_1(\epsilon) + \phi_2(\epsilon) \right ] </math>
 
 
where
 
 
:<math>\phi_1(\epsilon) = g_1(b) + g_0(\kappa)</math>
 
:<math>\phi_2(\epsilon) = g_2(b) + g_0(\kappa)</math>
 
:<math>g_1(b) = \frac{2}{3} - 2 \ln(1+b^2) - 6b \tan^{-1}(\frac{1}{b} - b^2 \left [ 4 -4b\tan^{-1}(\frac{1}{b} -3 \ln(1 + \frac{1}{b^2}\right]</math>
 
:<math>g_2(b) = \frac{11}{6} - 2 \ln(1+b^2) - 3b \tan^{-1}(\frac{1}{b} - \frac{b^2}{2} \left [ 4 -4b\tan^{-1}(\frac{1}{b} -3 \ln(1 + \frac{1}{b^2}\right]</math>
 
:<math>g_0(\kappa) =4 \ln ( \frac{R m_e c}{\hbar}) + 4f(Z) + F_0(\kappa,Z)</math>
 
:<math>F_0(\kappa,Z) = \left [ -0.1774 - 12.10 \alpha Z + 11.18 (\alpha Z)^2 \right ] \sqrt{\frac{2}{\kappa}} + \left [ 8.523 + 73.26 \alpha Z - 44.41(\alpha Z)^2 \right ] \frac{2}{\kappa}</math>
 
:<math>-  \left [ 13.52 + 121.1 \alpha Z - 96.41 (\alpha Z)^2 \right ] \left ( \frac{2}{\kappa} \right)^{2/3} + \left [ 8.946 + 62.05 \alpha Z - 63.41(\alpha Z)^2 \right ] \left ( \frac{2}{\kappa} \right )^2</math>= low energy coulomb correction
 
:<math>R</math> = screening radius (adjustable parameter)
 
:<math>\frac{\hbar}{m_ec} = 3.8616 \times 10^{-13} m</math> = Compton wavelength
 
:<math>b = \frac{Rm_e c}{\hbar} \frac{1}{2\kappa} \frac{1}{\epsilon(1-\epsilon)}</math>
 
:<math>\kappa = \frac{E_{\gamma}}{m_e c^2}</math>
 
 
An older simulation description with graphs of the cross-sections may be found below.
 
 
[[File:Bigg_Lighthill_conversion_1Mev-100MeV_SandiaPreprint_SC-RR-68-619.pdf]]
 
 
 
Below is the review of Modern Physics article basically summarizing all of the physics for pair production up to 1969
 
 
[[File:RevModPhys_PairPoduction_1969_MotzOlsonKoch.pdf]]
 
 
=== Contributions  as function of Z ===
 
The plot below shows the contributions of the three photon absorption physics processes as a function of the incident photon energy and the Z of the target material.
 
At low energy (keV), the photo-electric effect dominant while at high energies (> 1 MeV) pair production starts to dominate.  Compton scattering dominates in the intermediate energy region.
 
 
[[Image:SPIM_PhotoAbsorptionPhysicsProcess-vs-Z.jpg | 400 px]]
 
 
= Hadronic Interactions =
 
 
=== Proton Bremsstrahlung===
 
 
:<math>\sigma_{class}|_{90^o} = 2.1\times 10^{-31} cm^2/sterad</math> = cross section for dipole radiation emitted at 90 degrees with respect to incident beam of particles scattered in a Coulomb field.
 
 
[[File:ProtonBrem_Drell_Huang_PhysRev_v99_n3_1955_pg686.pdf]]
 
 
 
=== Pluto event generator===
 
 
[http://arxiv.org/abs/0708.2382 A ROOT based Hadronic Simulation package based on Pluto]
 
 
I installed Pluto V 5.14.1 on inca
 
 
I needed to set the environmental variables under tcsh
 
 
setenv ROOTSYS ~/src/ROOT/root
 
setenv PATH ${PATH}:${ROOTSYS}/bin
 
setenv LD_LIBRARY_PATH $ROOTSYS/lib
 
 
 
There is a subdirectory called "macros"
 
 
cd macros
 
 
Go to that subdirectory and type root, this will run the contents of the file "rootlogin.C"
 
 
cd macros
 
 
inca:~/src/Pluto/pluto_v5.14.1/macros> root
 
  *******************************************
 
  *                                        *
 
  *        W E L C O M E  to  R O O T      *
 
  *                                        *
 
  *  Version  5.17/03    30 August 2007  *
 
  *                                        *
 
  *  You are welcome to visit our Web site  *
 
  *          http://root.cern.ch            *
 
  *                                        *
 
  *******************************************
 
Compiled on 5 September 2007 for linux with thread support.
 
CINT/ROOT C/C++ Interpreter version 5.16.24, July 26, 2007
 
Type ? for help. Commands must be C++ statements.
 
Enclose multiple statements between { }.
 
  *********************************************************
 
  * The Pluto event generator                             
 
  * (C) HADES collaboration and all contributing AUTHORS 
 
  * www-hades.gsi.de/computing/pluto/html/PlutoIndex.html 
 
  * Version: 5.14.1
 
  * Compiled on 10 December 2008
 
  *********************************************************
 
Shared library Pluto.so loaded
 
 
to run a pp elastic model type
 
 
root [0] .x pp_elastic.C
 
 
a root ntuple is generate called "pp_elastic.root"
 
 
you can then analyze the data in the root file with
 
 
data->MakeClass();
 
 
the above command within root generates an analysis skeleton program.
 
 
using t.Show you can see the structure of the events within the ntuple.  A few functions are also stored in the root tree which you can use.
 
You can use the root file event to create an input file which GEANT4 can then use as its event generator.  GEANT4 then reads the events in and propagates them through your geometry.
 
 
=== Neutron Interactions ===
 
 
{| border="3"  cellpadding="20" cellspacing="0"
 
| Name || Energy
 
|-
 
|Cold Neutron ||  micro eV
 
|-
 
|Thermal || <math>\sim \frac{1}{40}</math> eV
 
|-
 
|epithermal || <math> \frac{1}{40}  eV \rightarrow 100 keV</math>
 
|-
 
|fast || <math>100 keV \rightarrow 100 MeV</math>
 
|-
 
|high energy || <math> > 100 MeV</math>
 
|-
 
|}<br>
 
 
'''Note:''' Interaction length for neutrons is ~<math>10^{-13}</math> .<br>
 
Neutrons are even better than photons for penetrating matter.<br>
 
 
==== Elastic scattering====
 
 
 
[[Image:Elastic_scattering_from_Nuclei.jpg]]<br>
 
 
<math>v_{CM} = \frac{m_n v_L + M(0)}{m_n + M} = \frac{v_0}{1 + \frac{M}{m_n}} = \frac{v_0}{1+A} =</math> velocity of CM frame<br>
 
 
<math>{v_L}^' = </math> Magnitude of Neutron velocity in CM frame before and after collision<br> <math>= v_c - v_{CM} = v_0 -\frac{v_0}{1+A} = \frac{(1+A)-1}{1+A} v_0 = \frac{A}{1+A} v_0</math><br>
 
 
<math> v = </math> Magnitude of Nucleus  velocity in CM frame before and after collision<br>
 
<math>= v_{CM} = \frac{v_0}{1+A}</math><br>
 
 
 
'''Note:'''  In elastic collision only the particles direction changes.<br>
 
 
<math>\vec{v}_L = {\vec{v}_L}^' + \vec{v}_{CM}</math> <br>
 
 
[[Image:Rule_of_cosines.jpg]]<br>
 
 
;;;<math>c^2 = a^2 + b^2 - 2abcos \theta</math><br>
 
 
<math>(v_L)^2 = ({v_L}^')^2 + (v_{CM})^2 - 2 v_{CM} {v_L}^' cos(\pi - {\theta}_{CM})=</math><br>
 
;;<math>= ({v_L}^')^2 + (V)^2 - 2 V {v_L}^' cos(\pi - {\theta}_{CM})=</math><br>
 
 
where<br>
 
 
;;<math>({v_L}^')^2 = (\frac{A}{1+A})^2 {v_0}^2</math><br>
 
;;<math>(V)^2 = (\frac{1}{1+A})^2 {v_0}^2</math><br>
 
 
After substitution we get following:<br>
 
 
<math>(\frac{v_L}{v_0})^2 = \frac{A^2 +1 - 2 A cos(\pi - {\theta}_{CM})}{(1+A)^2} = \frac{A^2 +1 + 2 A cos({\theta}_{CM})}{(1+A)^2}</math><br>
 
 
<math>
 
cos(A+/-B) = cosAcosB -/+ sinAsinB</math><br>
 
 
<math>\frac{E}{E_0} = \frac{A^2 + 1 + 2Acos({\theta}_{CM})}{(1+A)^2}</math>
 
 
when <math>{\theta}_{CM}=0</math>, <math>E_{max} = E_0</math>. <br>
 
 
 
<math>E_{min} = \frac{(A-1)^2}{(A+1)^2} E_0 = (\frac{A-1}{A+1}) E_0 =</math> Minimum energy of scattered Neutron in LAB frame.<br>
 
 
[[Image:Rule_of_cosines_1.jpg]]<br>
 
 
<math>({v_L}^')^2 = (v_L)^2 + (v_{CM})^2 - 2 v_{CM} v_L cos(\pi - {\theta}_{CM})=</math><br>
 
;;<math>= (v_L)^2 + (V)^2 - 2 V v_L cos({\theta}_{L})=</math><br>
 
 
<math>(\frac{Av_0}{1+A})^2 = {v_L}^2 + (\frac{v_0}{1+A})^2 - 2 v_L (\frac{v_0}{1+A}) cos({\theta}_L)</math><br>
 
 
After substituting <math>v_L</math><br>
 
 
<math>cos{\theta}_L = [\frac{A^2 + 1 + 2 A cos{\theta}_{CM}}{(1+A)^2} + (\frac{1}{1+A})^2 - (\frac{A}{1+A})^2] \times \frac{(1+A)^2}{2\sqrt{A^2 +1 + 2Acos{\theta}_{CM}}} = </math><br>
 
 
<math>= \frac{[A^2 +1 + 2Acos{\theta}_{CM} + 1 - A^2]}{2\sqrt{A^2 +1 + 2Acos{\theta}_{CM}}} = \frac{1 + Acos{\theta}_{CM}}{\sqrt{A^2 +1 + 2Acos{\theta}_{CM}}}
 
</math><br>
 
 
'''Note:'''  <math> {E_A}^{CM} = \frac{1}{2} M_A V^2 = \frac{1}{2} A m_n (\frac{v_0}{1+A})^2 = \frac{A}{(1+A)^2} \frac{m_n {v_0}^2}{2}= </math><br>
 
 
<math> = \frac{A}{(1+A)^2}E_0 = </math> Energy of recoil Nuclei in CM frame.<br>
 
 
'''Conservation of Energy:''' <math>E_0 = E + E_A</math> <br>
 
 
;;;<math> E_A = E_0 - E = E_0 - \frac{A^2 + 1 + 2Acos{\theta}_{CM}}{(1+A)^2} E_0 = </math><br>
 
 
;;;;<math>\frac{(1+A)^2 - (A^2 +1 + 2Acos{\theta}_{CM})}{(1+A)^2}E_0 =</math><br>
 
 
==== Inelastic Scattering====
 
 
==Resources==
 
[http://geant4.web.cern.ch/geant4/  GEANT4 Home Page]
 
 
[http://root.cern.ch ROOT Home page]
 
 
[http://conferences.fnal.gov/g4tutorial/g4cd/Documentation/WorkshopExercises/  Fermi Lab Example]
 
 
 
[http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html NIST Range Tables]
 
 
[http://ie.lbl.gov/xray/  X-ray specturm]
 
 
[[Installing_GEANT4.9.3_Fsim]]
 
 
===Building GEANT4.9.5===
 
 
[[TF_GEANT4.9.5]]
 
  
  
  
 
[[TF_SPIM_OLD]]
 
[[TF_SPIM_OLD]]

Latest revision as of 22:31, 15 January 2024

Class Admin

TF_SPIM_ClassAdmin

Homework Problems

HomeWork_Simulations_of_Particle_Interactions_with_Matter

Introduction

TF_SPIM_Intro

Energy Loss

TF_SPIM_StoppingPower

Ann. Phys. vol. 5, 325, (1930)

Interactions of Electrons and Photons with Matter

TF_SPIM_e-gamma

Hadronic Interactions

TF_SPIM_HadronicInteractions

Final Project

A final project will be submitted that will be graded with the following metrics:

1.) The document must be less than 15 pages.

2.) The document must contain references in a bibliography (5 points) .

3.) A comparison must be made between GEANT4's prediction and either the prediction of someone else or an experimental result(30 points).

4.) The graphs must be of publication quality with font sizes similar or larger than the 12 point font (10 points).

5.) The document must be grammatically correct (5 points).

6.) The document format must contain the following sections: An abstract of 5 sentences (5 points) , an Introduction(10 points), a Theory section (20 points) , if applicable a section describing the experiment that was simulated, a section delineating the comparisons that were made, and a conclusion( 15 points).

Resources

GEANT4 Home Page

ROOT Home page

Fermi Lab Example


NIST Range Tables

X-ray specturm

Installing_GEANT4.9.3_Fsim

Saving/restoring Random number seed

You save the current state of the random number generator with the command

/random/setSavingFlag 1

/run/beamOn 100

/random/saveThisRun

A file is created called

currentEvent.rndm

/control/shell mv currentEvent.rndm currentEvent10.rndm


You can restore the random number generator and begin generating random number from the last save time

/random/resetEngineFrom currentEvent.rndm

Building GEANT4.11

4.11.2

TF_GEANT4.11

Building GEANT4.10

4.10.02

TF_GEANT4.10.2

4.10.01

TF_GEANT4.10.1

Building GEANT4.9.6

TF_GEANT4.9.6

Building GEANT4.9.5

TF_GEANT4.9.5

An old version of Installation notes for versions prior to 9.5

Old Install Notes


Visualization Libraries:

OpenGL

DAWN


Coin3D

Compiling G4 with ROOT

These instruction describe how you can create a tree within ExN02SteppingVerbose to store tracking info in an array (max number of steps in a track is set to 100 for the desired particle)

G4CompileWRootforTracks

Using SLURM

http://slurm.schedmd.com/quickstart.html


https://rc.fas.harvard.edu/resources/documentation/convenient-slurm-commands/

simple batch script for one process job

create the file submit.sbatch below

#!/bin/sh
#SBATCH --time=1
cd src/PI
./PI_MC 100000000000000

the execute

sbatch submit.sbatch

check if its running with

squeue

to kill a batch job

scancel JOBID

On minerve

Sample script to submit 10 batch jobs.

the filename is minervesubmit and you run like

source minervesubmit
cd /home/foretony/src/GEANT4/geant4.9.5/Simulations/N02wROOT/batch
qsub submit10mil
qsub submit20mil
qsub submit30mil
qsub submit40mil
qsub submit50mil
qsub submit60mil
qsub submit70mil
qsub submit80mil
qsub submit90mil
qsub submit100mil

The file submit10mil looks like this

#!/bin/sh
#PBS -l nodes=1
#PBS -A FIAC
#PBS -M foretony@isu.edu
#PBS -m abe
#
source /home/foretony/src/GEANT4/geant4.9.5/geant4.9.6-install/bin/geant4.sh
cd /home/foretony/src/GEANT4/geant4.9.5/Simulations/N02wROOT/batch/10mil
../../exampleN02 run1.mac > /dev/null 


use

qstat

to check that the process is still running

use

qdel jobID#

if you want to kill the batch job, the jobID number shows up when you do stat.

for example

[foretony@minerve HW10]$ qstat
Job id                    Name             User            Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
27033.minerve             submit           foretony        00:41:55 R default        
[foretony@minerve HW10]$ qdel 27033
[foretony@minerve HW10]$ qstat

Definitions of Materials

File:MCNP Compendium of Material Composition.pdf

Minerve2 GEANT 4.10.1 Xterm error

On OS X El Capitan V 10.11.4 using XQuartz

~/src/GEANT4/geant4.10.1/Simulations/B2/B2a/exsmpleB2a


# Use this open statement to create an OpenGL view:
/vis/open OGL 600x600-0+0
/vis/sceneHandler/create OGL
/vis/viewer/create ! ! 600x600-0+0
libGL error: failed to load driver: swrast
X Error of failed request:  BadValue (integer parameter out of range for operation)
  Major opcode of failed request:  150 (GLX)
  Minor opcode of failed request:  3 (X_GLXCreateContext)
  Value in failed request:  0x0
  Serial number of failed request:  25
  Current serial number in output stream:  26



TF_SPIM_OLD