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

From New IAC Wiki
Jump to navigation Jump to search
 
(93 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===
+
[[TF_SPIM_ClassAdmin]]
−
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===
+
== Homework Problems==
−
http://wiki.iac.isu.edu/index.php/Forest_Class_Policies
+
[[HomeWork_Simulations_of_Particle_Interactions_with_Matter]]
−
 
 
−
=== 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 the theoretical description 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 which 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 output 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]]
 
−
 
 
−
== Overview ==
 
−
=== Particle Detection ===
 
−
A device detects a particle only after the particle transfers energy to the device.
 
−
 
 
−
Energy intrinsic to a device depends on the material used in a device
 
−
 
 
−
Consider a device made of some material with an average atomic number (<math>Z</math>)  at some temperature (<math>T</math>).  The material's atoms are in constant thermal motion (unless T = zero degrees Klevin).
 
−
 
 
−
Statistical Thermodynamics tells us that the canonical energy distribution of the atoms is given by the Maxwell-Boltzmann statistics such that
 
−
 
 
−
<math>P(E) = \frac{1}{kT} e^{-\frac{E}{kT}}</math>
 
−
 
 
−
<math>P(E)</math> represents the probability of any atom in the system having an energy <math>E</math> where
 
−
 
 
−
<math>k= 1.38 \times 10^{-23} \frac{J}{mole \cdot K}</math>
 
−
 
 
−
Note:  You may be more familiar with the Maxwell-Boltzmann distribution in the form
 
−
 
 
−
<math>N(\nu) = 4 \pi N \left ( \frac{m}{2\pi k T} \right ) ^{3/2} v^2  e^{-mv^2/2kT}</math>
 
−
 
 
−
where <math>N(v) \Delta v</math> would represent the molecules in the gas sample with speeds between <math>v</math> and <math>v + \Delta v</math>
 
−
 
 
−
==== Example 1: P(E=5 eV) ====
 
−
 
 
−
;What is the probability that an atom in a 12.011 gram block of carbon would have an energy of 5 eV?
 
−
 
 
−
First lets check that the probability distribution is Normailized; ie: does <math>\int_0^{\infty} P(E) dE =1</math>?
 
−
 
 
−
 
 
−
<math>\int_0^{\infty} P(E) dE = \int_0^{\infty} \frac{1}{kT} e^{-\frac{E}{kT}} dE = \frac{1}{kT} \frac{1}{\frac{1}{-kT}} e^{-\frac{E}{kT}} \mid_0^{\infty} = - [e^{-\infty} - e^0]= 1</math>
 
−
 
 
−
<math>P(E=5eV)</math> is calculated by integrating P(E) over some energy interval ( ie:<math> N(v) dv</math>).  I will arbitrarily choose 4.9 eV to 5.1 eV as a starting point.
 
−
 
 
−
 
 
−
<math>\int_{4.9 eV}^{5.1 eV} P(E) dE =  - [e^{-5.1 eV/kT} - e^{4.9 eV/kT}]</math>
 
−
 
 
−
<math>k= (1.38 \times 10^{-23} \frac{J}{mole \cdot K} )  =  (1.38 \times 10^{-23} \frac{J}{mole \cdot K} )(6.42 \times 10^{18} \frac{eV}{J})= 8.614 \times 10^{-5} \frac {eV}{mole \cdot K}</math>
 
−
 
 
−
assuming a room temperature of <math>T=300 K</math>
 
−
 
 
−
then<math>kT = 0.0258 \frac{eV}{mole}</math>
 
−
 
 
−
and
 
−
 
 
−
<math>\int_{4.9 eV}^{5.1 eV} P(E) dE =  - [e^{-5.1/0.0258} - e^{4.9/0.0258}] = 4.48 \times 10^{-83} - 1.9 \times 10^{-86} \approx 4.48 \times 10^{-83}</math>
 
−
 
 
−
or in other words the precise mathematical calculation of the probability may be approximated by just using the distribution function alone
 
−
 
 
−
<math>P(E=5eV) = e^{-5/0.0258} \approx 10^{-85}</math>
 
−
 
 
−
This approximation breaks down as <math>E \rightarrow 0.0258 eV</math>
 
−
 
 
−
Since we have 12.011 grams of carbon and 1 mole of carbon = 12.011 g = <math>6 \times 10^{23} </math>carbon atoms
 
−
 
 
−
We do not expect to see a 5 eV carbon atom in a sample size of <math>6 \times 10^{23} </math> carbon atoms when the probability of observing such an atom is <math>\approx 10^{-85}</math>.  Note: The mass of the earth is about <math>10^{27}</math> g <math>\approx 10^{50}</math> atoms, so a detector the size of the earth would also have a very small probability of having a carbon atom with an energy of 5 eV.
 
−
 
 
−
The energy we expect to see would be calculated by
 
−
 
 
−
<math><E> = \int_{0}^{\infty} E \cdot P(E) dE</math>
 
−
 
 
−
If you used this block of carbon as a detector you would easily notice an event in which a carbon atom absorbed 5 eV of energy as compared to the energy of a typical atom in the carbon block.
 
−
 
 
−
----
 
−
 
 
−
;Silicon detectors and Ionization chambers are two commonly used devices for detecting radiation.
 
−
 
 
−
approximately 1 eV of energy is all that you need to create an electron-ion pair in Silicon
 
−
 
 
−
<math>P(E=1 eV) = e^{-1/0.0258} \approx 10^{-17}</math>
 
−
 
 
−
approximately 10 eV of energy is needed to ionize an atom in a gas chamber
 
−
 
 
−
<math>P(E=10 eV) = e^{-10/0.0258} \approx 10^{-169}</math>
 
−
 
 
−
 
 
−
 
 
−
The low probability of having an atom with 10 eV of energy means that an ionization chamber would have a better Signal to Noise ratio (SNR) for detecting 10 eV radiation than a silicon detector
 
−
 
 
−
But if you cool the silicon detector to 200 degrees Kelvin (200 K) then
 
−
 
 
−
<math>P(E=1 eV) = e^{-1/0.0172} \approx 10^{-26} << 10^{-17}</math>
 
−
 
 
−
So cooling your detector will slow the atoms down making it more noticable when one of the atoms absorbs energy.
 
−
 
 
−
also, if the radiation flux is large, more electron-hole pairs are created and you get a more noticeable signal.
 
−
 
 
−
Unfortunately, with some detectore, like silicon, you can cause radiation damage that diminishes it's quantum efficiency for absorbing energy.
 
−
 
 
−
; What does this have to do with Simulations?
 
−
: You just did a SImulation.  Consider the following description of the Monte Carlo Method
 
−
 
 
−
=== The Monte Carlo method ===
 
−
; Stochastic
 
−
: from the greek word "stachos"
 
−
: a means of, relating to, or characterized by conjecture and randomness.
 
−
 
 
−
 
 
−
A stochastic process is one whose behavior is non-deterministic in that the next state of the process is partially determined.
 
−
 
 
−
The above particle detector was an example of describing a stochastic process using a probability distribution to determine the likely hood of finding an atom with a certain energy.
 
−
 
 
−
Quantun Mechanics is perhaps the best example such a non-deterministic systems.  The canonical systems in Thermodynamics is another example.
 
−
 
 
−
Basically the monte-carlo method uses a random number generator (RNG) to generate a distribution (gaussian, uniform, Poission,...) which is used to solve a stochastic process based on an astochastic description.
 
−
 
 
−
 
 
−
==== Example 2 Calculation of <math>\pi</math>====
 
−
 
 
−
;Astochastic description:
 
−
: <math>\pi</math> may be measured as the ratio of the area of a circle of radius <math>r</math> divided by the area of a square of length <math>2r</math>
 
−
 
 
−
[[Image:PI_from_AreaRatio.jpg]]<math>\frac{A_{circle}}{A_{square}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}</math>
 
−
 
 
−
You can measure the value of <math>\pi</math> if you physically measure the above ratios.
 
−
 
 
−
; Stochastic description:
 
−
: Construct a dart board representing the above geometry, throw several darts at it, and look at a ratio of the number of darts in the circle to the total number of darts thrown (assuming you always hit the dart board).
 
−
 
 
−
; Monte-Carlo Method
 
−
:Here is an outline of a program to calulate <math>\pi</math> using the Monte-Carlo method with the above Stochastic description
 
−
[[Image:MC_PI_fromAreaRatio.jpg]]
 
−
begin loop
 
−
  x=rnd
 
−
  y=rnd
 
−
  dist=sqrt(x*x+y*y)
 
−
  if dist <= 1.0 then numbCircHits+=1.0
 
−
  numbSquareHist += 1.0
 
−
end loop
 
−
  print PI = 4*numbCircHits/numbSquareHits
 
−
 
 
−
=== A Unix Primer ===
 
−
To get our feet wet using the UNIX operating system, we will try to solve example 2 above using a RNG under UNIX
 
−
 
 
−
==== List of important Commands====
 
−
 
 
−
# ls
 
−
# pwd
 
−
# cd
 
−
# df
 
−
# ssh
 
−
# scp
 
−
# mkdir
 
−
# printenv
 
−
# emacs, vi, vim
 
−
# make, gcc
 
−
# man
 
−
# less
 
−
# rm
 
−
 
 
−
----
 
−
Most of the commands executed within a shell under UNIX have command line arguments (switches) which tell the command to print information about using the command to the screen.  The common forms of these switches are "-h", "--h", or "--help"
 
−
 
 
−
ls --help
 
−
ssh -h
 
−
 
 
−
'' the switch deponds on your flavor of UNIX''
 
−
 
 
−
if using the switch doesn't help you can try the "man" (sort for manual) pages (if they were installed). 
 
−
Try
 
−
man -k pwd
 
−
 
 
−
the above command will search the manual for the key word "pwd"
 
−
 
 
−
==== Example 3: using UNIX to compile a RNG====
 
−
 
 
−
Step
 
−
# login to aztec.<br> [http://physics.isu.edu/~tforest/Classes/NucSim/Day1/XwindowsOnWindows.html click here for a description of logging in if using windows]
 
−
# mkdir src
 
−
# cd src
 
−
# cp -R ~tforest/NucSim/Day1 ./
 
−
# ls
 
−
# cd Day1
 
−
# make
 
−
#./rndtest
 
−
 
 
−
[http://physics.isu.edu/~tforest/Classes/NucSim/Day1/RNG/Marsaglia/noviceExample/ Here is a web link to the source files you can copy in case the above doesn't work]
 
−
 
 
−
=== A Root Primer ===
 
−
In bash shell do
 
−
 
 
−
export ROOTSYS=~tforest/src/ROOT/root
 
−
 
 
−
if chsh do
 
−
 
 
−
setenv ROOTSYS ~tforest/src/ROOT/root
 
−
 
 
−
To start the root program type
 
−
 
 
−
$ROOTSYS/bin/root
 
−
 
 
−
==== Example 1: Create Ntuple  and Draw Histogram====
 
−
 
 
−
=== Cross Sections ===
 
−
==== Definitions ====
 
−
;Total cross section
 
−
:<math>\sigma</math> = <math>\equiv \frac{\# particles\; scattered} {\frac{ \# incident \; particles}{Area}}</math>
 
−
 
 
−
;Differential cross section
 
−
:<math>\sigma(\theta)</math> = <math> \frac{d \sigma}{d \Omega} \equiv \frac{\frac{\# particles\; scattered}{solid \; angle}} {\frac{ \# incident \; particles}{Area}}</math>
 
−
 
 
−
; Solid Angle
 
−
:[[Image:SolidAngleDefinition.jpg]]
 
−
: <math>\Omega</math>= surface area of a sphere covered by the detector
 
−
: ie;the detectors area projected onto the surface of a sphere
 
−
:A= surface area of detector
 
−
:r=distance from interaction point to detector
 
−
:<math>\Omega = \frac{A}{r^2} </math>sterradians
 
−
: <math>A_{sphere} = 4 \pi r^2</math> if your detector was a hollow ball
 
−
:<math>\Omega_{max} = \frac{4 \pi r^2}{r^2} = 4\pi</math>sterradians
 
−
 
 
−
;Units
 
−
:Cross-sections have the units of Area
 
−
:1 barn = <math>10^{-28} m^2</math>
 
−
; [units of <math>\sigma(\theta)</math>] =<math>\frac{\frac{[particles]}{[sterradian]}} {\frac{ [ particles]}{[m^2]}} = m^2</math>
 
−
;Luminosity
 
−
:<math>L = \frac{Number of Scatterers}{Area \cdot time} \sim i_{beam} \rho_{target} l_{target}</math>
 
−
 
 
−
[[Image:FixedTargetScatteringCrossSection.jpg | 300 px]]
 
−
; Fixed target scattering
 
−
: <math>N_{in}</math>= # of particles in = <math>I \cdot A_{in}</math>
 
−
:: <math>A_{in}</math> is the area of the ring of incident particles
 
−
:<math>dN_{in} = I \cdot dA = I (2\pi b) db</math>= # particles in a ring of radius <math>b</math> and thickness <math>db</math>
 
−
 
 
−
You can measure <math>\sigma(\theta)</math> if you measure the # of particles detected <math>d N</math> in a known detector solid angle <math>d \Omega</math> from a known incident particle Flux (<math>I</math>)  as
 
−
 
 
−
<math>\sigma(\theta) = \frac{\frac{d N}{ d \Omega}}{I}</math>
 
−
 
 
−
Alternatively if you have a theory which tells you <math>\sigma(\theta)</math> which you want to test experimentally with a beam of flux <math>I</math> then you would measure counts (particles)
 
−
 
 
−
<math>dN = I \sigma(\theta) d \Omega = I \sigma(\theta)  \frac{d A}{r^2} = I \sigma(\theta) \frac{r^2 \sin(\theta) d \theta d \phi}{r^2}</math>
 
−
 
 
−
;Units
 
−
: <math>[d N] = [\frac {particles}{m^2}][m^2] [sterradian] </math> = # of particles
 
−
: or for a count rate divide both sides by time and you get beam current on the RHS
 
−
: integrate and you have the total number of counts
 
  
−
;Classical Scattering
+
=Introduction=
−
: In classical scattering you get the same number of particles out that you put in (no capture, conversion,..)
 
−
: <math>d N_{in} = dN</math>
 
−
:<math>d N_{in} = I dA = I (2\pi b) db</math>
 
−
: <math>d N = I \sigma(\theta) d \Omega =  I \sigma(\theta) \sin(\theta) d \theta d \phi = I \sigma(\theta) \sin(\theta) d \theta (2 \pi )</math>
 
−
:<math>  I (2\pi b) db =  I \sigma(\theta) \sin(\theta) d \theta (2 \pi )</math>
 
−
:<math>    b  db =  \sigma(\theta) \sin(\theta) d \theta </math>
 
−
:<math>\sigma(\theta) =  \frac{b}{\sin(\theta)}\frac{db}{d \theta}</math>
 
−
:<math>\frac{db}{d \theta}</math> tells you how the impact parameter <math>b</math> changes with scattering angle <math>\theta</math>
 
  
−
==== Example 4: Elastic Scattering ====
+
[[TF_SPIM_Intro]]
−
This example is an example of classical scattering.
 
  
−
Our goal is to find <math>\sigma(\theta)</math> for an elastic collision of 2 impenetrable spheres of diameter <math>a</math>.  We need to look for a relationship between the impact parameter<math> b</math> and the scattering angle<math> \theta</math>.  To find this relationship, let's solve this elastic scattering problem by  describing the collision using the Center of Mass (C.M.) coordinate system in terms of the reduced mass.  As we shall see, by using C.M. coordinate system the 2-body collision becomes a 1-body problem.  Then we will describe the motion of the reduced mass in the C.M. Frame.
+
= Energy Loss =
  
−
[[Image:SPIM_ElasCollis_Lab_CM_Frame.jpg | 500 px]]
+
[[TF_SPIM_StoppingPower]]
−
[[Media:SPIM_ElasCollis_Lab_CM_Frame.xfig.txt]]
 
  
−
; Variable definitions
 
−
:<math>b</math>= impact parameter ; distance of closest approach
 
−
:<math>m_1</math>= mass of incoming ball
 
−
:<math>m_2</math>= mass of target ball
 
−
:<math>u_1</math>= iniital velocity of  incoming ball in Lab Frame
 
−
:<math>v_1</math>= final velocity of  <math>m_1</math> in Lab Frame
 
−
:<math>\psi</math>= scattering angle of <math>m_1</math> in Lab frame after collision
 
−
:<math>u_1^{\prime}</math>= iniital velocity of  <math>m_1</math> in C.M.  Frame
 
−
:<math>v_1^{\prime}</math>= final velocity of  <math>m_1</math> in C.M. Frame
 
−
:<math>u_2^{\prime}</math>= iniital velocity of  <math>m_2</math> in C.M.  Frame
 
−
:<math>v_2^{\prime}</math>= final velocity of  <math>m_2</math> in C.M. Frame
 
−
:<math>\theta</math>= scattering angle of <math>m_1</math> in C.M. frame after collision
 
−
 
−
 
−
;Determining the reduced mass:
 
−
 
−
 
−
[[Image:SPIM_2Body-1BodyCoordSystem.jpg | 500 px]]
 
−
 
−
; vector definitions
 
−
:<math>\vec{r}_1</math> = a position vector pointing to the location of <math>m_1</math>
 
−
:<math>\vec{r}_2</math> = a position vector pointing to the location of <math>m_2</math>
 
−
:<math>\vec{R}</math> = a position vector pointing to the center of mass of the two ball system
 
−
:<math>\vec{r} \equiv \vec{r}_1 - \vec{r}_2</math> = the magnitude of this vector is the distance between the two masses
 
−
 
−
In the C.M. reference frame the above vectors have the following relationships
 
−
 
−
# <math>\vec{R} = 0 = \frac{m_1 \vec{r}_1 + m_2 \vec{r}_2}{m_1 + m_2} \Rightarrow m_1 \vec{r}_1 = -m_2 \vec{r}_2</math>
 
−
# <math>\vec{r}_1 - \vec{r}_2 = \vec{r}</math>
 
−
 
−
solving the above equations for <math>\vec{r_1}</math> and <math>\vec{r_2}</math> and defining the reduced mass <math>\mu</math> as
 
−
 
−
:<math>\mu = \frac{m_1 \cdot m_2}{m_1 + m_2} \equiv</math> reduced mass
 
−
 
−
leads to
 
−
 
−
: <math>\vec{r}_1 = \frac{\mu}{m_1} \vec{r}</math>
 
−
: <math>\vec{r}_2 = \frac{\mu}{m_2} \vec{r}</math>
 
−
 
−
We can use the above reduced mass relationships to construct the Lagrangian in terms of <math>\vec{r}</math> instead of <math>\vec{r}_1</math> and <math> \vec{r}_2</math> thereby reducing the problem from a 2-body problem to a 1-body problem.
 
−
 
−
; Construct the Lagrangian
 
−
 
−
The Lagrangian is defined as:
 
−
 
−
 
−
<math>\mathcal{L} = T - U</math>
 
−
 
−
where
 
−
 
−
<math>T \equiv</math> kinetic energy of the system
 
−
 
−
<math>U \equiv</math> Potential energy of the system which describes the interaction
 
−
 
−
 
−
<math>\mathcal{L} = \frac{1}{2} |\vec{\dot{r}}_1|^2 + \frac{1}{2} |\vec{\dot{r}}_2|^2 - U</math>
 
−
:= <math>\frac{1}{2} m_1 \left (\frac{m_2}{m1+m_2} \right )^2  |\vec{\dot{r}}|^2 + \frac{1}{2} m_2 \left (\frac{m_1}{m1+m_2} \right )^2  |\vec{\dot{r}}|^2 -U(\vec{r})</math>
 
−
:= <math>\frac{1}{2} \left ( m_2 + m_1 \right ) \left (\frac{m_1m_2}{(m1+m_2)^2} \right )  |\vec{\dot{r}}|^2 -U(\vec{r})</math>
 
−
 
−
after substituting derivative of the expressions for <math>\vec{r_1}</math> and <math>\vec{r}_2</math>
 
−
 
−
: = <math>\frac{1}{2} \mu |\vec{\dot{r}}|^2 -U(\vec{r})</math>
 
−
 
−
The 2-body problem is now described by a 1-body Lagrangian we need to determine which coordinate system (cartesian, spherical,..) to use to write an expression for (<math>|\vec{\dot{r}}|^2</math>).  Polar seems best unless there is a dependence in the azimuthal angle.
 
−
 
−
Lagranges equations of motion are given by
 
−
: <math>\frac{\partial \mathcal{L}}{\partial q} = \frac{d}{dt} \frac{\partial \mathcal{L}}{\dot{q}}</math>
 
−
where <math>q</math> represents one of the coordinate (cannonical variables).
 
−
 
−
To get the classical scattering cross section we are interested in finding an expression for the dependence of the impact parameter on the scattering angle,<math>\frac{d b}{d \theta}</math>. 
 
−
 
−
Now lets redraw the collision in terms of a reference frame fixed on <math>m_2</math> (before collision its the Lab Frame but not after collision).
 
−
 
−
[[Image:SPIM_ElasCollis_CMFrame.jpg]] [[Media:SPIM_ElasColls_CMFrame_xfig.txt]]
 
−
 
−
The C.M. Frame rides along the center of mass, the above coordinate system though has its origin on <math>m_2</math>.  The above drawing identifies <math>\theta</math> and <math>\phi</math> for the system at the point of the collision in which the CM frame is a distance <math>a</math> (the size of the ball) from the origin of the coordinate system fixed to <math>m_2</math>.  If <math>b > a</math> then there is no collision (<math>\theta=0</math>),  otherwise a collision happens when r=a (the distance between the balls is equal to their diameter).  A head on collision is defined as <math>b=0</math> (<math>\theta=\pi</math>).
 
−
 
−
;Observation
 
−
: as <math>\theta</math> gets smaller, <math>b</math> gets bigger
 
−
: <math>\frac{d b}{d \theta} < 0</math>
 
−
 
−
Using plane polar coordinates (<math>r, \phi</math>) we can describe the problem in the lab frame as:
 
−
 
−
<math>\vec{v} = \dot{R} \hat{e}_r + r \dot{\phi} \hat{e}_{\phi}</math>
 
−
 
−
<math>T = \frac{1}{2} \mu ( \dot{r)}^2 + r^2 \dot{\phi}^2)</math>
 
−
 
−
 
−
<math>U(r) = \left \{  {0 \; r > a \atop \infty \; r \le a} \right .</math>
 
−
 
−
<math>\mathcal{L} = T -U = \frac{1}{2} \mu ( \dot{r}^2 + r^2 \dot{\phi}^2) - U(r)</math>
 
−
 
−
Lagranges Equation of Motion:
 
−
 
−
<math>\frac{\partial \mathcal {L}}{\partial \phi} = \frac{d}{d t} \frac{\partial \mathcal{L}}{\partial \dot{\phi}}</math>
 
−
<math>0 = \frac{d}{d t} [ \mu r^2 \dot{\phi}] \Rightarrow</math>  there is a constant of motion ( Constant angular momentum)
 
−
 
−
<math>\ell \equiv \mu r^2 \dot{\phi} = \vec{r} \times \vec{p} = \vec{r} \times \mu \vec{v} = r^2 \mu \dot{\phi}</math>
 
−
 
−
substitute <math>\ell</math> into <math>\mathcal{L}</math>
 
−
 
−
<math>\mathcal{L} = \frac{1}{2} ( \mu  \dot{r}^2 + \frac{\ell}{\mu r^2} ) - U(r)</math>
 
−
 
−
The two equations above are in terms of <math>r</math> and <math>\phi</math> whereas our goal is to find an expression for <math>\frac{ d b}{ d \theta}</math>.  Since <math>r</math> is related to <math>b</math> and <math>\phi</math> is related to<math> \theta</math> (<math>\theta = \pi - 2\phi</math>; see figure above) we should try and find expressions for <math>d \phi</math> in terms of <math>r(b)</math>
 
−
 
−
;Trick
 
−
: <math>\dot{\phi} = \frac{d \phi}{d t} = \frac{d \phi}{d r} \frac{d r}{d t}</math>
 
−
:<math>\Rightarrow \ell = \mu r^2 \frac{d \phi}{d r} \dot{r}</math>
 
−
:or
 
−
: <math>d \phi = \frac{\ell}{\mu r^2 \dot{r}} dr</math>
 
−
 
−
We now need an expression for <math>\dot{r}</math> in order to integrate the above equation to determine the functional dependence of <math>\phi</math> and hence<math> \theta</math>.
 
−
 
−
Since Energy is conserved (Elastic Scattering), we may define the Hamiltonian as
 
−
 
−
<math>H = T + U = \frac{1}{2} (mu \dot{r}^2 + \frac{\ell}{\mu r^2}) + U(r) = constant \equiv E</math>
 
−
 
−
solving for <math>\dot{r}</math>
 
−
 
−
<math>\dot{r} = \pm \sqrt{\frac{2(E-U(r))}{\mu} - \frac{\ell^2}{\mu^2 r^2}}</math>
 
−
 
−
substituting the above into the equation for <math>d \phi</math> and integrating:
 
−
 
−
<math>\phi = \int d \phi = \int_{r_{min}}^{r_{max}} \frac{\ell}{\mu r^2 \dot{R}} dr</math>
 
−
 
−
<math>r_{min} = a  \; \; \;  r_{max}= \infty  \; \; \;  U(r) = 0 : a \le r \le \infty</math>
 
−
 
−
<math>\phi = \int_a^{\infty} \frac{\ell} {r^2 \sqrt{2 \mu E - \frac{\ell^2}{r^2}} }dr</math>
 
−
 
−
For <math>a \le r \le \infty</math> : <math>E = \frac{1}{2} \mu v^2_{cm} \Rightarrow v_{cm} = \sqrt{\frac{2E}{\mu}}</math>
 
−
 
−
<math>\vec{\ell} = \vec{r} \times \vec{p} \Rightarrow |\vec{\ell}| = |\vec{r}| |\vec{p}| \sin(\phi) = r \mu v_{cm} \sin(\phi) = r \mu \left ( \sqrt{\frac{2E}{\mu}} \right) \sin(\phi) = \sqrt{2 \mu E} r\sin(\phi) =\sqrt{2 \mu E} b</math>
 
−
 
−
substituting this expression for <math>\ell</math> into the last expression for <math>\phi</math> above :
 
−
 
−
<math>\phi =\int_a^{\infty}  \frac{b dr}{r\sqrt{(r^2-b^2)}}</math>
 
−
 
−
;Integral Table
 
−
: <math>\int  \frac{dx}{x\sqrt{(\alpha x^2+\beta x+\gamma)}} = \frac{-1}{\sqrt{-\gamma}} \sin^{-1} \left (\frac{\beta x+2\gamma}{|x|\sqrt{\beta^2-4\alpha \gamma}} \right )</math>
 
−
 
−
let <math>x=r \;\; \alpha=1 \;\; \beta=0 \;\; \gamma=-b^2</math>
 
−
 
−
then
 
−
 
−
<math>\phi = \left . b \frac{1}{\sqrt{-(-b^2)}} \sin^{-1} \left (\frac{-2b^2}{r\sqrt{0-4(1)(-b^2) } }\right ) \right |_a^{\infty} = \sin^{-1} (0)- \sin^{-1}(-\frac{b}{a})</math>
 
−
 
−
or
 
−
 
−
<math>\sin(\phi) = \frac{b}{a} = \sin \left ( \frac{\pi}{2} - \frac{\theta}{2} \right ) = \cos \left ( \frac{\theta}{2} \right )</math>
 
−
 
−
:<math>\Rightarrow b = a \cos \left( \frac{\theta}{2} \right)</math>
 
−
 
−
 
−
; Now substitute the above into the expression for <math>\sigma(\theta)</math>
 
−
 
−
<math>\sigma(\theta) = \frac{b}{\sin(\theta)} \frac{d b}{d \theta} = \frac{a \cos(\theta/2)}{sin(\theta)} a[-\sin(\theta/2)]\frac{1}{2}
 
−
= \frac{a^2}{2} \frac{\cos(\theta/2) \sin(\theta/2)}{\sin(\theta)}</math>
 
−
 
−
drop the negative sign, sqrt in denominator allows this, and use the trig identity
 
−
 
−
:<math>\sin \left (\frac{\theta}{2} + \frac{\theta}{2} \right ) = \cos \left (\frac{\theta}{2} \right) \sin \left (\frac{\theta}{2} \right ) + \cos \left ( \frac{\theta}{2} \right ) \sin \left (\frac{\theta}{2} \right )</math>
 
−
:<math>\sin(\theta) = 2 \cos \left (\frac{\theta}{2} \right ) \sin \left (\frac{\theta}{2} \right )</math>
 
−
 
−
<math>\sigma(\theta) = \frac{a^2}{2} \frac{1}{2} = \frac{a^2}{4}</math>
 
−
 
−
<math>\sigma = \int \sigma(\theta) d \Omega = \frac{a^2}{2} \frac{1}{2} 4 \pi  = \pi a^2</math>
 
−
 
−
 
−
;compare with result from definition
 
−
 
−
:<math>\sigma</math> = scattering cross-section <math>\equiv \frac{\# particles\; scattered} {\frac{ \# incident \; particles}{Area}}</math>
 
−
:number of particles scattered = number of incident particles
 
−
: Area = <math> \pi a^2</math> = The area profile in which a collision occurs( the ball diameter is <math>a</math>) [[Image:ClassicalEffectiveScatteringArea.jpg | 200 px]]
 
−
 
−
<math>\sigma = \frac{{N}}{\frac{ N}{\pi a^2}} = \pi a^2</math>
 
−
 
−
==== Lab Frame Cross Sections ====
 
−
 
−
The C.M. frame is often chosen to theoretically calculate cross-sections even though experiments are conducted in the Lab frame.  In such cases you will need to transform cross-sections between two frames.
 
−
 
−
The total cross-section should be frame independent
 
−
 
−
:<math>\sigma_{C.M.} = \sigma_{Lab}</math>
 
−
 
−
or
 
−
 
−
: <math>\sigma(\theta) d \Omega = \sigma(\psi) d \Omega^{\prime}</math>
 
−
 
−
where
 
−
 
−
<math>\theta</math> is in the CM frame and <math>\psi</math> is in the Lab frame.
 
−
 
−
 
−
;A non-relativistic transformation:
 
−
 
−
: <math>\sigma(\theta) d \Omega = \sigma(\psi) d \Omega^{\prime}</math>
 
−
: <math>\sigma(\theta) 2 \pi \sin(\theta) d \theta = \sigma(\psi) 2 \pi \sin (\psi) d \psi</math>
 
−
: <math>\Rightarrow \sigma(\psi) = \frac{\sin(\theta)}{\sin(\psi)} \frac{d \theta}{d \psi} \sigma(\theta)</math>
 
−
 
−
The transformation is governed by the dependence of <math>\theta</math> on <math> \psi</math> <math> \left( \frac{d \theta}{d \psi} \right )</math>
 
−
 
−
Lets return back to our picture of the scattering Process
 
−
 
−
[[Image:SPIM ElasCollis Lab CM Frame.jpg | 500 px]]
 
−
 
−
if we superimpose the vectors <math>\vec{v}_1</math> and <math>\vec{v}_1^{\prime}</math> we have
 
−
 
−
[[Image:SPIM ElasCollis Lab CM Frame_Velocities.jpg]]
 
−
 
−
Trig identities (non-relativistic Gallilean transformation) tell us
 
−
 
−
<math>v_1 \sin(\psi) = v_1^{\prime} \sin(\theta)</math>
 
−
 
−
 
−
<math>v_1 cos(\psi) = v_{cm} + v_1^{\prime} \cos(\theta)</math>
 
−
 
−
solving for <math>\psi</math>
 
−
 
−
<math>\tan(\psi) = \frac{\sin(\psi)}{\cos(\psi)} = \frac{v_1^{\prime} \sin(\theta)/v_1}{\frac{v_{CM}}{v_1} + \frac{v_1^{\prime} \cos(\theta)}{v_1} }
 
−
= \frac{\sin(\theta)}{\cos(\theta) + \frac{v_{CM}}{v_1^{\prime}}}</math>
 
−
 
−
For an elastic collision only the directions change in the CM Frame: <math>u_1^{\prime}= v_1^{\prime}</math>  & <math>u_1^{\prime}= v_2^{\prime}</math>
 
−
 
−
;From the definition of the C.M.
 
−
;<math>\vec{v}_{CM} = \frac{m_1 \vec{u}_1 + m_2 \vec{u}_2}{m_1+m_2} = \frac{m_1}{m_1+m_2} \vec{u}_1</math>
 
−
 
−
;conservation of momentum in CM Frame <math>\Rightarrow</math> :
 
−
:<math>m_1 u_1^{\prime} = - m_2 u_2{\prime}</math>
 
−
 
−
:<math> \Rightarrow v_1^{\prime} = u_1^{\prime} = \frac{-m_2}{m_1} u_2^{\prime}</math>
 
−
 
−
; Gallilean Coordinate transformation:
 
−
;<math>\vec{u}_1 = \vec{u}_1^{\prime} + \vec{v}_{CM} = \vec{u}_1^{\prime} + \frac{m_1}{m_1+m_2} \vec{u}_1</math>
 
−
:<math>\Rightarrow u_1{\prime} = \left [ 1 - \frac{m_1}{m_1 + m_2} \right ] u_1 = \frac{m_2}{m_1+m_2}u_1</math>
 
−
:<math>\Rightarrow v_1^{\prime} = u_1^{\prime}  =\frac{m_2}{m_1+m_2} u_1</math>
 
−
 
−
; another expression for <math>\psi</math>
 
−
 
−
using the above gallilean transformation we can do the following
 
−
 
−
:<math>\frac{v_{CM}}{v_1^{\prime}}= \frac{\frac{m_1}{m_1+m_2} u_1}{\frac{m_2}{m_1+m_2} u_1} = \frac{m_1}{m_2}</math>
 
−
 
−
or
 
−
 
−
: <math>\tan(\psi) = \frac{\sin(\theta)}{\cos(\theta) + \frac{m_1}{m_2}}</math>
 
−
 
−
after a little trig substitution
 
−
 
−
<math>\Rightarrow \frac{m_1}{m_2} = \frac{sin(\theta - \psi)}{\sin(\psi)} =</math> constant
 
−
 
−
now use the chain rule to find <math>\frac{d \theta}{d \psi}</math>
 
−
 
−
: <math>f \equiv  \frac{sin(\theta - \psi)}{\sin(\psi)} =</math> constant
 
−
:<math>df = 0 = \frac{ \partial f}{\partial \psi} d \psi  + \frac{ \partial f}{\partial \theta} d \theta </math>
 
−
: <math>\Rightarrow \frac{d \theta}{d \psi} = \frac{-\frac{ \partial f}{\partial \psi} }{\frac{ \partial f}{\partial \theta} }</math>
 
−
 
−
:<math>-\frac{ \partial f}{\partial \psi} = \frac{\cos(\theta - \psi)}{\sin(\psi)} + \frac{\sin(\theta - \psi)}{\sin(\psi)}</math>
 
−
:<math>\frac{ \partial f}{\partial \theta }= 1 + \frac{\sin(\theta - \psi) \cos(\psi)}{\cos(\theta - \psi) \sin(\psi)}</math>
 
−
 
−
after substitution:
 
−
: <math>\sigma(\psi) = \frac{\sin(\theta)}{\sin(\psi)} \frac{d \theta}{d \psi} \sigma(\theta)</math>
 
−
: <math>=\frac{\sin(\theta)}{\sin(\psi)} \left [ 1 + \frac{\sin(\theta - \psi) \cos(\psi)}{\cos(\theta - \psi) \sin(\psi)} \right ] \sigma(\theta)</math>
 
−
 
−
For the above equation to be more useful one would prefer to recast it in terms of only <math>\psi</math> and masses.
 
−
 
−
:<math>\sigma(\psi) = \frac{\left [ \frac{m_1}{m_2}\cos(\psi) + \sqrt{1-\left ( \frac{m_1 \sin(\psi) }{m_2} \right )^2 }\right ]}{\sqrt{1 - \left ( \frac{m_1 \sin(\psi)}{m_2}\right )^2 }}\sigma(\theta)</math>
 
−
 
−
== Stopping Power ==
 
 
Ann. Phys. vol. 5, 325, (1930)
 
Ann. Phys. vol. 5, 325, (1930)
−
=== Bethe Equation ===
 
−
====Classical Energy Loss ====
 
  
−
Consider the energy lost when a particle of charge (<math>ze</math>) traveling at speed <math>v</math> is scattered by a target of charge (<math>Ze</math>).  Assume only the coulomb force causes the particle to scatter from the target as shown below.
+
=Interactions of Electrons and Photons with Matter=
  
−
[[Image:SPIM_Bethe_ClassCoulScat.jpg]]
+
[[TF_SPIM_e-gamma]]
  
−
; Notice:
+
= Hadronic Interactions =
−
: as <math>ze</math> is scattered the horizontal component of the coulomb force (<math>F</math>) flips direction; ie net horizontal force for the scattering
 
  
−
:<math>F_{vertical} = k \frac{zZe^2}{r^2} \sin(\theta) = k \frac{zZe^2}{r^2}  \frac{b}{r}</math>
+
[[TF_SPIM_HadronicInteractions]]
  
−
where
+
= Final Project=
−
: k =<math>\frac{1}{4 \pi \epsilon_0}</math>
 
−
: r = distance between incident projectile and target atom
 
−
: b= impact parameter of collision
 
  
 +
A final project will be submitted that will be graded with the following metrics:
  
−
Using the definition of Impulse one can determine the momentum change of <math>ze</math> as
+
1.) The document must be less than 15 pages.
  
−
: <math>\Delta p = \int F dt</math>
+
2.) The document must contain references in a bibliography (5 points) .
  
−
Let's assume that the energy lost by the incident particle <math>ze</math> is absorbed by an electron in the target atom. This energy may be cast in terms of the incident particles momentum change as
+
3.) A comparison must be made between GEANT4's prediction and either the prediction of someone else or an experimental result(30 points).
  
−
:<math>\frac{(\Delta p)^2}{2m_e}</math>
+
4.) The graphs must be of publication quality with font sizes similar or larger than the 12 point font (10 points).
  
−
By calculating the change in momentum (<math>\Delta p</math>) of the incident particle we can infer that the energy lost by the incident particle is absorbed by one of the target material's atomic electrons.
+
5.) The document must be grammatically correct (5 points).
  
−
:<math>\Delta P = \int F dt = \int k \frac{zZe^2b}{r^3}  dt</math>
+
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).
  
−
using  <math>dt = \frac{dx}{v} = \frac{d x}{\beta c}</math> we have
+
=Resources=
  
−
: <math>= k \frac{zZe^2b}{\beta c} \int_{-\infty}^{+\infty} \frac{ dx}{(x^2+b^2)^{3/2}}</math>
+
[http://geant4.web.cern.ch/geant4/ GEANT4 Home Page]
−
: <math>=\frac{kzZe^2b}{\beta c b^2} \int_{-\infty}^{+\infty} \frac{ dx/b}{(1+\frac{x^2}{b^2})^{3/2}}</math>
 
  
−
:<math>\int_{-\infty}^{+\infty} \frac{ dx/b}{(1+\frac{x}{b^2})^{3/2}}=2</math>
+
[http://root.cern.ch ROOT Home page]
  
−
:<math> \Delta p = \frac{2kzZe^2b}{\beta c b^2}</math>
+
[http://conferences.fnal.gov/g4tutorial/g4cd/Documentation/WorkshopExercises/  Fermi Lab Example]
  
−
casting this in terms of the classical atomic electron radius <math>r_e</math>
 
  
−
: <math>r_e = \frac{k e^2}{m_e v^2} \sim \frac{k e^2}{m_e c^2}</math>  just equate <math>F = \frac{ke^2}{r_e^2} = m \frac{v^2}{r_e}</math>
+
[http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html NIST Range Tables]
  
−
Then
+
[http://ie.lbl.gov/xray/  X-ray specturm]
  
−
:<math> \Delta p = \frac{2zZr_e m_e c}{\beta  b}</math>
+
[[Installing_GEANT4.9.3_Fsim]]
  
−
and
+
== Saving/restoring Random number seed==
  
−
: <math>\Delta E = \frac{(\Delta p)^2}{2m_e} = 2 \left ( \frac{r_e m_e}{\beta b}\right )^2 \frac {z^2 Z^2 c^2}{m_e}</math> : <math>Z</math> = 1 here because I shall assume the energy is lost to just the electron and the Atom is a spectator
+
You save the current state of the random number generator with the command
  
−
Now let's calculate an expression representing the AVERAGE energy lost for an incident particle traversing a material of some thickness.
+
/random/setSavingFlag 1
  
−
Let
+
/run/beamOn 100
  
−
: <math>P(\Delta E)</math> = Probability of an interaction taking place which results in an energy loss <math>\Delta E</math>
+
/random/saveThisRun
  
−
If we let
+
A file is created called
  
−
Z  = Atomic Number = # electrons in target Atom = number of protons in an Atom
+
currentEvent.rndm
−
 
 
−
N = Avagadros number  = <math>6.022 \times 10^{23}  \frac{Atoms}{mol}</math>
 
−
 
 
−
A = Atomic mass =  <math>\frac{g}{mole}</math>
 
−
 
 
−
<math>dP(\Delta E)</math>  = probability of hitting an atomic electron in the area of an annulus of radius (<math>b + db</math>) with an energy transfer between <math>\Delta E</math> and <math>\Delta E + d(\Delta E)</math>
 
−
 
 
−
Then
 
−
 
 
−
:  <math>\frac{-dE }{dx}= \int_0^{\infty} dP(\Delta E) \Delta E</math> = energy lost by the incident particle per distance traversed through the material
 
−
 
 
−
I am just adding up all the energy losses weighted by the probability of the energy loss to find the total energy loss.
 
−
 
 
−
;What is <math>dP(\Delta E)</math> :
 
−
: <math>dP(\Delta E)</math> = probability of an energy transfer taking place = probability of an interaction = <math>\frac{N}{A} d \sigma</math>  [ Atoms <math>m^2</math>/g]
 
−
 
 
−
;<math>dP(\Delta E) = \frac{N}{A} d \sigma =\frac{N}{A} (2 \pi b db) Z</math>
 
−
:classically <math>\sigma = \pi b^2 ; d \sigma = 2\pi b db</math>
 
−
:In practice <math> \sigma</math> is a measured cross-section which is a function of energy.
 
−
 
 
−
<math>\Rightarrow \frac{-dE}{dx} =  \int_0^{\infty} \frac{N}{A} (2 \pi b db) Z \Delta E = \frac{2 \pi N Z}{A} \int_0^{\infty} \Delta E b db</math>
 
−
:= <math>\frac{2 \pi N Z}{A} \int_0^{\infty} \left [ \frac{2 r_e^2 m_e c^2 z^2}{\beta^2 b^2}\right ] b db</math>
 
−
:= <math>4 \pi N r_e^2 m_e c^2 \frac{z^2 Z}{A \beta^2} \int_0^{\infty} \frac{db}{b}</math>
 
−
:=<math>\frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \int_0^{\infty} \frac{db}{b}</math>
 
−
 
 
−
where
 
−
<math>\frac{\mathcal{K}}{A} = \frac{4 \pi N r_e^2 m_e c^2}{A} = 0.307 \frac{MeV cm^2}{g}</math>''' if A=1'''
 
−
 
 
−
The limits of the above integral should be more physical in order to reflect the limits of the physics interaction.  Let b_{min}  and b_{max} represent the minimum and maximum possible impact parameter where the physics is discribed,  as shown above,  by the coulomb force.
 
−
 
 
−
;What is <math>b_{min}</math>?
 
−
 
 
−
if <math>b \rightarrow 0</math> then <math>\frac{d E}{dx}</math> diverges and the energy transfer <math>\rightarrow \infty : \Delta E \sim \frac{1}{b}</math>.  Physically there is a maximum energy that may be transferred before the physics of the problem changes (ie: nuclear excitation, jet production, ...).  The de Borglie wavelength of the atom is used to estimate a value for <math>b_{min}</math> such that
 
−
 
 
−
: <math>b_{min} \sim  \frac{1}{2} \lambda_{de Broglie} = \frac{h}{2p} = \frac{h}{2 \gamma m_e \beta c}</math>
 
−
 
 
−
;What is <math>b_{max}</math>?
 
−
 
 
−
As <math>b</math> gets bigger the interaction is "softer" and longer.  If the interaction time (<math>\tau_i</math>) is so long that it is equivalent to an electron orbit (<math>\tau_R</math>) then the atom looks more like it is neutrally charged.  You move from an interaction in which the electron orbit is perturbed adiabatically such that there is no orbit change and the minimum amount of energy is transferred to no interaction taking place because the atom is neutral.
 
−
 
 
−
Let
 
−
 
 
−
: <math>\tau_i = \frac{b_{max}}{v} (\sqrt{1-\beta^2})</math>  : fields at high velocities get Lorentz contracted
 
−
: <math>\tau_R \equiv \frac{h}{I}</math> : I <math>\equiv</math> mean excitation energy of target material ( <math>E = h \nu = h/ \tau</math>)
 
−
 
 
−
Condition for <math>b_{max}</math> :
 
−
:<math>\tau_i = \tau_R</math>
 
−
<math>\Rightarrow b_{max} = \frac{h \gamma \beta c}{I}</math>
 
−
 
 
−
<math>-\frac{dE}{dx} = \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \int_0^{\infty} \frac{db}{b}</math>
 
−
: <math>= \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \ln \frac{b_{max}}{b_{min}}</math>
 
−
: <math>= \frac{\mathcal{K} }{A} \frac{z^2 Z}{\beta^2} \ln \frac{2 \gamma^2 m_e \beta^2 c^2}{I}</math>
 
−
 
 
−
=====Example 5: Find <math>\frac{dE}{dx}</math> for a 10 MeV proton hitting a liquid hydrogen (<math>LH_2</math>) target=====
 
−
A = Z=z=1<br>
 
−
<math>m_e c^2</math> = 0.511 MeV <br>
 
−
I = 21.8 eV : see gray data point for Liquid <math>H_2</math> From  Figure 27.5 on pg 6 of [http://pdg.lbl.gov/2007/reviews/passagerpp.pdf PDG] below.<br>
 
−
[[Image:PDG_IonizationPotential.jpg | 500 px]]
 
−
 
 
−
Just need to know <math>\gamma</math> and <math>\beta</math>
 
−
 
 
−
"a 10 MeV proton" <math>\Rightarrow</math> Kinetic Energy (K.E.) = 10 MeV = <math>(\gamma - 1) mc^2</math>
 
−
 
 
−
:<math>\Rightarrow \gamma = \frac{K.E.}{mc^2} + 1 = \frac{10 MeV}{938 MeV} + 1 \sim 1 = \frac{1}{\sqrt{1-\beta^2}}</math>
 
−
 
 
−
Proton is not relativistic
 
−
 
 
−
: <math>v^2 = \frac{2 K.E.}{m} = \frac{2 \cdot 10 MeV}{938 MeV/c^2} = 2 \times 10^{-2} c^2 \Rightarrow \beta^2 = \frac{v^2}{c^2} = 2\times 10^{-2}</math>
 
−
 
 
−
Plugging in the numbers:
 
−
: <math>\frac{dE}{dx} = \left ( 0.307 \frac{MeV \cdot cm^2}{g}\right ) (1)^2 (1) \frac{1}{2 \times10^{-2}} \ln \left( \frac{2 (1) (0.511 MeV) (2 \times10^{-2})}{21.8 eV} \frac{10^6 eV}{MeV}\right)</math>
 
−
: <math>= 105 \frac{MeV cm^2}{g}</math>
 
−
 
 
−
;How much energy is lost after 0.3 cm?
 
−
 
 
−
'''Notice that the units for energy loss are normalized by the density of the material'''
 
−
<math>\rho_{LH_2}</math> = 0.07 <math>\frac{g}{cm^3}</math><br>
 
−
 
 
−
To get the actual energy lost I need to multiply by the density.  So for any given atom the energy loss will depend on the state (solid, gas, liqid) of the atom as this effects the density of the material. 
 
−
 
 
−
:<math>\Delta E = (105 \frac{MeV cm^2}{g}) (0.07 \frac{g}{cm^3}) (0.3 cm)</math> = 2.2 MeV
 
−
 
 
−
[[Image:SPIM_HydrogenStoppingPower.pdf]] Compare with Triumf Kinematics Handbook, 2nd edition, September 1987, L.G. Greeniaus
 
−
 
 
−
====Bethe-Bloch Equation ====
 
−
 
 
−
While the classical equation above works in a limited kinematic regime, the Bethe-Bloch equation includes the corrections needed to cover most kinematic regimes for heavy particle energy loss.
 
−
 
 
−
:<math>-\frac{dE}{dx} = \mathcal{K} z^2 \frac{Z}{A} \frac{1}{\beta^2} \left [ \frac{1}{2} \ln \left (\frac{2 m_2 c^2 \beta^2 \gamma^2 T_{max}}{I^2} \right) - \beta^2 - \frac{\delta}{2}\right ]</math>[http://pdg.lbl.gov/2007/reviews/passagerpp.pdf PDG reference Eq 27.1 pg 1]
 
−
 
 
−
where
 
−
 
 
−
; <math>T_{max} = \frac{2 m_e c^2 \beta^2 \gamma^2}{1+ \frac{2 \gamma m_e}{M} + \frac{m_e}{M}}</math>
 
−
:= Max K.E. transferable to the Target of mass <math>M</math> in a single collision.
 
−
 
 
−
;<math>-\beta^2</math>
 
−
: =  correction for electron spin and very distant collisions which deforms the electron atomic orbits each process reducing dE/dx by <math>\frac{\beta}{2}</math>
 
−
 
 
−
; <math>\frac{\delta}{2}</math>
 
−
:= density correction term: in the classical derivation the material is treated as just a system of <math>N</math> atoms uniformly distributed in space.  These Atoms, however, give the material polarizability which can reduce the electric field (dielectric).
 
−
 
 
−
==== GEANT 4 implementation ====
 
−
 
 
−
The GEANT4 file  (version 4.8.p01)
 
−
 
 
−
source/processes/electromagnetic/standard/src/G4BetheBlockModel.cc
 
−
 
 
−
is used to calculate hadron energy loss.
 
−
 
 
−
line 132 <math>\Rightarrow</math>
 
−
 
 
−
:<math>-\frac{dE}{dx} =  \log \left ( \frac{2 m_e c^2 \tau (\tau +2) E_{min}}{I^2}\right) - \left (1 - \frac{E_{min}}{E_{max}} \right ) \beta^2</math>
 
−
 
 
−
where
 
−
 
 
−
:<math> \tau = \frac{K.E.}{M}</math>
 
−
 
 
−
line 143 <math>\Rightarrow</math>
 
−
 
 
−
:<math>\frac{dE}{dx} -= \log  ( \tau (\tau + 2) ) -cden</math> = density corection = <math>\frac{\delta}{2}</math>
 
−
 
 
−
line 148 <math>\Rightarrow</math>
 
−
 
 
−
:<math>\frac{dE}{dx} -= \frac{2c}{Z_{target}}</math> = shell correction, corrects for the classical asumption that the atomic electron's velocity is initially zero; or the the incident particles velocity is far greater than the atomic electron's velocity.
 
−
 
 
−
line 154 <math>\Rightarrow</math>
 
−
 
 
−
: <math>\frac{dE}{dx} *= \frac{2 \pi m_e c^2 r_e^2 z^2}{\beta^2} \rho_e \;\;\;\; \rho_e \propto \frac{NZ}{A}</math>
 
−
 
 
−
==== Energy Dependence ====
 
−
 
 
−
[[Image:SPIM_EnergyLoss_EnergyDependence.jpg | 500 px]]
 
−
 
 
−
The above curve shows the energy loss per distance traveled (<math>\frac{dE}{dx}</math>) as a function of the incident particles energy.  There are three basic regions.  At low incident energies ( < 10^5 eV) the incident particle tends to excite or even ionize the atoms in the material it is penetrating.  The maximum amount of energy loss per distance traveled is defined at as the Bragg peak.  The region after the Bragg peak in which the energy loss per distance traveled reaches its smallest value is refered to as the point of minimum ionizing.  Minimimum ionizing particles will have incident energies corresponding to this value or larger.  The characteristic of the minimum ionizing particles is that their energy loss per distance traveled is essentially constant making simulations easier until the particle's energy drops below the minimum ionizing energy level as it passes through the material.
 
−
 
 
−
In general the Bethe-Bloch equation breaks down at low energies (below the Bragg peak)  and is a good description (to within 10%) for
 
−
 
 
−
:<math>10 \frac{MeV}{a.m.u.} < E < 2 \frac{GeV}{a.m.u.}</math>  and <math>Z</math> < 26 (Iron)  : a.m.u = Atomic Mass Unit
 
−
 
 
−
the <math>\frac{1}{\beta^2}</math> term in the Bethe-Bloch equation dominates between the Bragg peak and the minimum ionization region.
 
−
 
 
−
the <math>\ln</math> term and its corrections influence the dependence of  <math>\frac{dE}{dx}</math> as you move up in energy beyond the minimum ionization point.
 
−
 
 
−
=== Energy Straggling ===
 
−
 
 
−
While the Bethe-Bloch formula gives you a way to quantify the amount of energy a heavy charged particle  loses as a function of the distance traveled, you should realize that when you calculate the total energy lost via
 
−
 
 
−
:<math> \Delta E = \int_{E_i}^{E_f} \left ( \frac{dE}{dx} \right ) dx</math>
 
−
 
 
−
you are only determining the AVERAGE energy loss.  In other words, Bethe-Bloch is the Astochastic process describing energy loss.
 
−
 
 
−
In reality the energy loss process is a stochastic process because of the statistical fluctuations which occur in the actual number of collisions which take place.
 
−
 
 
−
 
 
−
==== Thick Absorber ====
 
−
 
 
−
A thick absorber is one in which a large number of collisions takes place.  In this situation the central limit theorem from statistics tells you that the larger number of random variables <math>N</math> involved will result in observables which are distributed in a Gaussian manner. The Gaussian distribution is a good approximation to the binomial distribution when the number of trials is large.
 
−
 
 
−
[[Forest_ErrAna_StatDist#Binomial_with_Large_N_becomes_Gaussian]]
 
−
 
 
−
, and to a Poisson distribution when the mean is a lot larger than 1.
 
−
 
 
−
[[Forest_ErrAna_StatDist#Gaussian_approximation_to_Poisson_when]]
 
−
 
 
−
The gaussian probability function is defined as
 
−
 
 
−
: <math>P(x,\Delta) \propto e^{\frac{(\Delta - \bar{\Delta})^2}{ \sigma^2}}</math>
 
−
 
 
−
where the Full Width at Half Max (FWHM) of the distribution = <math>\left ( 2 \sqrt{2 \ln 2} \right ) \sigma</math>
 
−
 
 
−
In the case of energy loss, the variance using the Bethe-Bloch equation should be
 
−
 
 
−
:<math>\sigma_0^2 = 4 \pi N r_e^2 (m_e c^2)^2 \rho \frac{Z}{A} x</math>
 
−
 
 
−
the realitivistic variance is
 
−
 
 
−
: <math>\sigma^2 = [\frac{1-\beta^2/2}{1-\beta^2} ]\sigma_0^2</math>
 
−
 
 
−
for very thick absorbers see
 
−
 
 
−
C. Tschaler, NIM '''64''', (1968) 237 ; ''ibid'', '''61''', (1968) 141
 
−
 
 
−
When simulating energy loss of heavy charged particles the Bethe-Bloch equation may be used to calculate a <math>\frac{dE}{dx}</math> which can determine the average energy loss at the given kinetic energy of the particle.  This average is then smeared according to a gaussian distribution of variance
 
−
 
 
−
: <math>\sigma^2 =4 \pi N r_e^2 (m_e c^2)^2 \rho \frac{Z}{A} x [\frac{1-\beta^2/2}{1-\beta^2} ]</math>
 
−
 
 
−
====Thin Absorbers====
 
−
 
 
−
In thin absorbers the number of collisions is small preventing the use of the central limit theorem to describe the stochastic process of energy loss in terms of a Gaussian distribution.  The large energy transfers that are possible cause the energy loss distribution to look like a Gaussian  with a high energy tail (or foot).
 
−
 
 
−
The skewness of the resulting energy loss distribution is quantified as
 
−
 
 
−
:<math>\kappa = \frac{\bar{\Delta}}{W_{max}}</math>
 
−
 
 
−
: <math>\Delta  \equiv 2 \pi N_a r_e^2 m_e c^2 \rho \frac{Z}{A} \left ( \frac{z}{\beta}\right)^2 x </math> = lead term in Bethe Bloch equation
 
−
 
 
−
<math>\rho</math> = density of absorbing material.
 
−
 
 
−
:<math>W_{max} = \frac{(pc)^2}{\frac{1}{2} \left [  m_e c^2 + \left ( \frac{M^2 c^2}{m_e} \right ) \right ] + \sqrt{(pc)^2 + (Mc^2)^2}}</math> = max energy transfered in 1 collision (headon / knock out collision)
 
−
 
 
−
This comes from the relativistic kinematics of an Elastic Collision.<br>
 
−
[[Image:SPIM_ThinAbsorbers_Scatering.jpg | 400 px]]
 
−
 
 
−
: <math>\gamma = \frac{E_{tot}}{Mc^2} = \frac{ \sqrt{(pc)^2 + (Mc^2)^2}}{Mc^2}</math>
 
−
:<math>\beta= \frac{pc}{\gamma Mc^2} = \frac{pc}{E_{tot}}</math>
 
−
: <math>E_k = E_{tot} - Mc^2 = \gamma Mc^2 - Mc^2 = (\gamma - 1 ) Mc^2</math>
 
−
: <math>E_k = \sqrt{(pc)^2 + (Mc^2)^2} - Mc^2 </math>
 
−
:<math>  (p^{\prime}c)^2 = E_k^2 + 2E_km_ec^2</math>
 
−
 
 
−
 
 
−
Conservation of Momentum <math>\Rightarrow</math> :
 
−
 
 
−
: <math>\vec{p} = \vec{p}^{\; \prime \prime} + \vec{p}^{\; \prime}</math>
 
−
 
 
−
Conservation of Energy <math>\Rightarrow</math> :
 
−
 
 
−
: <math>E_{tot} + m_ec^2 = E_{tot}^{\prime \prime} + E_{tot}^{\prime}</math>
 
−
: <math>\sqrt{(pc)^2 + (Mc^2)^2} + m_ec^2 = \sqrt{(p^{\; \prime \prime} c)^2 + (Mc^2)^2}  + E_k +  m_e c^2</math>
 
−
 
 
−
 
 
−
using conservation of E & P as well as substituting for <math>p^{\prime}</math> you can show
 
−
 
 
−
:<math>(p^{\; \prime \prime}c)^2 = (pc)^2 - 2E_k\sqrt{(pc)^2 +(Mc^2)^2} + E_k^2</math> : cons of E
 
−
:<math>= (pc)^2 + E_k^2 + 2E_km_ec^2 -2pc\sqrt{E_k^2+2E_km_ec^2} \cos(\theta)</math> : cons of P
 
−
 
 
−
<math>\Rightarrow</math>
 
−
 
 
−
: <math>pc \cos(\theta) \sqrt{1+\frac{2m_ec^2}{E_k}} = \sqrt{(pc)^2+(Mc^2)^2} + m_ec^2</math>
 
−
 
 
−
solving for <math>E_k</math>
 
−
 
 
−
:<math>E_k = \frac{2m_ec^2(pc)^2\cos^2 (\theta)}{[\sqrt{(pc)^2 + (Mc^2)^2} +m_ec^2]^2 - (pc)^2 \cos^2 (\theta)}</math>
 
−
 
 
−
===== (Landau Theory) =====
 
−
<math>\kappa \leq 0.01</math>
 
−
 
 
−
Landau assumed
 
−
:# <math>W_{max} = \infty</math> is max energy transfer
 
−
:# electrons are free (energy transfer is so large you can neglect binding)
 
−
:# incident particle maintains velocity (large momentum transfer from big mass to small mass) (bowling ball hits ping pong ball)
 
−
 
 
−
L. Landau, "On the Energy Loss of Fast Particles by Ionization", J. Phys., vol 8 (1944), pg 201
 
−
 
 
−
instead of a gaussian distribution Landau used
 
−
 
 
−
: <math>P(x,\Delta) \propto \frac{1}{\bar{\Delta}\pi} \int_0^{\infty} e^{-u \ln u - u \lambda} \sin(\pi u) du</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>\lambda = \frac{1}{\bar{\Delta}} \left [ \Delta - \bar{\Delta} \ln \bar{\Delta} - \ln \epsilon + 1 -C \right ]</math>
 
−
: <math>\bar{\Delta} = 2\pi N_a r_e^2 m_e c^2 \rho \frac{Zz^2}{A \beta^2}x</math>
 
−
:<math>\ln \epsilon = \ln \left [ \frac{(1-\beta^2)I^2}{2m_ec^2 \beta^2} \right ]</math>
 
−
: <math>C = 0.577</math>
 
−
 
 
−
[[Image:SPIM_Landau_ThinAbsorberDist.jpg]]
 
−
 
 
−
===== (Vavilou's Theory) =====
 
−
 
 
−
Vavilous paper
 
−
 
 
−
P.V. Vavilou, "Ionization losses of High Energy Heavy Particles", Soviet Physics JETP, vol 5 (1950? )pg 749
 
−
 
 
−
describe the physics for the case
 
−
 
 
−
;<math>0.01 < \kappa < \infty </math>
 
−
 
 
−
The distribution function derived is shown below as well as a conceptual overlay of Vavilou's and Landau's distributions.  (The <math>\zeta f(x,\Delta)</math> in the picture should be a <math>\bar{\Delta}P(x,\Delta)</math> )
 
−
 
 
−
 
 
−
: <math>P(x,\Delta) = \frac{1}{\bar{\Delta}\pi} x e^{x(1+\beta^2C)} \int_0^{\infty} e^{xf_1} \cos(y \lambda_1 + xf_2) dy</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>f_1 = \beta^2 \left [ \ln(y) - C_i(y)\right ] - \cos(y) - y S_i(y)</math>
 
−
: <math>f_2 = y\left [ \ln(y) - C_i(y)\right ] + \sin(y) + \beta^2  S_i(y)</math>
 
−
: <math>C_i(y) \equiv - \int_y^{\infty} \frac{\cos(t)}{t} dt</math>
 
−
: <math>S_i(y) \equiv \int_0^{y} \frac{\sin(t)}{t} dt</math>
 
−
: <math>C = 0.577</math>
 
−
 
 
−
[[Image:SPIM_Vavilou_Landau_ThinAbsorber.jpg | 400 px]]
 
−
 
 
−
==== GEANT4's implementation ====
 
−
 
 
−
GEANT 4 uses the skewness parameter <math>\kappa</math> to determine if it will use a "fluctuations model" to calculate energy straggling or the gaussian model described in section 3.2.1.
 
−
 
 
−
===== kappa > 10 =====
 
−
If
 
−
: <math>\kappa \equiv \frac{ \bar{\Delta}}{W_{max}}</math> > 10
 
−
and we have a thick absorber ( large step size) then the Gausian function in 3.2.1 is used to calculate energy straggling.
 
−
 
 
−
What happens is <math>\Delta E</math> is calculated via <math>\int_{E_i}^{E_f} \frac{dE}{dx} dx</math> then the actual energy loss predicted by the simulation is chosen from a Gaussian distribution to account for energy straggling such that the <math>\sigma</math> of this Gaussian distribution is given by:
 
−
 
 
−
:<math>\sigma^2 = 2 \pi r_e^2m_ec^2N_{el} \frac{Z_h}{\beta^2} T_C s (1 - \frac{\beta^2}{2})</math>
 
−
 
 
−
where
 
−
 
 
−
:<math>N_{el}</math> = electron density of the medium
 
−
:<math>Z_h</math> = charge of the incident particle
 
−
:<math>s</math> = step size
 
−
:<math> T_C</math> = cutoff kinetic energy for <math>\delta </math>-electrons
 
−
 
 
−
<math>T_C</math> tells GEANT where to put the cutoff for using the Gaussian distribution for energy straggling.  This tells the simulation the low energy cutoff where Bethe-Bloch starts to fail due to ionization.
 
−
 
 
−
=====Delta-electrons =====
 
−
What is a <math>\delta</math> - electron?
 
−
 
 
−
<math>\delta</math> - electrons are also known as "knock -on" electrons and delta rays.
 
−
 
 
−
As heavy particles traverse a medium they can ionize electrons from atoms.  The ejected electrons can be given enough energy to ionize as well.
 
−
 
 
−
In a cloud chamber (a supercooled volume of super saturated water vapor which ionizes as charged particles pass through)  such and event would look like:
 
−
 
 
−
[[Image:SPIM_DeltaRay_CloudChamber.jpg | 400 px]]
 
−
 
 
−
The blue spiral in the above gas chamber picture is a high energy electron from the collision.  The B-field is directed out of the picture.
 
−
 
 
−
The physics of ionization is different from the physics used to calculate Bethe-Bloch energy loss.  Remember Bethe-Bloch  starts to break down at low energies below the Bragg peak. 
 
−
 
 
−
Because of this GEANT 4 sets the cutoff for this process to be
 
−
 
 
−
: <math>T_{cut}</math> > 1 keV
 
−
 
 
−
 
 
−
Note:  The BE energies of an electron in Hydrogen is 13.6 ev and the electrons in Argon have binding energies between 15.7 eV and 3.2 keV.
 
−
 
 
−
===== Fluctuations Model: kappa < 10=====
 
−
 
 
−
If <math>\kappa \equiv \frac{ \bar{\Delta}}{W_{max}} < \frac{\Delta E}{T_C}</math>
 
−
 
 
−
Then GEANT 4 uses a "Fluctuations Model" to determine energy loss instead of Bethe-Bloch.
 
−
 
 
−
; Fluctuations Model
 
−
:# the atom is assumed to have on 2 energy levels <math>E_1</math> and <math>E_2</math>
 
−
:# you can excite the atom and lose either <math>E_1</math> or <math>E_2</math> or you can ionize the atom and lose energy according to a <math>\frac{1}{E^2}</math> function <math>u_j</math>.
 
−
 
 
−
The total energy loss in a step will be
 
−
 
 
−
: <math>\Delta E = \Delta E_{exc} + \Delta E_{ion}</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>\Delta E_{exc} = \eta_1 E_1 + \eta_2 E_2</math>
 
−
 
 
−
: <math>\Delta E_{ion} = \sum_{j=1}^{\eta_3} \frac{I}{1 - u_j \frac{T_{up}-I}{T_{up}}}</math>
 
−
 
 
−
:<math>\eta_1</math>, <math>\eta_2</math>, and <math>\eta_3</math> are the number of collisions which are sampled from a poison distribution
 
−
 
 
−
:<math>u_j = \int_{I}^{E_j} \frac{I T_{up}}{T_{up} - I} \frac{dx}{x^2}</math>
 
−
 
 
−
: <math>E_j = \frac{I}{1- rand  \frac{T_{up}-1}{T_{up}}}</math> : rand = random number between 0 and 1
 
−
 
 
−
:<math>T_{up} = \left \{  {~ 1 keV \;  threshold \;energy \;for \; \delta- ray \; production \atop T_{max} \; \;\;\; if \; T_{max} < 1 keV} \right .</math>
 
−
 
 
−
: <math>I</math> = mean ionization energy
 
−
 
 
−
:<math>E_2 \approx (10 eV) Z^2</math>
 
−
 
 
−
:<math>\ln E_1 = \frac{\ln (I) - f_2 \ln (E_2)}{f_1}</math>
 
−
 
 
−
:<math>f_1 + f_2 =1</math>
 
−
 
 
−
:<math>f_2 =\left \{  {0 \; z=1 \atop \frac{2}{z} \; z \ge 2} \right .</math>
 
−
 
 
−
The fluctuation model was comparted with data in
 
−
 
 
−
K. Lassila-Perini and L. Urban, NIM, A362 (1995) pg 416
 
−
 
 
−
The cross sections used for excitation and ionization may be found in
 
−
 
 
−
H. Bichel, Rev. Mod. Phys., vol 60 (1988) pg 663
 
−
 
 
−
=== Range Straggling===
 
−
 
 
−
;Def of Range (R):
 
−
: The distance traveled before all the particles energy is lost.
 
−
 
 
−
:<math>R \equiv \int_0^T \frac{dE}{\frac{dE}{dx}}</math>
 
−
:  = theoretical calculation of the path length traveled by a particle of incident energy <math>T</math>
 
−
 
 
−
: Note units: <math>\left [ R \right ] =  \frac{g}{cm^2} ; \left [ \frac{dE}{dx} \right ] = \frac{MeV \cdot cm^2}{g}</math>
 
−
 
 
−
the Energy Straggling introduced in the previous section can explain why identical particles penetrate material to different depths.  The energy straggling results in Range straggling.
 
−
 
 
−
If we do a shielding experiment where we have a source of incident particles of energy E and we count how many "punch" through a material of thickness (x) we would see a transmission coefficient <math>\left ( \frac{N_{out}}{N_{in}} \right) </math> which would look like
 
−
 
 
−
[[Image:SPIM_RangeStraggling.jpg | 400 px]]
 
−
 
 
−
====Fractional Range Straggling ====
 
−
 
 
−
<math>\frac{\sigma_R}{R} \equiv</math> fractional range straggling
 
−
 
 
−
Assuming the energy loss of a non-relativistic heavy ion through matter follows a Gaussian (thick absorber)
 
−
 
 
−
Then it can be shown that
 
−
 
 
−
<math>\frac{\sigma_R}{R} \approx \frac{1}{2} \sqrt{\frac{M}{A}}</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>M</math> = mass of the target electrons
 
−
 
 
−
: <math>A</math> = atomic mass of the Projectile
 
−
 
 
−
since
 
−
 
 
−
:<math>m_e = 9.11 \times 10^{-31}</math> kg
 
−
 
 
−
and
 
−
 
 
−
: 1 a.m.u. = <math>1.66 \times 10^{-27}</math> kg
 
−
 
 
−
then
 
−
 
 
−
: <math>\frac{\sigma_R}{R} \approx \frac{1}{2} \sqrt{\frac{9.11 \times 10^{-31}}{1.66 \times 10^{-27}A}}</math>
 
−
: = 1.17 % if using a proton (A=1)
 
−
 
 
−
The above is a "back of the envelope" estimate.  The experimentally measured values for Cu, Al, and Be target using a proton projectile are
 
−
 
 
−
[[Image:SPIM_RangeStrag_SigmaR_overR.jpg | 400 px]]
 
−
 
 
−
If the incident projectile is an electron then <math>\frac{\sigma_R}{R}  \approx \frac{1}{2}</math> making electron range straggling a vague concept.
 
−
 
 
−
There are several definitions of electron range
 
−
 
 
−
;1.) Maximum Range (<math>R_0</math>):
 
−
:This range is defined using the continuous slowing down approximation (CSDA) in the electrons are assumed to have many collisions over very small dimensions making it appear to be continuous energy loss instead of discrete.  The range is then calculated by integrating over these average energy losses <math>\frac{dE}{dx} \cdot s</math>.
 
−
 
 
−
;2.) Practical Range (<math>R_P</math>):
 
−
: This stopping distance is defined by extrapolating the electron transmission curve to zero (see below).
 
−
 
 
−
[[Image:SPIM_PracticalRangStraggline_4Electrons.jpg | 400 px]]
 
−
 
 
−
=== Electron Capture and Loss ===
 
−
====Bohr Criterion====
 
−
:"A rapidly moving nucleus is fully ionized if its velocity exceeds that of its most tightly bound electron"
 
−
 
 
−
The Bohr Model:
 
−
:<math>\Rightarrow E = \frac{mz^2e^4}{8 \epsilon_0^2 h^2 n^2}</math>
 
−
 
 
−
for the inner most electron (<math>n=1</math>)
 
−
 
 
−
:Electron K.E. = <math>\frac{1}{2} mv^2 = \frac{mz^2e^4}{2(4\pi \epsilon_0)^2 \hbar^2} \Rightarrow v = \frac{z e^2}{4 \pi \epsilon_0 \hbar}</math>
 
−
 
 
−
 
 
−
:the fine structure constant <math>\alpha \equiv \frac{e^2}{4 \pi \epsilon_0 \hbar c} = \frac{1}{137}</math>
 
−
 
 
−
:<math> v = zc \alpha</math>
 
−
 
 
−
If <math>v > zc \alpha</math> the nucleus is fully ionized
 
−
 
 
−
or
 
−
 
 
−
if <math>\frac{z}{v/c} = \frac{z}{\beta} < \frac{1}{\alpha} = 137</math>
 
−
 
 
−
alternatively if the ion is moving through a material with a speed such that
 
−
 
 
−
:<math>\frac{z}{\beta} > \frac{1}{\alpha} =137</math>
 
−
 
 
−
 
 
−
Then electrons may be captured by the projectile and lost by the target.
 
−
 
 
−
==== Z-effective====
 
−
Describing the charge state of your heavy ion traveling through matter at a velocity below the Bohr criterion is very complicated.  There is a competition between electron capture and loss.  Accurate cross sections are needed to simulate the process reliably.
 
−
 
 
−
Some insight into this process can be found using the Thomas-Fermi model
 
−
 
 
−
: <math>V \propto \frac{Ze^{-r/a}}{r}</math>
 
−
 
 
−
to describe an atom moving slow enough so it has captured many electrons but fast enough so its not neutral.  In the Thomas-Fermi model the distribution of electrons in an atom is described as being uniformly distributed such that there are 2 electrons in each discrete volume of phase space( the space in which all possible states of a system are represented)  defined using planks constant as <math>h^3</math>. 
 
−
 
 
−
For the purpose of simulations you would like a relationship for <math>Z_{eff}</math> in terms of <math>\beta</math> and <math>Z</math>. 
 
−
 
 
−
It is usually adequate to use fits for empirical data as long as we know that we are in the kinematic range in which those fits are valid.
 
−
 
 
−
 
 
−
when <math>E < 10</math> MeV the data indicates that
 
−
 
 
−
: <math>Z_{eff} = Z(1 - e^{-\beta\frac{B}{Z^{2/3}}})</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>B = 130 \pm 5</math>
 
−
:<math>Z_{eff} \equiv</math> effective charge f the projectile = <math>Z - \bar{q}_c</math>
 
−
: <math>Z</math> = number of protons
 
−
:<math>\bar{q}_c</math> = average number of captured electrons
 
−
 
 
−
 
 
−
 
 
−
'''When calculating stopping power for E < 10 MeV you use <math>Z_{eff}</math> in the Bethe-Bloch equation.'''
 
−
 
 
−
Note:  As the ions charge state fluctuates while it slows down (or if accelerated through materials) you will need to recalculate the energy loss, and as a result you will get larger energy loss fluctuations in this energy range.
 
−
 
 
−
For thin absorber you will look for stripping and loss cross sections.
 
−
 
 
−
: Here a thin absorber is one whose thickness is less than the charge equilibrium distance defined as the distance traveled until the projectile's velocity is <math> v \ll zc\alpha</math>
 
−
 
 
−
A rule of thumb is that a thin absorber for low energy ions has a thickness <math>\le \frac{5 \frac{\mu g}{cm^2}}{\rho}</math>
 
−
 
 
−
For thick absorbers:  The experimentally determined expression for the change in <math>Z_{eff}</math> from <math>Z</math> is
 
−
 
 
−
:<math>\Delta Z_{eff} = \frac{1}{2} \sqrt{ \left [  Z_{eff} \left (1 - \frac{Z_{eff}}{Z} \right )^{1.67}\right ] }</math>
 
−
 
 
−
=== Multiple Scattering ===
 
−
 
 
−
The Bethe-Bloch equation tells us how much energy is lost and GEANT4s calculation of this energy is described above.
 
−
 
 
−
Now we need to know which direction the scattered particle goes after it has lost this energy.
 
−
 
 
−
The work of Moliere describes the angular deflection of the particle which lost the energy thereby leading to a prediction of the Cross-section.  GEANT4 though uses the more complete Lewis theory to describe Multiple Couloumb Scattering (MCS) sometimes generically referred to as multiple scattering.
 
−
 
 
−
There are 3 regions in which coulomb scattering is calculated
 
−
 
 
−
; 1.) Single Scattering:
 
−
: For thin materials.
 
−
: If the probability of more than 1 coulomb scattering is small
 
−
:The use the Rutherford formula for <math>\frac{d \sigma}{d \Omega}</math>
 
−
 
 
−
;2.)Multiple Scattering:
 
−
: In this case the number of independent scatterings is large (N > 20) and the energy loss is small such that the problem can be treated statisticaly to obtain a probability distribution for the net deflection angle <math> [P(\theta)]</math> as a function of the material thickness that is traversed.
 
−
 
 
−
 
 
−
;3.) Plural Scattering:
 
−
: If 1< N <math>\le</math> 20 then you can't use Rutherford to describe the scattering nor use a normal random statistical description.
 
−
 
 
−
see E. Keil, Z. Naturforsch, vol 15 (1960), pg 1031
 
−
 
 
−
 
 
−
Reviews of rigorous multiple scattering calculations may be found in
 
−
: P.C. Hemmer, et. al., Phys. Rev, vol 168 (1968), pg 294
 
−
 
 
−
==== GEANT4's implementation of MSC (N>20) ====
 
−
 
 
−
GEANT4 models MSC when N>20 using model functions to determine the angular and spatial distributions chosen to give the same moments of these distributions as the Lewis theory.
 
−
 
 
−
:H.W. Lewis, Phys. Rev., vol 78 (1950), pg 526
 
−
 
 
−
modern versions of the above are at
 
−
 
 
−
: J.M. Fernandez-Varea, et. al., NIM, B73 (1993), pg 447
 
−
: I. Kawrakow, et. al., NIM, B142 (1998) pg 253
 
−
 
 
−
When N>20 multiple scattering can be described as a statistical process using a modified version of the Boltzman transport equation from statistical mechanics. 
 
−
 
 
−
;Note: The simulation step size is chosen such that (N>20),  If you have materials so thin that N < 20 then GEANT4 will likely skip the material.  (one way around this is to increase the thickness and change the density).  If the material thickness can't be increased because its sandwhiched between two other materials then you will need to write a special step algorithm for the volume and have GEANT4 use it for the step.
 
−
 
 
−
 
 
−
Let <math>f(s,\vec{x},\hat{v}) \equiv</math> the distribution function for a system of incident particles traveling through a material.
 
−
 
 
−
where
 
−
 
 
−
:<math>s =</math> arc length of the particle's path through the material
 
−
:<math>\vec{x} =</math> position of a charged particle
 
−
: <math>\hat{v} =</math> direction of motion of the particle <math>\frac{\vec{v}}{|\vec{v}|}</math>
 
−
 
 
−
The multiple scattering experienced by a single charged particle traveling through the material is then simulated by sampling from the distribution <math>f(s,\vec{x},\hat{v} )</math>
 
−
 
 
−
The governing transport/diffusion equation is based on the continuity equation but with a "sink" term representing the possibility of collisions ejecting particles out of the volume.
 
−
 
 
−
[[Image:SPIM_MultScatDiffEq.jpg | 400 px]]
 
−
 
 
−
:<math>\frac{\partial f(s,\vec{x},\hat{v} ) }{\partial s} + \hat{v} \bullet \vec{\nabla}f(s,\vec{x},\hat{v} ) = N \int \sigma(\hat{v} \bullet\hat{v}^{\prime} )\left [ f(s,\vec{x},\hat{v}^{\prime} )  - f(s,\vec{x},\hat{v} ) \right ] d \hat{v}^{\prime}</math>
 
−
 
 
−
where
 
−
 
 
−
:<math>N</math> = number of atoms per volume
 
−
:<math>\sigma(\hat{v} \bullet\hat{v}^{\prime} )</math> = cross sections for elastic scattering per Solid angle <math>\left ( \frac{d \sigma}{d \Omega} \right )</math>
 
−
 
 
−
To solve the above diffusion equation the distribution function, <math>f(s,\vec{x},\hat{v} )</math> is expanded in Spherical Harmonics ( <math>Y_{\ell}^m(\theta,\phi)</math> ) and <math>\sigma</math> expand  in Legendre Polynomials (<math>P_N(cos \theta)</math>) since it has no <math>\phi</math> angle dependence.
 
−
 
 
−
;Note: For Coulomb Scattering in polar coordinates you can write the potential in terms of Legendre Polynomials such that:
 
−
 
 
−
:<math>U=k \frac{q}{r}</math>
 
−
:= <math>k\frac{q}{\sqrt{r^2-a^2-2ar \cos \theta}}</math>  in polar coordinates
 
−
:= <math>k\frac{q}{r} \sum_{n=0}^{\infty} P_n(\cos \theta) \left ( \frac{a}{r}\right )^n</math> (the sqrt term above is expanded using binomial series
 
−
 
 
−
:  <math>f(s,\vec{x},\hat{v} ) = \sum_{\ell,m} f_{\ell,m}(\vec{x},s) Y_{\ell}^m(\hat{v})</math>
 
−
 
 
−
after substituting into the diffusion equation and doing the integral on the righ hand side you get
 
−
 
 
−
:<math>\frac{\partial f_{\ell,m}(\vec{x},s) }{\partial s} + \frac{f_{\ell,m}(s,\vec{x},\hat{v} }{\lambda_{\ell}} = - \sum_{\lambda, i, j} \vec{\nabla} f_{i,j}(\vec{x},s ) \bullet \int Y_{\ell,m}^{\star} \hat{v} Y_{i,j} d \hat{v} \; \; \; \; \; \; \; \;\hat{v} = f(\theta.\phi)</math>
 
−
 
 
−
where
 
−
: <math>\frac{1}{\lambda_{\ell}} = 2 \pi N \int_0^{\pi} \left [ 1-P_{\ell}(\cos \theta)\right ] \sigma(\theta) \sin(\theta) d \theta</math> = <math>\ell^{th}</math> transport mean free path for the <math>f_{\ell}</math> distribution function ( <math>\phi</math> symmetry is assumed making it <math>m</math> independent)
 
−
 
 
−
From the above one can find the average distances traveled and the average deflection angle of the distribution.  Again, see :
 
−
 
 
−
: J.M. Fernandez-Varea, et. al., NIM, B73 (1993), pg 447
 
−
 
 
−
 
 
−
The "moments" of <math>f(s,\vec{x},\hat{v}) </math>  are defined as
 
−
 
 
−
:<math><z> = 2 \pi \int z f(s,\vec{x},\hat{v}) \sin(\theta) d \theta d |\vec{x}| = \lambda_1 \left [ 1-e^{-s/\lambda_1}\right ]</math> = mean geometrical path length
 
−
:<math><\cos(\theta)> = 2 \pi \int_{-1}^1 \sum_{\ell} P_{\ell}(\cos \theta) \int f(s,\vec{x},\hat{v}) \sin(\theta) d \theta d |\vec{x}| = e^{-s/\lambda_1}
 
−
</math>
 
−
:<math>\frac{1}{\lambda_1} = 2 \pi N \int_0^{\pi} \left [ 1-P_1(\cos \theta)\right ] \sigma(\theta) \sin(\theta) d \theta</math>
 
−
 
 
−
Notice there are 3 lengths
 
−
 
 
−
[[Image:SPIM_MultScatDiffEq_PathLength.jpg | 400 px]]
 
−
 
 
−
:<math>s</math> = geometrical path length between endpoints of the step =<math> \left \{  {line \; if \; \vec{B} = 0 \atop arc \; if \; \vec{B} \ne 0 } \right .</math>
 
−
:<math>t</math> = true path length = actual length of the path taken by particle
 
−
:<math><z></math> - mean geometrical path length along the z-axis
 
−
 
 
−
In GEANT4 the <math>\lambda_{\ell}</math>'s are taken from
 
−
 
 
−
If 100 eV < K.E. of electron or positron < 10 MeV
 
−
 
 
−
:D. Liljequist, J. Applied Phys, vol 62 (1987), 342
 
−
:J. Applied Phys, vol 68 (1990), 3061
 
−
 
 
−
If K.E. > 10 MeV
 
−
 
 
−
:R. Mogol, Atomic Data, Nucl, Data tables, vol 65 (1997) pg 55
 
−
 
 
−
 
 
−
 
 
−
with <z> now known GEANT will try to determine "<math>t</math>" for the energy loss and scattering calculations.
 
−
 
 
−
A model is used for this where
 
−
 
 
−
:<math>t=\frac{1}{\alpha} \left [ 1 - (1- \alpha \omega z)^{\frac{1}{\omega}})\right ]</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>\omega = 1 + \frac{1}{\alpha \lambda_{10}}</math>
 
−
: <math>\alpha =\left \{  {\frac{\lambda_{10} - \lambda_{11}}{s \lambda_{10}}\;\;\;\; K.E. \ge M_{particle} \atop \frac{1}{R}\;\;\;\; K.E. < M_{particle}} \right .</math>
 
−
: <math>s</math> = stepsize
 
−
: <math>\lambda_{10} - \frac{\lambda_1}{1-\alpha s}</math>
 
−
:<math>\lambda_{11} = \lambda_1</math> at end of strep
 
−
 
 
−
while <math><cos \theta ></math> is calculable, GEANT4 evaluates <math>\cos (\theta)</math> from a probability distribution whose general form is
 
−
 
 
−
:<math>g[\cos(\theta)] - p \left ( qg_1[\cos(\theta)] + (1=q)g_3[\cos(\theta)] \right ) + (1-p)g_2[\cos(\theta)]</math>
 
−
 
 
−
where
 
−
 
 
−
:<math>g_1(x) = C1e^{-a(1-x)}</math>
 
−
:<math>g_2(x) = \frac{C_2}{(b-x)^d}</math>
 
−
:<math>g_3(x) = C_3</math>
 
−
 
 
−
:<math>C_1, C_2, C_3</math> are normalization constants
 
−
:<math>p,q,a,b,d</math> are parameters which follow the work reported in
 
−
 
 
−
:V.L. Highland, NIM, vol 219 (1975) pg497
 
−
 
 
−
The GEANT4 files in version 4.8 were located in
 
−
 
 
−
/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc
 
  
−
and
+
/control/shell mv currentEvent.rndm currentEvent10.rndm
  
−
/source/processes/electromagnetic/standard/src/G4MscModel.cc
 
  
−
/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc
+
You can restore the random number generator and begin generating random number from the last save time
  
−
== Interactions of Electrons and Photons with Matter==
+
/random/resetEngineFrom currentEvent.rndm
−
=== 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
+
==Building GEANT4.11==
  
−
;Note: Bethe & Heitler first calculated this radiation in 1934 which is why you will sometimes hear Bremsstrahlung radiation refererd to as Bethe-Heitler.
+
===4.11.2===
 +
[[TF_GEANT4.11]]
  
−
:<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>
+
==Building GEANT4.10==
  
−
where
 
  
−
: <math>E_0</math> = initial total energy of the electron
+
===4.10.02===
−
:<math>E</math> = final total energy of the electron
 
−
: <math>\nu = \frac{E-E_0}{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
 
−
:<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>
+
[[TF_GEANT4.10.2]]
 +
===4.10.01===
  
−
:<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>
+
[[TF_GEANT4.10.1]]
−
:<math>\phi_2(\gamma) = \phi_1(\gamma) - \frac{2}{3}(1+6.5 \gamma + 6 \gamma^2)</math>
 
  
 +
==Building GEANT4.9.6==
  
−
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]
+
[[TF_GEANT4.9.6]]
  
−
:if <math>3 \ge Z < 5</math> use Equation 3.46 and 3.47
+
==Building GEANT4.9.5==
  
−
:if <math> Z < 2</math> use Equation 3.25 and 3.26
+
[[TF_GEANT4.9.5]]
  
−
;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)
+
An old version of Installation notes for versions prior to 9.5
  
−
;We now know 2 ways charged particles can loose energy when passsing through matter.
+
[http://brems.iac.isu.edu/~tforest/NucSim/Day3/ Old Install Notes]
  
−
;Energy loss: <math>\left ( \frac{dE}{dx} \right )_{tot} = \left ( \frac{dE}{dx} \right )_{rad} + \left ( \frac{dE}{dx} \right )_{col}</math>
 
−
:<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>
+
Visualization Libraries:
  
−
where
+
[http://www.opengl.org/ OpenGL]
  
−
: <math>N = \frac{ number\; atoms}{cm^3} = \frac{\rho N_a}{A} = \frac{density \;\times\; Avagadros\;\;Number}{Atomic number}</math>
+
[http://geant4.kek.jp/~tanaka/DAWN/About_DAWN.html  DAWN]
−
:<math>\left ( h \nu \right )</math> = Energy of emmitted photon
 
−
:<math>\left ( \frac{d \sigma}{d \nu} \right )</math> = Probabitlity of Energy loss
 
  
  
−
The quantity<math> \Phi_{rad} </math> is defined such that
+
[http://doc.coin3d.org/Coin/  Coin3D]
  
−
:<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>
+
==Compiling G4 with ROOT==
  
−
<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>
+
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>\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>
+
[[G4CompileWRootforTracks]]
  
−
where
+
==Using SLURM==
  
−
: <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>
+
http://slurm.schedmd.com/quickstart.html
−
: <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
 
  
−
: <math>- \frac{dE}{dx} = N E_0 \Phi_{rad}</math>
+
https://rc.fas.harvard.edu/resources/documentation/convenient-slurm-commands/
  
−
;Note: for intermediate value of <math>\gamma</math> you need to integrate numerically
+
===simple batch script for one process job===
−
:<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. 
+
create the file submit.sbatch below
  
−
[[Image:SPIM_Bethe-Brem_Eloss-vs-Energy.jpg | 600 px]]
+
<pre>
 +
#!/bin/sh
 +
#SBATCH --time=1
 +
cd src/PI
 +
./PI_MC 100000000000000
 +
</pre>
  
 +
the execute
  
−
;Critical Energy<math> E_C</math>:
+
:sbatch submit.sbatch
−
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
+
check if its running with
  
−
:<math>E_C \sim \frac{800 MeV}{Z+1.2}</math>
+
:squeue
  
−
;Examples:
+
to kill a batch job
−
{| 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 ====
+
:scancel JOBID
−
; 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]]
+
===On minerve===
  
−
as a result
+
Sample script to submit 10 batch jobs.
  
−
<math>d \sigma_{tot} = \frac{Z(Z+1)}{Z^2} d \sigma_{Brem}</math>
+
the filename is minervesubmit and you run like
−
: = <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>
+
  source minervesubmit
  
−
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).
+
<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>
  
 +
The file submit10mil looks like this
 +
<pre>
 +
#!/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>
  
−
==== Coherent Bremsstrahlung ====
 
  
 +
use
  
−
The above equations represent the physics of incoherent bremsstrahlung production.
+
qstat
  
 +
to check that the process is still running
  
−
http://137.99.79.133/halld/tagger/references/inelasticBrems-94.pdf
+
use
  
−
==== Radiation Length (Xo)====
+
qdel jobID#
  
−
;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>.
+
if you want to kill the batch job, the jobID number shows up when you do stat.
  
−
in the high energy limit where <math>\left ( \frac{dE}{dx} \right )_{col}</math> can be ignored <math>( E > E_C )</math>
+
for example
−
 
+
<pre>
−
:<math>\frac{dE}{dx} = NE_0 \Phi_{rad} \Rightarrow \int_{E_0}^{E} \frac{dE}{E} = -\int_0^X N \Phi_{rad} dx</math>
+
[foretony@minerve HW10]$ qstat
−
: <math>\Rightarrow \ln(\frac{E}{E_0}) =-N \Phi_{rad}X</math>
+
Job id                    Name            User            Time Use S Queue
−
 
+
------------------------- ---------------- --------------- -------- - -----
−
or
+
27033.minerve            submit          foretony        00:41:55 R default       
−
 
+
[foretony@minerve HW10]$ qdel 27033
−
: <math>E = E_0 e^{-N \Phi_{rad}X} = E_0 e^{-\frac{X}{X_0}}</math>
+
[foretony@minerve HW10]$ qstat
−
 
+
</pre>
−
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"
 
−
|colspan="2"|
 
−
 
 
−
====Table of Radiation Lengths for several materials====
 
−
|-
 
−
|-
 
−
|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>
 
−
: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
 
−
 
 
−
: <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
 
−
:<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
 
−
: <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
 
−
 
 
−
;<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
 
−
 
 
−
;<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
 
−
 
 
−
;<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
 
−
 
 
−
; After 2.3 radiation lengths the electron energy is down by a factor of 10 from its original value.
 
−
 
 
−
====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.
 
−
 
 
−
:<math>T_C</math> = incident particle K.E. cutof = secondary particle production threshold
 
−
:<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
 
−
:<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)
 
−
: reference: J. Tuli, "Evaluated Nuclear Structure Data File", BNL-NCS - 51655 -Rev 87, 1987 from Brookhaven Nat. Lab
 
−
: see  [http://www.nndc.bnl.gov National Nuclear Data Center]:
 
−
 
 
−
; 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.
 
−
 
 
−
To improve simulation speed though, GEANT 4 actually uses a fit to the above cross sections such that
 
−
 
 
−
:<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>
 
−
 
 
−
where
 
−
 
 
−
: <math>m</math> = electron mass
 
−
:<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
 
−
 
 
−
: The energy of the emitted photon is determine by sampling a probability distribution from
 
−
:S. Seltzer and M. Berger, Atomic Data & Nucle. Data tables, vol 35 (1986) pg 345-418
 
−
 
 
−
and
 
−
 
 
−
: the angular distribution (<math>\cos(\theta)</math>)  is sampled according to
 
−
:E. Acosta, Appl. Phys. Letter, vol 80 # 17 (2002) pg 3228-3230
 
−
 
 
−
 
 
−
Note
 
−
:* The MC program PENELOPE was used to generate the energy distributions that are sampled
 
−
:* GEANT4 uses a modified version of  base equations for <math>e^+e^-</math> bremsstrahlung with model corrections for <math>e^+</math>
 
−
 
 
−
; 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 ====
 
−
 
 
−
The radiation length for compounds and mixtures is determined by parallel weighting (resistors in parallel)
 
−
 
 
−
:<math>\frac{1}{L_{rad}} = \omega_1 \left ( \frac{1}{L_{rad}} \right )_1 + \omega_2 \left ( \frac{1}{L_{rad}} \right )_2</math>
 
−
 
 
−
where
 
−
 
 
−
: <math>\omega_1</math> = fraction , by wieght, of each element in the mixture/compound.
 
−
: <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===
 
−
 
 
−
The photo-electric effect identifies the physics process by which bound electrons in an atom are liberated by an interaction with an incident photon.
 
−
 
 
−
[[Image:SPIM_PhotoElectricEffectSchematic.jpg | 400 px]]
 
−
 
 
−
:<math>E_f = E-E_B = h \nu - E_B</math>
 
−
 
 
−
where
 
−
 
 
−
:<math>h\nu</math> = incident photon energy
 
−
:<math>E_B</math> = electron binding energy
 
−
 
 
−
==== Moseley's Law ====
 
−
 
 
−
Moseley's law approximates the binding energies of electrons in atoms as
 
−
 
 
−
: <math>E_B = 13.605 \frac{\left ( Z-k_s \right )^2 }{n^2}</math> (eV)
 
−
 
 
−
X-ray electron shells are labeled K,L,M
 
−
 
 
−
 
 
−
{| 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 produced 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.
 
−
 
 
−
=== 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} = \frac{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 wihtin 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 =====
 
−
; Triples 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 be found 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}{\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>
 
−
 
 
−
[[File:Bigg_Lighthill_conversion_1Mev-100MeV_SandiaPreprint_SC-RR-68-619.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====
 
−
 
 
−
== Homework Problems==
 
−
[[HomeWork_Simulations_of_Particle_Interactions_with_Matter]]
 
−
 
 
−
== Changing the Random number seed in GEANT4 ==
 
−
 
 
−
You can set the random number seed by adding the relevant calls in the vis.mac file
 
−
 
 
−
typing the command below turns on the ability to save the random number seed
 
−
 
 
−
/random/setSavingFlag 1
 
−
 
 
−
then create a file to save the random seed to
 
−
 
 
−
/random/saveThisRun
 
−
 
 
−
is will create the files
 
−
 
 
−
currentEvent.rndm
 
−
currentRun.rndm
 
−
 
 
−
 
 
−
then use the command
 
−
 
 
−
/random/resetEngineFrom  currentRun.rndm
 
−
 
 
−
to restart the run using the random number seed stored in the file "currentRun.rndm"
 
−
 
 
−
The easiest way to use this is to create several of these random number seed files by running geant for X events and saving the seed to a different filename each time.
 
−
 
 
−
==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]]
 
  
 +
==Definitions of Materials==
  
−
===Running on Inca===
+
[[File:MCNP_Compendium_of_Material_Composition.pdf]]
  
−
To use my version of software for this class try
+
==Minerve2 GEANT 4.10.1 Xterm error==
  
−
;ROOT
 
−
:setenv ROOTSYS /home/tforest/src/ROOT/root
 
−
:$ROOTSYS/bin/root
 
  
−
;GEANT4
+
On OS X El Capitan V 10.11.4 using XQuartz
−
:source /home/tforest/src/GEANT4/geant/env.csh
 
  
−
=Assessment=
+
~/src/GEANT4/geant4.10.1/Simulations/B2/B2a/exsmpleB2a
−
==learning outcomes==
 
  
−
#Student will learn how to compile and run the GEANT4 software.
+
<pre>
−
#Students will learn how the alter target materials and incident projectile kinematics
 
−
#Student will learn how to extract trajectory information for analysis
 
  
−
==assess if outcomes are meet==
+
# 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>
  
−
A series of simulation projects will be completed during the course to assess the students ability to meet outcome goals
 
  
−
==Shortfall Identification==
 
  
−
If outcomes are not met the shortfall will be identified by comparing available learning materials to students goal achievement level.
 
  
−
==Remediation==
 
  
−
Projects will be adjusted to overcome identified shortfall.
+
[[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

Retrieved from "https://wiki.iac.isu.edu/index.php?title=Simulations_of_Particle_Interactions_with_Matter&oldid=129967"