Exploratory Data Analysis: Gamma Spectroscopy in Python

, data science and data analysis can be closely related to physics and hardware. Not everything can run in the cloud, and some applications require the use of real things. In this article, I will show how to collect the data from a Radiacode 103G radiation detector and gamma spectrometer, and we will see what kind of information we can get from it. We will do this analysis in Python, and I will also show the gamma spectra of different objects that I got in the second-hand store. In the next part, I will use machine learning methods to detect object types and isotopes automatically.

All files collected for this article are available on Kaggle, and readers are also welcome to run all tests on their own. The link is added to the end of the page.

Let’s get started!

1. Hardware

As readers can guess from the top picture, we will talk about the radiation. Which is, interestingly, always around us, and the radiation level is never zero. Increased levels of radiation can be found in many places and objects, from vintage watches in the thrift store to airplane flights (because of cosmic rays, the radiation level during the flight is about 10x higher compared to ground level).

Simply speaking, there are mostly two kinds of radiation detectors:

• A Geiger-Müller tube. As its name suggests, it’s a tube filled with a mixture of gases. When a charged particle reaches the tube, the gas is ionized, and we can detect the short pulse. The higher the radiation level, the more pulses per minute we get. Radiation detectors often show values in CPM (counts per minute), which can be converted into Sieverts or other units. The Geiger tube is cheap and reliable; it is used in many radiation detectors.
• A scintillation detector. This detector is based on a special type of crystal, which generates light when a charged particle is detected. A scintillator has an important property—the intensity of the light is proportional to the particle’s energy. Because of that, we can not only detect the particles but can also determine which type of particles we get.

Obviously, from a hardware perspective, it’s easier said than done. In a scintillation crystal, only several photons can be emitted when a particle is detected. Earlier, those detectors had a 5-6 digit price, and were used only in labs and institutions. Nowadays, because of the progress in electronics, we can buy a scintillation detector for the price of a mid-level smartphone. This makes gamma spectroscopy analysis doable even for science enthusiasts with a moderate budget.

Let’s get into it and see how it works!

2. Collecting the Data

As was mentioned in the beginning, I will use a Radiacode 103G. It’s a portable device that can be used as a radiation detector and a gamma spectrometer — we can get both CPS (counts per second) and gamma spectrum values. Radiacode needs only the USB connection and has the size of a USB stick:

To get the data, I will use a radiacode open-source library. As a simple example, let’s collect the gamma spectrum within 30 seconds:

from radiacode import RadiaCode


rc = RadiaCode()
rc.spectrum_reset()
time.sleep(30)

spectrum = rc.spectrum()
print(spectrum.duration)
#> 0:00:30

print(spectrum.a0, spectrum.a1, spectrum.a2)
#> 24.524023056030273 2.2699732780456543 0.00043278629891574383

print(len(spectrum.counts))
#> 1024

print(spectrum.counts)
#> [0, 0, 0, 0, 2, 6, … 0, 0, 1]

I will explain the meaning of these fields in the following chapter when we start the analysis.

If we use the Radiacode in the Jupyter Notebook, we also need to close the connection at the end of the cell, otherwise the next run will get a USB “resource busy” error:

usb.util.dispose_resources(rc._connection._device)

We can also make a logger for both CPS (counts per second) and spectrum values. For example, let’s log gamma spectra into the CSV file every 60 seconds:

from radiacode import RadiaCode, RawData, Spectrum


SPECTRUM_READ_INTERVAL_SEC = 60
spectrum_read_time = 0

def read_forever(rc: RadiaCode):
   """ Read data from the device """
   while True:
       self.process_device_data(rc)
       time.sleep(0.3)

       t_now = time.monotonic()
       if t_now - spectrum_read_time >= SPECTRUM_READ_INTERVAL_SEC:
           self.process_spectrum_data(rc)
           spectrum_read_time = t_now

def process_device_data(rc: RadiaCode):
    """ Get CPS (counts per second) values """
    data = rc.data_buf()
    for record in data:
        if isinstance(record, RawData):
            dt_str = record.dt.strftime(self.TIME_FORMAT)
            log_str = f"{dt_str},,{int(record.count_rate)},"
            logging.debug(log_str)

def process_spectrum_data(rc: RadiaCode):
    """ Get spectrum data from the device """
    spectrum: Spectrum = rc.spectrum()
    save_spectrum_data(spectrum)

def save_spectrum_data(spectrum: Spectrum):
    """ Save spectrum data to the log """
    spectrum_str = ';'.join(str(x) for x in spectrum.counts)
    s_data = f"{spectrum.a0};{spectrum.a1};{spectrum.a2};{spectrum_str}"
    t_now = datetime.datetime.now()
    dt_str = t_now.strftime(self.TIME_FORMAT)
    log_str = f"{dt_str},,,{s_data}"
    logging.debug(log_str)


if __name__ == '__main__':
   rc = RadiaCode()
   rc.spectrum_reset()
   read_forever(rc)

Here, I get two types of data. A RawData contains CPS values, which can be used to get radiation levels in µSv/hour. I use them only to see that the device works; they are not used for analysis. A Spectrum data contains what we’re interested in – gamma spectrum values.

Collecting the data every 60 seconds allows me to see the dynamics of how the data is changing. Usually, the spectrum must be collected within several hours (the more the better), and I prefer to run this app on a Raspberry Pi. The output is saved as CSV, which we will process in Python and Pandas.

3. Data Analysis

3.1 Gamma Spectrum

To understand what kind of data we have, let’s print the spectrum again:

spectrum = rc.spectrum()
print(spectrum.duration)
#> 0:00:30

print(spectrum.a0, spectrum.a1, spectrum.a2)
#> 24.524023056030273 2.2699732780456543 0.00043278629891574383

print(len(spectrum.counts))
#> 1024

print(spectrum.counts)
#> [0, 0, 0, 0, 2, 6, … 0, 0, 1]

What did we get here? As was mentioned before, a scintillation detector gives us not only the number of particles but also their energy. The energy of the charged particle is measured in keV (kiloelectronvolts) or MeV (megaelectronvolts), where the electronvolt is the amount of kinetic energy of the particle. A Radiacode detector has 1024 channels, and the energy for each channel can be found using a simple formula:

Here, ch is a channel number, and a0, a1, and a2 are calibration constants, stored in the device.

During the specified time (in our case, 30 seconds), the Radiacode detects the particles and saves the result into channels as a simple arithmetic sum. For example, the value “2” in the 5th place means that 2 particles with the 33.61 keV energy were detected in channel 5.

Let’s draw the spectrum using Matplotlib:

def draw_simple_spectrum(spectrum: Spectrum):
   """ Draw spectrum from the Radiacode """
   a0, a1, a2 = spectrum.a0, spectrum.a1, spectrum.a2
   ch_to_energy = lambda ch: a0 + a1 * ch + a2 * ch**2

   fig, ax = plt.subplots(figsize=(12, 4))
   fig.gca().spines["top"].set_color("lightgray")
   fig.gca().spines["right"].set_color("lightgray")
   counts = spectrum.counts
   energy = [ch_to_energy(x) for x in range(len(counts))]
   # Bars
   plt.bar(energy, counts, width=3.0, label="Counts")
   # keV label values
   ticks_x = [
      ch_to_energy(ch) for ch in range(0, len(counts), len(counts) // 20)
   ]
   labels_x = [f"{ch:.1f}" for ch in ticks_x]
   ax.set_xticks(ticks_x, labels=labels_x)
   plt.xlim(energy[0], energy[-1])
   plt.title("Gamma-spectrum, Radiacode 103")
   plt.xlabel("Energy, keV")
   plt.legend()
   fig.tight_layout()
   fig.show()

draw_simple_spectrum(spectrum)

The output looks like this:

I collected the data in 30 seconds, and even within this short interval, we can already see a probability distribution! As usual in statistics, the longer the better, and by collecting the data within several hours, we can get good-visible results. It is also convenient that Radiacode collects the spectrum automatically. We can keep the device running autonomously for several hours or even days, then run the code and retrieve the spectrum.

For those readers who don’t have a Radiacode detector, I made two functions to save and load the spectrum to a JSON file:

def save_spectrum(spectrum: Spectrum, filename: str):
   """ Save spectrum to a file """
   duration_sec = spectrum.duration.total_seconds()
   data = {
       "a0": spectrum.a0,
       "a1": spectrum.a1,
       "a2": spectrum.a2,
       "counts": spectrum.counts,
       "duration": duration_sec,
   }
   with open(filename, "w") as f_out:
       json.dump(data, f_out, indent=4)
       print(f"File '{filename}' saved, duration {duration_sec} sec")

def load_spectrum(filename: str) -> Spectrum:
   """ Load spectrum from a file """
   with open(filename) as f_in:
       data = json.load(f_in)
       return Spectrum(
           a0=data["a0"], a1=data["a1"], a2=data["a2"],
           counts=data["counts"],
           duration=datetime.timedelta(seconds=data["duration"]),
       )

We can use it to get the data:

spectrum = load_spectrum("sp_background.json")
draw_simple_spectrum(spectrum)

As mentioned in the beginning, all collected files are available on Kaggle.

3.2 Isotopes

Now, we are approaching the fun part of the article. What is the purpose of the gamma spectrum? It turned out that different elements emit gamma rays with different energies. That allows us to tell what kind of object we have by observing the spectrum. This is a crucial difference between a simple Geiger-Müller tube and a scintillation detector. With a Geiger-Müller tube, we can tell that the object is radioactive and see that the level is, let’s say, 500 counts per minute. With a gamma spectrum, we can not only see the radiation level but also see why the object is radioactive and how the decay process is going.

How does it work in practice? Let’s say I want to measure radiation from a banana. Bananas naturally contain Potassium-40, which emits gamma rays with a 1.46 MeV energy. In theory, if I take a lot of bananas (really a lot because each banana has less than 0.5g of potassium) and place them in an isolated lead chamber to block the background radiation, I will see a 1.46 MeV peak on a spectrum.

Practically, it is often more complicated. Some materials, like uranium-238, may have a complex decay chain like this:

As we can see, uranium has a series of radioactive decays. It decays to thorium, thorium to radium, and so on. Every element has its gamma spectrum peak and a half-life time. As a result, all these peaks will be present on a spectrum at the same time (but with different intensities).

We can show isotope peaks with Matplotlib as well. As an example, let’s draw the K40 isotope line:

isotopes_k40 = [ ("K-40", 1460, "#0000FF55") ]

We need to add only two lines of code in Matplotlib to draw it on a spectrum:

# Isotopes
for name, en_kev, color in isotopes:
   plt.axvline(x = en_kev, color = color, label=name)
# Spectrum
plt.bar(energy, counts, width=3.0, label="Counts")

Now, let’s see how it works in action and do some experiments!

4. Experiments

I don’t work in a nuclear institution, and I don’t have access to official test sources like Cesium-137 or Strontium-90, used in a “big science.” However, it is not required, and some objects around us can be slightly radioactive.

Ideally, a test object must be placed in a lead chamber to reduce the background radiation. A thick layer of lead will reduce the background radiation to a much lower level. However, lead is heavy, the shipment can be expensive, and it is a toxic material that requires a lot of safety precautions. Instead, I will collect two spectra—the first for the object itself and the second for the background (without the object). As we will see, it is enough, and the difference will be visible.

Let’s combine all parts and draw both spectra and isotopes:

def draw_spectrum(
   spectrum: Spectrum, background: Spectrum, isotopes: List, title: str
):
   """ Draw the spectrum, the background, and the isotopes """
   counts = np.array(spectrum.counts) / spectrum.duration.total_seconds()
   counts_b = np.array(background.counts) / background.duration.total_seconds()

   a0, a1, a2 = spectrum.a0, spectrum.a1, spectrum.a2
   ch_to_energy = lambda ch: a0 + a1 * ch + a2 * ch**2

   # X-range
   x1, x2 = 0, 1024
   channels = list(range(x1, x2))
   energy = [ch_to_energy(x) for x in channels]

   fig, ax = plt.subplots(figsize=(12, 8))
   fig.gca().spines["top"].set_color("lightgray")
   fig.gca().spines["right"].set_color("lightgray")
   # Isotopes
   for name, en_kev, color in isotopes:
      plt.axvline(x = en_kev, color = color, label=name)
   # Bars
   plt.bar(energy, counts[x1:x2], width=3.0, label="Counts")
   plt.bar(energy, counts_b[x1:x2], width=3.0, label="Background")
   # X labels
   ticks_x = [ch_to_energy(ch) for ch in range(x1, x2, (x2 - x1) // 20)]
   labels_x = [f"{ch:.1f}" for ch in ticks_x]
   ax.set_xticks(ticks_x, labels=labels_x)
   plt.xlim(energy[0], energy[-1])
   plt.title(f"{title}, gamma-spectrum, Radiacode 103G")
   plt.xlabel("Energy, keV")
   plt.legend()
   fig.tight_layout()
   fig.show()

Different spectra can be collected during different time intervals. Because of that, I normalised the graph by dividing the count values by the total time, so we always see the counts per second on a graph.

Now, let’s start with experiments! As a reminder, all data files used in the tests are available on Kaggle.

4.1 Bananas

First, let’s answer the most asked question in nuclear physics – how radioactive are the bananas? A banana is a surprisingly difficult object to measure – it contains less than 1g of potassium-40, and it also contains 74% of water. As we can guess, the radiation from a banana is very low. First, I tried with a regular banana, and it did not work. Then I bought dried bananas in a supermarket, and only after that, I could see a small difference.

Using the methods created before, it’s easy to load a spectrum and see the result:

spectrum = load_spectrum("sp_bananas_dried.json")
background = load_spectrum("sp_bananas_background.json")
isotopes_k40 = [ ("K-40", 1460, "#0000FF55") ]

draw_spectrum(spectrum, background, isotopes_k40, title="Dried bananas")

The output looks like this:

As we can see, the difference is minuscule, and it is barely visible. Let’s calculate the difference caused by the Potassium-40 (1460 keV energy, which corresponds to a Radiacode channel number N=570):

ch, width = 570, 5
counts = np.array(spectrum.counts) / spectrum.duration.total_seconds()
counts_b = np.array(background.counts) / background.duration.total_seconds()
sum1 = sum(counts[ch - width:ch + width])
sum2 = sum(counts_b[ch - width:ch + width])
diff = 100*(sum1 - sum2)/sum(counts)

print(f"Diff: {diff:.3f}%")
#> Diff: 0.019%

Here, I compared the difference to the total number of particles, caused by background radiation. As we can see, the banana is only 0.019% more radioactive than the background! Obviously, this number is too small and cannot be seen with the naked eye just by watching the radiation counter.

4.2 Vintage Watch with a Radium Dial

Now, let’s test some “stronger” objects. Radium-226 was used for painting watch hands and dials from the 1910s to the 1960s. The mix of radium and a special lume allowed watches to glow in the dark. 
Many manufacturers were producing watches with radium paint, from cheap noname models to luxury ones like the Rolex Submariner with a >$40K modern price tag.

I bought this watch for testing in the vintage shop for about $20:

I also used a UV light to show how the watch was glowing in the dark when it was made (nowadays, the lume is depleted, and without UV, it does not glow anymore).

Let’s draw the spectrum:

spectrum = load_spectrum("sp_watch.json")
background = load_spectrum("sp_background.json")
isotopes_radium = [
   ("Ra-226", 186.21, "#AA00FF55"),
   ("Pb-214", 242.0,  "#0000FF55"),
   ("Pb-214", 295.21, "#0000FF55"),
   ("Pb-214", 351.93, "#0000FF55"),
   ("Bi-214", 609.31, "#00AAFF55"),
   ("Bi-214", 1120.29, "#00AAFF55"),
   ("Bi-214", 1764.49, "#00AAFF55"),
]

draw_spectrum(spectrum, background, isotopes_radium, title="Vintage watch")

The data was collected within 3.5 hours, and the output looks like this:

Lots of processes are going on here, and we can see the parts of the radium decay chain:

Radium-226 (1602-year half-life) yields an alpha particle and the radon gas (3.82-day half-life), but the radon itself is also an alpha emitter and does not produce gamma rays. Instead, we can see some daughter products as bismuth and lead.

As an aside question: is it dangerous to have a watch with a radium dial? Generally not, these watches are safe if the case and the glass are not damaged. We can see many decay products on the spectrum, but in the same way as with bananas, this spectrum was collected within several hours, and the real radiation level from this watch is relatively small.

4.3 Uranium Glass

Uranium glass has a surprisingly long history. Uranium ore was used for glass manufacturing since the 19th century, long before even the word “radioactivity” became known. Uranium glass was produced in huge amounts and can now be found in almost every vintage store. The production stopped in the 1950s; after that, all the uranium was mostly used for industrial and military purposes.

Uranium glass mostly contains only a tiny fraction of uranium; these objects are safe to keep in the cabinet (still, I would not use it for eating:). The amount of radiation produced by the glass is also small.

This Victorian glass was made in the 1890s, and it will be our test object:

Let’s see the spectrum:

spectrum = load_spectrum("sp_uranium.json")
background = load_spectrum("sp_background.json")
isotopes_uranium = [
   ("Th-234", 63.30, "#0000FF55"),
   ("Th-231", 84.21, "#0000FF55"),
   ("Th-234", 92.38, "#0000FF55"),
   ("Th-234", 92.80, "#0000FF55"),
   ("U-235", 143.77, "#00AAFF55"),
   ("U-235", 185.72, "#00AAFF55"),
   ("U-235", 205.32, "#00AAFF55"),
   ("Pa-234m", 766.4, "#0055FF55"),
   ("Pa-234m", 1000.9, "#0055FF55"),
]

draw_spectrum(spectrum, background, isotopes_uranium, title="Uranium glass")

As was mentioned, the radiation level from the uranium glass is not high. Most of the peaks are small, and I used a log scale to see them better. The result looks like this:

Here, we can see some peaks from thorium and uranium; these elements are present in the uranium decay chain.

4.4 “Quantum” Pendant

Readers may think that radioactive objects were produced only in the “dark ages,” a long time ago, when people did not know about radiation. Funny enough, it is not always true, and many “Negative Ion” related products are still available today.

I bought this pendant for $10 on eBay, and according to the seller’s description, it can be used as a negative ions source to “improve health”:

An ionizing radiation indeed can produce ions. The medical properties of those ions are out of the scope of this article. Anyway, it is a good test object for gamma spectroscopy – let’s draw the spectrum and see what is inside:

spectrum = load_spectrum("sp_pendant.json")
background = load_spectrum("sp_background.json")
isotopes_thorium = [
   ("Pb-212", 238.63, "#0000FF55"),
   ("Ac-228", 338.23, "#0000FF55"),
   ("TI-208", 583.19, "#0000FF55"),
   ("AC-228", 911.20, "#0000FF55"),
   ("AC-228", 968.96, "#0000FF55"),
]

draw_spectrum(
   spectrum, background, isotopes_thorium, title="'Quantum' pendant"
)

The result looks like this:

According to gammaspectacular.com, this spectrum belongs to Thorium-232. Obviously, it was not written in the product description, but with a gamma spectrum, we can see it without object disassembly or any chemical tests!

5. Matplotlib Animation (Bonus)

As a bonus for readers patient enough to read up to this point, I will show how to make an animated graph in Matplotlib.

As a reminder, at the beginning of the article, I collected spectra from a Radiacode device every minute. We can use this data to animate a spectrum collection process.

First, let’s load the data into the dataframe:

def get_spectrum_log(filename: str) -> pd.DataFrame:
   """ Get dataframe from a CSV-file """
   df = pd.read_csv(
       filename, header=None, index_col=False,
       names=["datetime", "temp_c", "cps", "spectrum"],
       parse_dates=["datetime"]
   )
   return df.sort_values(by='datetime', ascending=True)

def get_spectrum_data(spectrum: str) -> np.array:
   """ Spectrum data: a0, a1, a2, 1024 energy values
   Example:
   9.572;2.351;0.000362;0;2;0;0; ... """
   values = spectrum.split(";")[3:]
   return np.array(list(map(int, values)))

def process_spectrum_file(filename: str) -> pd.DataFrame:
   df = get_spectrum_log(filename)
   df = df[
       df["spectrum"].notnull()
   ][["datetime", "spectrum"]].copy()

   df["spectrum"] = df["spectrum"].map(get_spectrum_data)
   return df

df = process_spectrum_file("spectrum-watch.log")
display(df)

The output looks like this:

As we can see, each row in the observation contains the same fields as we used before for a single Spectrum object.

The animation process is relatively straightforward. First, we need to save a Matplotlib’s “bar” object in the variable:

plt.bar(energy, counts_b[x1:x2], width=3.0, label="Background")
bar = plt.bar(energy, counts[x1:x2], width=3.5, label="Counts")

Then we can use an update function that takes spectra from a dataframe:

def update(frame: int):
   """ Frame update callback """
   index = frame * df.shape[0] // num_frames
   counts = df["spectrum"].values[index]
   # Update bar
   for i, b in enumerate(bar):
       b.set_height(counts[x1:x2][i])
   # Update title 
   ax.set_title(...)

Now, we can save the animation into the GIF file:

import matplotlib.animation as animation

# Animate
anim = animation.FuncAnimation(
  fig=fig, func=update, frames=num_frames, interval=100
)
writer = animation.PillowWriter(fps=5)
anim.save("spectrum-watch.gif", writer=writer)

The output looks like this:

Here, we can see how the number of particles increases during the collection time.

6. Conclusion

In this article, I described the basics of gamma spectroscopy analysis with a Radiacode scintillation detector, Python, and Matplotlib. As we can see, the data itself is simple, however, a lot of invisible processes are going “under the hood.” Using different test objects, we were able to see spectra of different elements of radioactive decay as Bismuth-214 or Protactinium-234. And, amazingly, these tests can be done at home, and detectors like Radiacode are available for science enthusiasts for the price of a mid-level smartphone.

This basic analysis allows us to go further. If there is some public interest in this topic, in the next article I will use machine learning to automate the isotopes detection. Stay tuned.

For those readers who would like to perform the tests on their own, all datasets are available on Kaggle. If someone would like to do the tests with their own hardware, a Radiacode manufacturer kindly provided a 10% discount code “DE10” for buying the device (disclaimer: I don’t have any profit or other commercial interest from these sales).

All readers are also welcome to connect via LinkedIn, where I periodically publish smaller posts that are not big enough for a full article. If you want to get the full source code for this and other articles, feel free to visit my Patreon page.

Thanks for reading.