Exploratory Data Analysis: Gamma Spectroscopy in Python

, knowledge science and knowledge evaluation may be carefully associated to physics and {hardware}. Not every thing can run within the cloud, and a few functions require using actual issues. On this article, I’ll present tips on how to accumulate the info from a Radiacode 103G radiation detector and gamma spectrometer, and we are going to see what sort of info we will get from it. We are going to do that evaluation in Python, and I can even 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 mechanically.

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 discuss concerning the radiation. Which is, curiously, all the time round us, and the radiation degree is rarely zero. Elevated ranges of radiation may 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 in the course of the flight is about 10x greater 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 crammed with a combination of gases. When a charged particle reaches the tube, the fuel is ionized, and we will detect the brief pulse. The upper the radiation degree, the extra pulses per minute we get. Radiation detectors usually present values in CPM (counts per minute), which may 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 kind of crystal, which generates mild when a charged particle is detected. A scintillator has an vital property—the depth of the sunshine is proportional to the particle’s vitality. Due to that, we cannot solely detect the particles however can even decide which sort of particles we get.

Clearly, from a {hardware} perspective, it’s simpler mentioned than carried out. In a scintillation crystal, solely a number of photons may 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 will purchase a scintillation detector for the value of a mid-level smartphone. This makes gamma spectroscopy evaluation doable even for science fans with a average funds.

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

2. Amassing 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 will get each CPS (counts per second) and gamma spectrum values. Radiacode wants solely the USB connection and has the dimensions of a USB stick:

To get the info, I’ll use a radiacode open-source library. As a easy instance, let’s accumulate 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 will 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 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 """
    knowledge = rc.data_buf()
    for document in knowledge:
        if isinstance(document, RawData):
            dt_str = document.dt.strftime(self.TIME_FORMAT)
            log_str = f"{dt_str},,{int(document.count_rate)},"
            logging.debug(log_str)

def process_spectrum_data(rc: RadiaCode):
    """ Get spectrum knowledge from the gadget """
    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 sorts 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 gadget works; they aren’t used for evaluation. A Spectrum knowledge incorporates what we’re focused on – gamma spectrum values.

Amassing the info each 60 seconds permits me to see the dynamics of how the info is altering. Normally, the spectrum should 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. Knowledge Evaluation

3.1 Gamma Spectrum

To grasp what sort of knowledge we now 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 offers us not solely the variety of particles but in addition 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 may be discovered utilizing a easy formulation:

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

In the course of 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 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 appears to be like like this:

I collected the info in 30 seconds, and even inside this brief interval, we will already see a likelihood distribution! As common in statistics, the longer the higher, and by gathering the info inside a number of hours, we will get good-visible outcomes. It is usually handy that Radiacode collects the spectrum mechanically. We will preserve 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 capabilities 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 will use it to get the info:

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 objective of the gamma spectrum? It turned out that completely different parts emit gamma rays with completely different energies. That permits us to inform what sort of object we now 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 will 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 degree but in addition see why the article is radioactive and the way the decay course of goes.

How does it work in follow? Let’s say I need to measure radiation from a banana. Bananas naturally comprise Potassium-40, which emits gamma rays with a 1.46 MeV vitality. In concept, if I take numerous bananas (actually lots 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 usually extra difficult. Some supplies, like uranium-238, might have a posh decay chain like this:

As we will see, uranium has a collection of radioactive decays. It decays to thorium, thorium to radium, and so forth. Each aspect has its gamma spectrum peak and a half-life time. Because of this, all these peaks can be current on a spectrum on the identical time (however with completely different intensities).

We will 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 traces of code in Matplotlib to attract it on a spectrum:

# Isotopes
for identify, en_kev, shade in isotopes:
   plt.axvline(x = en_kev, shade = shade, 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 take a look at sources like Cesium-137 or Strontium-90, utilized in a “large science.” Nevertheless, it isn’t required, and a few objects round us may 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 scale back the background radiation to a a lot decrease degree. Nevertheless, lead is heavy, the cargo may be costly, and it’s a poisonous materials that requires numerous security precautions. As a substitute, I’ll accumulate 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 can be seen.

Let’s mix all elements 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, shade in isotopes:
      plt.axvline(x = en_kev, shade = shade, 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("Power, keV")
   plt.legend()
   fig.tight_layout()
   fig.present()

Completely different spectra may be collected throughout completely different time intervals. Due to that, I normalised the graph by dividing the rely values by the full 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 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 troublesome object to measure – it incorporates lower than 1g of potassium-40, and it additionally incorporates 74% of water. As we will 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 might see a small distinction.

Utilizing the strategies created earlier than, it’s straightforward 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 appears to be like like this:

As we will 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 full variety of particles, brought on by background radiation. As we will see, the banana is just 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 palms and dials from the 1910s to the Nineteen Sixties. The combination of radium and a particular lume allowed watches to glow at the hours of darkness. 
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 indicate how the watch was glowing at the hours of darkness 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 appears to be like like this:

A lot of processes are occurring right here, and we will see the elements 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 a substitute, we will see some daughter merchandise as bismuth and lead.

As an apart query: is it harmful to have a watch with a radium dial? Usually not, these watches are protected if the case and the glass usually are not broken. We will see many decay merchandise on the spectrum, however in the identical approach 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 for the reason that nineteenth century, lengthy earlier than even the phrase “radioactivity” turned recognized. Uranium glass was produced in big quantities and may now be present in virtually each classic retailer. The manufacturing stopped within the Nineteen 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 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 degree from the uranium glass isn’t excessive. Many of the peaks are small, and I used a log scale to see them higher. The outcome appears to be like like this:

Right here, we will 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 folks didn’t find out about radiation. Humorous sufficient, it isn’t all the time true, and plenty of “Detrimental Ion” associated merchandise are nonetheless obtainable as we speak.

I purchased this pendant for $10 on eBay, and in response to the vendor’s description, it may be used as a damaging 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 outcome appears to be like 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 will 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 so far, I’ll present tips on how to make an animated graph in Matplotlib.

As a reminder, at the start of the article, I collected spectra from a Radiacode gadget each minute. We will 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 appears to be like like this:

As we will see, every row within the remark incorporates 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 will use an replace operate 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 will 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 appears to be like like this:

Right here, we will 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 will see, the info itself is easy, nevertheless, numerous invisible processes are going “below the hood.” Utilizing completely different take a look at objects, we had been in a position to see spectra of various parts of radioactive decay as Bismuth-214 or Protactinium-234. And, amazingly, these assessments may be carried out at dwelling, and detectors like Radiacode can be found for science fans for the value 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 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 offered 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 by way of LinkedIn, the place I periodically publish smaller posts that aren’t large enough 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.