(The package files and example scripts have all been updated to v1.1 on August 7, 2017.)

The files required to use the basic package are given below. All scripts require that these files be placed in the same directory where the script is running.

FILES:

• radioactive.py
• al.txt
• pb.txt
• nai.txt
• ge.txt

As an example of how the package functions, we provide a bare-bones sample script, 'runit.py'.

The function source() takes as inputs:

• act, the activity of the source in μCi. The default is _act _= 1 μCi.
• Tf the total time over which the source emits radiation. The default is Tf= 1s.
• axisand cone, which allow for emitted rays to be accepted only if they fall inside a cone around an axis. cone should be an angle between 0 and π, describing the maximum angle between an emitted ray and the indicated axis. If these are not given, rays will come out isotropically in every direction. If only cone is given, the angle will be relative to the z-axis.
• posn, the position of the source. The default is posn = [0, 0, 0].
• isotope, the radioactive isotope indicated as element followed by mass number, for example ‘Na22’. If the isotope indicated has not been preprogrammed or no isotope is given, the default is isotope = ‘Na22’.

The output is a list of events, organized as an array. Each element of this array corresponds to an emitted gammaray, and contains the information of time elapsed between the beginning of the information and emission of this ray, energy, position and direction of propagation.

The function source_compton() has the same inputs and outputs, but increased efficiency when there is only one decay mode and we are interested in gamma rays emitted in a small cone, such as for simulation of a Compton scattering experiment. This is achieved by rescaling the activity by

and generating propagating directions only inside the intended cone instead of checking after generating the direction. As currently implemented, source_compton()does not work for sources with more than one decay mode.

This function takes the following inputs:

• solid, indicating the shape of the scatterer. The default (and at the moment only supported option) is solid = ‘cylinder’.
• center and axis, two vectors indicating the position of the cylinder. The defaults are center = [0,0,5] and axis = [0,0,1].
• radius and length, two scalars indicating the size of the cylinder. Defaults are radius = 1 cm and length = 1 cm.
• eventsin, the list of incoming events, an array with the same description as the output of source.
• material, a string indicating the chemical symbol of the material the scatterer is made of. The default is material = ‘Al’. A text file with the corresponding interaction cross sections should be present in the same folder as the scripts, as detailed in Sec. 2.1.
• rho, the density of the scattering material. This is necessary only if the material is not Al or Pb.

For each gamma ray in eventsin, we start by adding it to an empty array named liverays. We then find the intersection of the path of that gammawith the volume of the scatterer. If there is no intersection, the gammaray is simply added to eventsout. If there is an intersection, the program finds the length of the intersection d and obtains the cross section σ for the required energy from the input data.

A uniform random number $u\in (0,1)$ is then drawn. If $u < e^{-\sigma\rho d}$ , there is no interaction, and the gamma is propagated to the final intersection point and added to eventsout. Otherwise, an interaction happens after traveling a distance of $-(\log u) / \sigma\rho$ inside the scatterer. The gammais propagated to that point. A new random number is then drawn to determine if the interaction is Compton, photoelectric or pair production. If it is photoelectric, the gammais absorbed, and we move on to the next ray. If Compton, see Sec. 2.2 for an explanation of the algorithm followed; the outgoing ray is then added back to liverays. If pair production, we draw a random direction and add two rays of energy 511 keV, one moving in that directionand one in the opposite.

Once liveraysis empty, we move on to the next element of eventsin. Once all gammarays have been scattered, the array eventsout is the result.

It should be noted that before any other calculations, we find the transformation from the lab frame to a frame where the axis of the scatterer points along the z-axis and the center of the detector is at the origin, and transform all positions and velocities from the lab to the scatterer frame,

and

Before returning eventsout, we perform the inverse transformation on all the vectors back to the lab frame,

and

This function has a similar structure to scatterer(). The inputs it takes are exactly the same, with the difference that the default material is material= “NaI”.

The algorithm is also similar, with the difference that for each ray in eventsin we create a variable corresponding to how much energy is dumped into the detector. In a photoelectric interaction we suppose all energy from the gamma ray is absorbed, in Compton we take all the energy difference to be absorbed, and for pair production we assume $E -1.022$ MeV is absorbed.

For each gammaray which dumps any amount of energy, we now calculate the energy detected depending on properties of the detector. For a Ge solid state detector, energy deposited in the crystal goes into the production of electron-hole pairs. Each electron hole pair takes an energy $\epsilon$ to create, so if an amount of energy $E$ is deposited in the crystal, on average $N = E/\epsilon$ pairs are created. In scintillation detectors such as NaI, energy goes into producing scintillation light, which is composed of photons of energy $\epsilon$ . Once again, $N = E/\epsilon$ photons are created. Note these $\epsilon$ are different quantities, but that will not be important in our derivation. Since these processes are approximately Poissonian, the variance in $N$ is proportional to $N$ :$\sigma_N^2 = kN$ for some constant $k$ which depends on the specific detector. (Note that the process in solid state detectors is, strictly speaking, not Poissonian for which $k=1$ , but the variance has the same $N$ -dependence with a different $k$ ). The standard deviation of the distribution of detected energy, given the amount of deposited energy $E$ , then depends on $E$ as

We can then measure $\sigma_0$  at some energy $E_0$  experimentally, and from there get the general distribution. The program uses values obtained for $E_0 = 662 \;\mathrm{keV}$  of $\sigma_0 = 20.4 \;\mathrm{keV}$  for NaI (measured on the 1.5“ x 1.5” SpecTech crystals) and $\sigma_0 = 1.99 \;\mathrm{keV}$  keV for Ge. The detected energy is then obtained from the deposited energy by multiplying $E$  by a draw of a Gaussian with mean 1 and standard deviation obtained as described above. This array is returned by the function, which can be then be plotted via histogram or otherwise analyzed.

The program draws the cross sections for each interaction process of a gamma ray with matter from a table provided. A good source for these tables is NIST. The program assumes the order of the columns is as follows:

• Energy in MeV;
• Cross section for Rayleigh scattering in cm2/g;
• Cross section for Compton scattering in cm2/g;
• Cross section for photoelectric effect in cm2/g;
• Cross section for pair production in nuclear field in cm2/g;
• Cross section for pair production in electron field in cm2/g;
• Total interaction cross section including Rayleigh scattering in cm2/g; and
• Total interaction cross section not including Rayleigh scattering in cm2/g.

At the moment, the two different pair production cross sections are simply summed, and Rayleigh scattering is ignored. The program interpolates linearly between the energies on the table, so the table should be filled closely enough for this to be accurate.

The input file should be called 'x.txt', where ‘x’ is given to the scatterer and detector functions with the material option. For example, if scatterer is called with _material _= ‘Al’, the program will look for the cross section data in the file 'Al.txt'. Unless the density of the material is preprogrammed, the function should also be called with the rho option specifying the density of the material in g/cm3. If the density is not preprogrammed and not specified, a warning will be printed. Note the program skips the first ten lines of the input file, so the data should start on the eleventh line of the file.

If a gamma ray is ever analyzed with energy below the energy of the first line in the data file, the program prints a warning:

Warning: Energy is lower than the lowest value on the table provided.

If this happens often, more data should be added for lower energies. If the data starts at low enough energies, this should never happen, as the photoelectric cross section dominates before we get to a very low energy gammaray. 

In this section, we will describe how the script implements the Klein-Nishina formula.

First, recall the relevant formulas for Compton scattering. When photons of energy $E$ are scattered by an angle $\theta$, the ratio between the final and initial energies is given by

where $m_ec^2 = 511 \;\mathrm{keV}$  is the electron rest mass. The differential cross section for scattering is given by 

where in the second line we just grouped all the multiplicative constants into a normalization $A$ that we will not need. (Therefore, we will not be too careful about keeping the meaning of $A$ invariant from equation to equation – it is simply the required normalization constant.) For low energies $E \ll m_ec^2$, $r(E,\theta ) \approx 1$  and the differential cross section simplifies into 

In this package, we take the total Compton cross section from existing data, and we only need to use the equations above to draw an angle$\theta$ according to the correct probability distribution and from $\theta$ calculate the correct final energy $E_f$ . We will do this by obtaining a $\theta_L$ obeying the distribution from the low-energy approximation, Eq. (7), by inverse transform sampling and then using that $\theta_L$ and rejection sampling to end with a $\theta$ drawn from the full Klein-Nishina distribution, Eq. (6).

The principle of inverse transform sampling is that given a probability distribution $P(x)$, we start by calculating the cumulative distribution function$f(x) = P(X < x)$. We then invert $f(x)$ to obtain $f^{-1}(u)$. We can now use this to transform a uniform distribution in the interval (0,1) to a variable with the distribution $P$ by setting $x = f^{-1}(u)$ . In brief, this works because the probability of this algorithm giving an $x$ smaller than any given value $x_0$ is given by $f(x_0)$ , which is the correct distribution. Note that we cannot use this method directly on the Klein-Nishina formula because we cannot invert the cumulative distribution function (CDF) in closed form. To use this method, the first step is therefore to find the CDF corresponding to Eq. (7) as a probability distribution for the scattering angle $\theta$. Note that $\mathrm{d}\Omega = \sin\theta \;\mathrm{d}\theta\;\mathrm{d}\phi$ , so the probability distribution as a function of $\theta$ is given by

Integrating this distribution to find the CDF, we have

We will take this as the CDF for $\cos\theta$, $f(x) = 4 - 3x - x^3)/8$ . Inverting this function (and keeping in mind that x should be real), we get

where

We can then implement the inverse transform sampling for the distribution Eq. (7) by drawing a uniform random number $u \in (0,1)$ and setting $\theta_L = f^{-1}(u)$.

Given this $\theta_L$ , we now use rejection sampling to arrive at a $\theta$ drawn from Eq. (6). This method transforms between a random variable drawn from a distribution $P_1(x)$ to a random variable drawn from a distribution $P_2(x)$ by adding an additional step where we accept the drawing with probability$P_2(x)/kP_1(x)$, where $k$ is a constant such that $P_2(x) \lt kP_1(x)$ for all $x$. If the drawing is rejected, we draw a new sample until we the drawing is accepted. Clearly this method is more efficient if the distributions $P_1$ and $kP_2$ are close to each other. In this instance, the probability of acceptance for an angle $\theta_L$ drawn from Eq. (7) is 

Note that for any $E$ and $\theta$ this is smaller than 1, as required.

gamma_cross_sections.py or gamma_cross_sections.ipynb

compton_scattering.py or compton_scattering.ipynb

compton_edge.py or compton_edge.ipynb

In the Relativistic Electrons experiment, one of the tasks is to measure the energy of Compton edges in spectra. Unlike peaks (which are well defined Gaussians), the “position” of the edge is somewhat ambiguous. In an ideal detector, the edge would be sharp, but a realistic detector smears that edge out. In addition, the edge is asymmetric; to the left of the edge, the count rate in the Compton plateau naturally rises just before falling off (due to the Klein-Nishina cross section), while to the right, there is influence due to some multiple Compton scatter (without photoelectric absorption) events. Though we can state the causes, we do not have a good mathematical model for the shape of this resulting edge.

Therefore, we can use this program to simulate spectra at different energies and for different detector types in order to understand which point along any observed Compton edge shape corresponds to the actual edge energy. In this way, we can justify (or modify) our estimates for the experimentally-measured Compton edge values and uncertainties despite the fact that the edge is typically an extended, spread-out feature.

The program creates simulated spectra and overlays these with lines indicating the exact energies of peak centers and Compton edges. One can zoom in on the features to see broadening effects due to detector resolution. The default settings are 100 sec runs for 1 μCi sources placed 3 cm from the face of a 1.5“ x 1.5” cylindrical detector placed behind a 1 mm Al “window” (simulating a typical thin Al crystal housing). The program investigates Cs-137 and Na-22 in NaI and Ge detectors, but is customizable to other radioactive sources.

Figure x: Example Monte Carlo output plots showing peaks and edges. On the left is the full spectrum of Na-22 in a NaI detector. On the right is a detail plot of the Compton edge of the 662 keV line in Cs-137 for a Ge detector.