# Generate a histogram of the number of neutrons emitted at
# angle "theta" to the direction of a proton beam striking a
# lithium-7 target. Run the code in the LWDAQ Toolmaker. It
# will take a few minutes.
# Set up the histogram with counts zero and zero for the number
# of neutrons emitted in an omnidirectional distribution and the
# number in the distribution when accounting for the momentum
# of the proton.
set pi 3.14159265359
set histogram [list]
for {set theta 0} {$theta <= 180} {incr theta} {lappend histogram "$theta 0 0"}
# Generate neutrons at random.
for {set i 0} {$i < 500000} {incr i} {
# Send a neutron off in a random direction with respect to the
# proton direction. The probability density function here is
# uniform.
set a [LWDAQ_random 0.0 180.0]
# We want to generate the probability distribution of angle
# when we have uniform distribution with respect to area
# on a sphere. This angular distribution is like sin(a),
# so we ignore our neutron if a random 0-1 is greater than
# sin(a). Our histogram of neutrons at angle a will look
# like a sinusoid.
if {[LWDAQ_random 0.0 1.0] > sin($a*$pi/180)} {continue}
# The proton is slowing down as it passes through the
# lithium target. Its energy when it reacts with a
# lithium atom can be anything from 0-4 MeV.
set Ep [LWDAQ_random 0.0 4.0]
# The reaction energy is the increase in mass. We calculated
# this value using a table of atomic masses.
set Er 1.64
# The parameter alpha is the ratio of the speed of the
# beryllium-8 atom to the speed of the neutron it later
# ejects, when it becomes beryllium-7. This fraction
# depends upon the energy of the proton and the kinetic
# energy of the neutron and beryllium-7 combined. We start
# by calculating the energy available for the movement of
# the beryllium-7 and neutron in the beryllium-8 frame. If
# this is greater than zero (slightly greater, actually), we
# proceed with our calculation. Otherwise, we continue to
# the next collision.
set Ex [expr $Ep - 64.0*$Er/56.0]
if {$Ex <= 0.1 } {continue}
set alpha [expr sqrt($Ep/$Ex)/7.0]
# The tangent of the momentum-corrected direction of the
# neutron we obtain by adding alpha to the component of its
# velocity in the direction of the proton.
set x [expr sin($a * $pi / 180.0) / (cos($a * $pi / 180.0) + $alpha)]
# We obtain the arctangent to get the corrected neutron
# direction, and handle negative angles.
set b [expr 180.0 / $pi * atan($x)]
if {$b < 0.0} {set b [expr 180.0 + $b]}
# We round both a and by for our histogram. Each bin is
# one degree wide.
set b [expr round($b)]
set a [expr round($a)]
# We use a and b to obtain, increment, and save our
# counts of neutrons with angle a when uncorrected and
# angle b when corrected.
lset histogram $a 1 [expr [lindex $histogram $a 1] + 1]
lset histogram $b 2 [expr [lindex $histogram $b 2] + 1]
# Print an occasional dot.
if {[LWDAQ_random 0 1000] == 1} {LWDAQ_print -nonewline $t "."}
# We break out of this loop if the execution window closes.
LWDAQ_support
if {![winfo exists $t]} {break}
}
# Print the histogram to the execution window. We will cut
# and paste into a spreadsheet to obtain plots.
LWDAQ_print $t
foreach c $histogram {LWDAQ_print $t $c}