Exploratory Data Analysis: Gamma Spectroscopy in Python

, information science and information evaluation will be intently associated to physics and {hardware}. Not the whole lot can run within the cloud, and a few purposes require the usage of actual issues. On this article, I’ll present how one can gather the information from a Radiacode 103G radiation detector and gamma spectrometer, and we are going to see what sort of data we are able to get from it. We’ll 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 varieties and isotopes robotically.

All recordsdata collected for this text can be found on Kaggle, and readers are additionally welcome to run all assessments on their very own. The hyperlink is added to the tip 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, apparently, at all times round us, and the radiation stage is rarely zero. Elevated ranges of radiation will be discovered in lots of locations and objects, from classic watches within the thrift retailer to airplane flights (due to cosmic rays, the radiation stage in the course of the flight is about 10x greater in comparison with floor stage).

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

• A Geiger-Müller tube. As its title suggests, it’s a tube stuffed with a mix of gases. When a charged particle reaches the tube, the fuel is ionized, and we are able to detect the brief pulse. The upper the radiation stage, the extra pulses per minute we get. Radiation detectors typically present values in CPM (counts per minute), which will be transformed into Sieverts or different items. The Geiger tube is affordable and dependable; it’s utilized in many radiation detectors.
• A scintillation detector. This detector is predicated on a particular sort of crystal, which generates mild when a charged particle is detected. A scintillator has an essential property—the depth of the sunshine is proportional to the particle’s vitality. Due to that, we cannot solely detect the particles however may decide which kind of particles we get.

Clearly, from a {hardware} perspective, it’s simpler mentioned than achieved. In a scintillation crystal, solely a number of photons will be emitted when a particle is detected. Earlier, these detectors had a 5-6 digit value, 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 reasonable price range.

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

2. Accumulating the Knowledge

As was talked about at first, I’ll use a Radiacode 103G. It’s a conveyable gadget 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 information, I’ll use a radiacode open-source library. As a easy instance, let’s gather the gamma spectrum inside 30 seconds:

from radiacode import RadiaCode


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

spectrum = rc.spectrum()
print(spectrum.period)
#> 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 which means of those fields within the following chapter after 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 subsequent 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 information from the gadget """
   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 """
    information = rc.data_buf()
    for report in information:
        if isinstance(report, RawData):
            dt_str = report.dt.strftime(self.TIME_FORMAT)
            log_str = f"{dt_str},,{int(report.count_rate)},"
            logging.debug(log_str)

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

def save_spectrum_data(spectrum: Spectrum):
    """ Save spectrum information 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 sorts of information. A RawData comprises CPS values, which can be utilized to get radiation ranges in µSv/hour. I exploit them solely to see that the gadget works; they don’t seem to be used for evaluation. A Spectrum information comprises what we’re taken with – gamma spectrum values.

Accumulating the information each 60 seconds permits me to see the dynamics of how the information is altering. Often, the spectrum should be collected inside a number of hours (the extra the higher), and I favor 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. Knowledge Evaluation

3.1 Gamma Spectrum

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

spectrum = rc.spectrum()
print(spectrum.period)
#> 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 offers 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 will be discovered utilizing a easy method:

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

Through the specified time (in our case, 30 seconds), the Radiacode detects the particles and saves the end result into channels as a easy arithmetic sum. For instance, the worth “2” within the fifth place signifies 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("Power, keV")
   plt.legend()
   fig.tight_layout()
   fig.present()

draw_simple_spectrum(spectrum)

The output seems to be like this:

I collected the information in 30 seconds, and even inside this brief interval, we are able to already see a chance distribution! As common in statistics, the longer the higher, and by gathering the information inside a number of hours, we are able to get good-visible outcomes. Additionally it is handy that Radiacode collects the spectrum robotically. We are able to maintain the gadget operating 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.period.total_seconds()
   information = {
       "a0": spectrum.a0,
       "a1": spectrum.a1,
       "a2": spectrum.a2,
       "counts": spectrum.counts,
       "period": duration_sec,
   }
   with open(filename, "w") as f_out:
       json.dump(information, f_out, indent=4)
       print(f"File '{filename}' saved, period {duration_sec} sec")

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

We are able to use it to get the information:

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

As talked about at first, all collected recordsdata 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 components 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 article is radioactive and see that the extent is, let’s say, 500 counts per minute. With a gamma spectrum, we cannot solely see the radiation stage but additionally see why the article is radioactive and the way the decay course of goes.

How does it work in follow? 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 principle, if I take a number of bananas (actually loads 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 difficult. Some supplies, like uranium-238, could have a posh decay chain like this:

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

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

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

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

# Isotopes
for title, en_kev, shade in isotopes:
   plt.axvline(x = en_kev, shade = shade, label=title)
# 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 take a look at sources like Cesium-137 or Strontium-90, utilized in a “huge science.” Nonetheless, it isn’t required, and a few objects round us will be barely radioactive.

Ideally, a take a look at object should be positioned in a lead chamber to scale back the background radiation. A thick layer of lead will cut back the background radiation to a a lot decrease stage. Nonetheless, lead is heavy, the cargo will be costly, and it’s a poisonous materials that requires a number of security precautions. As an alternative, I’ll gather two spectra—the primary for the article itself and the second for the background (with out the article). As we are going to see, it’s sufficient, and the distinction might be seen.

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

def draw_spectrum(
   spectrum: Spectrum, background: Spectrum, isotopes: Listing, title: str
):
   """ Draw the spectrum, the background, and the isotopes """
   counts = np.array(spectrum.counts) / spectrum.period.total_seconds()
   counts_b = np.array(background.counts) / background.period.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 = listing(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 title, en_kev, shade in isotopes:
      plt.axvline(x = en_kev, shade = shade, label=title)
   # 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("Power, keV")
   plt.legend()
   fig.tight_layout()
   fig.present()

Completely different spectra will be collected throughout completely different time intervals. Due to that, I normalised the graph by dividing the depend values by the whole time, so we at all times see the counts per second on a graph.

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

4.1 Bananas

First, let’s reply probably the most requested query in nuclear physics – how radioactive are the bananas? A banana is a surprisingly tough object to measure – it comprises lower than 1g of potassium-40, and it additionally comprises 74% of water. As we are able to guess, the radiation from a banana could be very low. First, I attempted with a daily 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 end 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 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.period.total_seconds()
counts_b = np.array(background.counts) / background.period.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 whole variety of particles, brought on by background radiation. As we are able to see, the banana is simply 0.019% extra radioactive than the background! Clearly, this quantity is just 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 take a look at some “stronger” objects. Radium-226 was used for painting watch arms and dials from the 1910s to the Sixties. The combination of radium and a particular lume allowed watches to glow at nighttime. 
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 mild to point out how the watch was glowing at nighttime 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 info was collected inside 3.5 hours, and the output seems to be like this:

Plenty 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 fuel (3.82-day half-life), however the radon itself can also be an alpha emitter and doesn’t produce gamma rays. As an alternative, 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 protected 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 manner as with bananas, this spectrum was collected inside a number of hours, and the true radiation stage 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 for the reason that nineteenth century, lengthy earlier than even the phrase “radioactivity” turned identified. Uranium glass was produced in big quantities and may now be present in virtually 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 comprises solely a tiny fraction of uranium; these objects are protected to maintain within the cupboard (nonetheless, I’d 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 take a look at 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 stage from the uranium glass is just not excessive. A lot of the peaks are small, and I used a log scale to see them higher. The end result seems to be like this:

Right here, we are able to see some peaks from thorium and uranium; these components 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 folks didn’t learn about radiation. Humorous sufficient, it isn’t at all times true, and plenty of “Adverse Ion” associated merchandise are nonetheless obtainable right now.

I purchased this pendant for $10 on eBay, and in accordance with the vendor’s description, it may be used as a destructive 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 take a look at 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 end result seems to be like this:

In response to 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 assessments!

5. Matplotlib Animation (Bonus)

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

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

First, let’s load the information 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 information: a0, a1, a2, 1024 vitality values
   Instance:
   9.572;2.351;0.000362;0;2;0;0; ... """
   values = spectrum.cut up(";")[3:]
   return np.array(listing(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 comprises the identical fields as we used earlier than for a single Spectrum object.

The animation course of is comparatively easy. 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 in the course of 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 information itself is easy, nevertheless, a number of invisible processes are going “beneath the hood.” Utilizing completely different take a look at objects, we had been capable of see spectra of various components of radioactive decay as Bismuth-214 or Protactinium-234. And, amazingly, these assessments will be achieved at house, 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 matter, within the subsequent article I’ll use machine studying to automate the isotopes detection. Keep tuned.

For these readers who want to carry out the assessments on their very own, all datasets are available on Kaggle. If somebody want to do the assessments 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 business curiosity from these gross sales).

All readers are additionally welcome to attach through 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 happy to go to my Patreon page.

Thanks for studying.