Exploratory Data Analysis: Gamma Spectroscopy in Python

, knowledge science and knowledge evaluation could be intently associated to physics and {hardware}. Not all the pieces can run within the cloud, and a few purposes require the usage of actual issues. On this article, I’ll present acquire the info from a Radiacode 103G radiation detector and gamma spectrometer, and we are going to see what sort of info we are able to get from it. We are going to do that evaluation in Python, and I may also present the gamma spectra of various objects that I bought within the second-hand retailer. Within the subsequent half, I’ll use machine studying strategies to detect object sorts and isotopes robotically.

All information collected for this text can be found on Kaggle, and readers are additionally welcome to run all checks on their very own. The hyperlink is added to the top of the web page.

Let’s get began!

1. {Hardware}

As readers can guess from the highest image, we are going to speak in regards to the radiation. Which is, curiously, all the time round us, and the radiation degree is rarely zero. Elevated ranges of radiation could be discovered in lots of locations and objects, from classic watches within the thrift retailer to airplane flights (due to cosmic rays, the radiation degree throughout the flight is about 10x increased in comparison with floor degree).

Merely talking, there are largely two sorts of radiation detectors:

• A Geiger-Müller tube. As its identify suggests, it’s a tube full of a mix of gases. When a charged particle reaches the tube, the gasoline is ionized, and we are able to detect the brief pulse. The upper the radiation degree, the extra pulses per minute we get. Radiation detectors typically present values in CPM (counts per minute), which could be transformed into Sieverts or different models. The Geiger tube is reasonable and dependable; it’s utilized in many radiation detectors.
• A scintillation detector. This detector is predicated on a particular sort of crystal, which generates gentle when a charged particle is detected. A scintillator has an necessary property—the depth of the sunshine is proportional to the particle’s vitality. Due to that, we can’t solely detect the particles however may also decide which sort of particles we get.

Clearly, from a {hardware} perspective, it’s simpler mentioned than executed. In a scintillation crystal, solely a number of photons could be emitted when a particle is detected. Earlier, these detectors had a 5-6 digit worth, and had been used solely in labs and establishments. These days, due to the progress in electronics, we are able to purchase a scintillation detector for the worth of a mid-level smartphone. This makes gamma spectroscopy evaluation doable even for science fanatics with a average finances.

Let’s get into it and see the way it works!

2. Gathering the Information

As was talked about to start with, I’ll use a Radiacode 103G. It’s a transportable system that can be utilized as a radiation detector and a gamma spectrometer — we are able to get each CPS (counts per second) and gamma spectrum values. Radiacode wants solely the USB connection and has the scale of a USB stick:

To get the info, I’ll use a radiacode open-source library. As a easy instance, let’s acquire the gamma spectrum inside 30 seconds:

from radiacode import RadiaCode


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

spectrum = rc.spectrum()
print(spectrum.length)
#> 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’ll clarify the that means of those fields within the following chapter once we begin the evaluation.

If we use the Radiacode within the Jupyter Pocket book, we additionally want to shut the connection on the finish of the cell, in any other case the following run will get a USB “useful resource busy” error:

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

We are able to additionally make a logger for each CPS (counts per second) and spectrum values. For instance, let’s log gamma spectra into the CSV file each 60 seconds:

from radiacode import RadiaCode, RawData, Spectrum


SPECTRUM_READ_INTERVAL_SEC = 60
spectrum_read_time = 0

def read_forever(rc: RadiaCode):
   """ Learn knowledge from the system """
   whereas 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 """
    knowledge = rc.data_buf()
    for file in knowledge:
        if isinstance(file, RawData):
            dt_str = file.dt.strftime(self.TIME_FORMAT)
            log_str = f"{dt_str},,{int(file.count_rate)},"
            logging.debug(log_str)

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

def save_spectrum_data(spectrum: Spectrum):
    """ Save spectrum knowledge to the log """
    spectrum_str = ';'.be a part of(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)

Right here, I get two kinds of knowledge. A RawData incorporates CPS values, which can be utilized to get radiation ranges in µSv/hour. I take advantage of them solely to see that the system works; they don’t seem to be used for evaluation. A Spectrum knowledge incorporates what we’re all in favour of – gamma spectrum values.

Gathering the info each 60 seconds permits me to see the dynamics of how the info is altering. Normally, the spectrum have to be collected inside a number of hours (the extra the higher), and I desire to run this app on a Raspberry Pi. The output is saved as CSV, which we are going to course of in Python and Pandas.

3. Information Evaluation

3.1 Gamma Spectrum

To know what sort of knowledge now we have, let’s print the spectrum once more:

spectrum = rc.spectrum()
print(spectrum.length)
#> 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 right here? As was talked about earlier than, a scintillation detector provides us not solely the variety of particles but additionally their vitality. The vitality of the charged particle is measured in keV (kiloelectronvolts) or MeV (megaelectronvolts), the place the electronvolt is the quantity of kinetic vitality of the particle. A Radiacode detector has 1024 channels, and the vitality for every channel could be discovered utilizing a easy system:

Right here, ch is a channel quantity, and a0, a1, and a2 are calibration constants, saved within the system.

Throughout the specified time (in our case, 30 seconds), the Radiacode detects the particles and saves the outcome into channels as a easy arithmetic sum. For instance, the worth “2” within the fifth place implies that 2 particles with the 33.61 keV vitality had been detected in channel 5.

Let’s draw the spectrum utilizing 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
   vitality = [ch_to_energy(x) for x in range(len(counts))]
   # Bars
   plt.bar(vitality, 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(vitality[0], vitality[-1])
   plt.title("Gamma-spectrum, Radiacode 103")
   plt.xlabel("Vitality, keV")
   plt.legend()
   fig.tight_layout()
   fig.present()

draw_simple_spectrum(spectrum)

The output seems to be like this:

I collected the info in 30 seconds, and even inside this brief interval, we are able to already see a chance distribution! As ordinary in statistics, the longer the higher, and by accumulating the info inside a number of hours, we are able to get good-visible outcomes. It is usually handy that Radiacode collects the spectrum robotically. We are able to hold the system working autonomously for a number of hours and even days, then run the code and retrieve the spectrum.

For these readers who don’t have a Radiacode detector, I made two features to save and cargo the spectrum to a JSON file:

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

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

We are able to use it to get the info:

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

As talked about to start with, all collected information can be found on Kaggle.

3.2 Isotopes

Now, we’re approaching the enjoyable a part of the article. What’s the function of the gamma spectrum? It turned out that completely different parts emit gamma rays with completely different energies. That enables us to inform what sort of object now we have by observing the spectrum. It is a essential distinction between a easy Geiger-Müller tube and a scintillation detector. With a Geiger-Müller tube, we are able to inform that the item is radioactive and see that the extent is, let’s say, 500 counts per minute. With a gamma spectrum, we can’t solely see the radiation degree but additionally see why the item is radioactive and the way the decay course of goes.

How does it work in observe? Let’s say I wish to measure radiation from a banana. Bananas naturally include Potassium-40, which emits gamma rays with a 1.46 MeV vitality. In concept, if I take loads of bananas (actually so much as a result of every banana has lower than 0.5g of potassium) and place them in an remoted lead chamber to dam the background radiation, I’ll see a 1.46 MeV peak on a spectrum.

Virtually, it’s typically extra sophisticated. Some supplies, like uranium-238, might have a posh decay chain like this:

As we are able to see, uranium has a sequence of radioactive decays. It decays to thorium, thorium to radium, and so forth. Each factor has its gamma spectrum peak and a half-life time. In consequence, all these peaks will likely be current on a spectrum on the identical time (however with completely different intensities).

We are able to present isotope peaks with Matplotlib as nicely. For instance, let’s draw the K40 isotope line:

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

We have to add solely two traces of code in Matplotlib to attract it on a spectrum:

# Isotopes
for identify, en_kev, coloration in isotopes:
   plt.axvline(x = en_kev, coloration = coloration, label=identify)
# Spectrum
plt.bar(vitality, counts, width=3.0, label="Counts")

Now, let’s see the way it works in motion and do some experiments!

4. Experiments

I don’t work in a nuclear establishment, and I don’t have entry to official check sources like Cesium-137 or Strontium-90, utilized in a “large science.” Nevertheless, it isn’t required, and a few objects round us could be barely radioactive.

Ideally, a check object have to be positioned in a lead chamber to cut back the background radiation. A thick layer of lead will scale back the background radiation to a a lot decrease degree. Nevertheless, lead is heavy, the cargo could be costly, and it’s a poisonous materials that requires loads of security precautions. As a substitute, I’ll acquire two spectra—the primary for the item itself and the second for the background (with out the item). As we are going to see, it’s sufficient, and the distinction will likely be seen.

Let’s mix all components and draw each spectra and isotopes:

def draw_spectrum(
   spectrum: Spectrum, background: Spectrum, isotopes: Checklist, title: str
):
   """ Draw the spectrum, the background, and the isotopes """
   counts = np.array(spectrum.counts) / spectrum.length.total_seconds()
   counts_b = np.array(background.counts) / background.length.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 = record(vary(x1, x2))
   vitality = [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 identify, en_kev, coloration in isotopes:
      plt.axvline(x = en_kev, coloration = coloration, label=identify)
   # Bars
   plt.bar(vitality, counts[x1:x2], width=3.0, label="Counts")
   plt.bar(vitality, 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(vitality[0], vitality[-1])
   plt.title(f"{title}, gamma-spectrum, Radiacode 103G")
   plt.xlabel("Vitality, keV")
   plt.legend()
   fig.tight_layout()
   fig.present()

Completely different spectra could be collected throughout completely different time intervals. Due to that, I normalised the graph by dividing the rely values by the entire time, so we all the time see the counts per second on a graph.

Now, let’s begin with experiments! As a reminder, all knowledge information used within the checks can be found on Kaggle.

4.1 Bananas

First, let’s reply essentially the most requested query in nuclear physics – how radioactive are the bananas? A banana is a surprisingly troublesome object to measure – it incorporates lower than 1g of potassium-40, and it additionally incorporates 74% of water. As we are able to guess, the radiation from a banana may be very low. First, I attempted with an everyday banana, and it didn’t work. Then I purchased dried bananas in a grocery store, and solely after that, I may see a small distinction.

Utilizing the strategies created earlier than, it’s simple to load a spectrum and see the outcome:

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 seems to be like this:

As we are able to see, the distinction is minuscule, and it’s barely seen. Let’s calculate the distinction brought on by the Potassium-40 (1460 keV vitality, which corresponds to a Radiacode channel quantity N=570):

ch, width = 570, 5
counts = np.array(spectrum.counts) / spectrum.length.total_seconds()
counts_b = np.array(background.counts) / background.length.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%

Right here, I in contrast the distinction to the entire variety of particles, brought on by background radiation. As we are able to see, the banana is just 0.019% extra radioactive than the background! Clearly, this quantity is simply too small and can’t be seen with the bare eye simply by watching the radiation counter.

4.2 Classic Watch with a Radium Dial

Now, let’s check some “stronger” objects. Radium-226 was used for painting watch palms and dials from the 1910s to the Sixties. The combination of radium and a particular lume allowed watches to glow at midnight. 
Many producers had been producing watches with radium paint, from low-cost noname fashions to luxurious ones just like the Rolex Submariner with a >$40K trendy price ticket.

I purchased this look ahead to testing within the classic store for about $20:

I additionally used a UV gentle to point out how the watch was glowing at midnight when it was made (these days, the lume is depleted, and with out UV, it doesn’t 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="Classic watch")

The information was collected inside 3.5 hours, and the output seems to be like this:

A lot of processes are occurring right here, and we are able to see the components of the radium decay chain:

Radium-226 (1602-year half-life) yields an alpha particle and the radon gasoline (3.82-day half-life), however the radon itself can also be an alpha emitter and doesn’t produce gamma rays. As a substitute, we are able to see some daughter merchandise as bismuth and lead.

As an apart query: is it harmful to have a watch with a radium dial? Typically not, these watches are secure if the case and the glass will not be broken. We are able to see many decay merchandise on the spectrum, however in the identical means as with bananas, this spectrum was collected inside a number of hours, and the actual radiation degree from this watch is comparatively small.

4.3 Uranium Glass

Uranium glass has a surprisingly lengthy historical past. Uranium ore was used for glass manufacturing because the nineteenth century, lengthy earlier than even the phrase “radioactivity” grew to become recognized. Uranium glass was produced in big quantities and may now be present in nearly each classic retailer. The manufacturing stopped within the Fifties; after that, all of the uranium was largely used for industrial and army functions.

Uranium glass largely incorporates solely a tiny fraction of uranium; these objects are secure to maintain within the cupboard (nonetheless, I might not use it for consuming:). The quantity of radiation produced by the glass can also be small.

This Victorian glass was made within the Eighteen Nineties, and will probably be our check 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 talked about, the radiation degree from the uranium glass will not be excessive. A lot of the peaks are small, and I used a log scale to see them higher. The outcome seems to be like this:

Right here, we are able to see some peaks from thorium and uranium; these parts are current within the uranium decay chain.

4.4 “Quantum” Pendant

Readers might imagine that radioactive objects had been produced solely within the “darkish ages,” a very long time in the past, when individuals didn’t find out about radiation. Humorous sufficient, it isn’t all the time true, and lots of “Unfavorable Ion” associated merchandise are nonetheless accessible at present.

I purchased this pendant for $10 on eBay, and in line with the vendor’s description, it may be used as a adverse ions supply to “enhance well being”:

An ionizing radiation certainly can produce ions. The medical properties of these ions are out of the scope of this text. Anyway, it’s a good check object for gamma spectroscopy – let’s draw the spectrum and see what’s 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 outcome seems to be like this:

In accordance with gammaspectacular.com, this spectrum belongs to Thorium-232. Clearly, it was not written within the product description, however with a gamma spectrum, we are able to see it with out object disassembly or any chemical checks!

5. Matplotlib Animation (Bonus)

As a bonus for readers affected person sufficient to learn up up to now, I’ll present make an animated graph in Matplotlib.

As a reminder, at first of the article, I collected spectra from a Radiacode system each minute. We are able to use this knowledge to animate a spectrum assortment course of.

First, let’s load the info 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 knowledge: a0, a1, a2, 1024 vitality values
   Instance:
   9.572;2.351;0.000362;0;2;0;0; ... """
   values = spectrum.break up(";")[3:]
   return np.array(record(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")
show(df)

The output seems to be like this:

As we are able to see, every row within the commentary incorporates the identical fields as we used earlier than for a single Spectrum object.

The animation course of is comparatively simple. First, we have to save a Matplotlib’s “bar” object within the variable:

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

Then we are able to use an replace perform that takes spectra from a dataframe:

def replace(body: int):
   """ Body replace callback """
   index = body * df.form[0] // num_frames
   counts = df["spectrum"].values[index]
   # Replace bar
   for i, b in enumerate(bar):
       b.set_height(counts[x1:x2][i])
   # Replace title 
   ax.set_title(...)

Now, we are able to save the animation into the GIF file:

import matplotlib.animation as animation

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

The output seems to be like this:

Right here, we are able to see how the variety of particles will increase throughout the assortment time.

6. Conclusion

On this article, I described the fundamentals of gamma spectroscopy evaluation with a Radiacode scintillation detector, Python, and Matplotlib. As we are able to see, the info itself is easy, nevertheless, loads of invisible processes are going “below the hood.” Utilizing completely different check objects, we had been capable of see spectra of various parts of radioactive decay as Bismuth-214 or Protactinium-234. And, amazingly, these checks could be executed at residence, and detectors like Radiacode can be found for science fanatics for the worth of a mid-level smartphone.

This fundamental evaluation permits us to go additional. If there may be some public curiosity on this subject, within the subsequent article I’ll use machine studying to automate the isotopes detection. Keep tuned.

For these readers who wish to carry out the checks on their very own, all datasets are available on Kaggle. If somebody wish to do the checks with their very own {hardware}, a Radiacode producer kindly supplied a ten% low cost code “DE10” for buying the device (disclaimer: I don’t have any revenue or different industrial curiosity from these gross sales).

All readers are additionally welcome to attach by way of LinkedIn, the place I periodically publish smaller posts that aren’t sufficiently big for a full article. If you wish to get the complete supply code for this and different articles, be at liberty to go to my Patreon page.

Thanks for studying.