embedded systems and applications in physics

miree

Introduction

People have built DIY gamma spectrometers and information can be found in the internet:

http://physicsopenlab.org/2016/01/26/diy-gamma-spectrometry/

http://www.gammaspectacular.com/pra-spectrometry-software

http://www.theremino.com/en/blog/gamma-spectrometry

http://www.gammaspectacular.com/index.php?route=common/home

All these projects use the sound card of a PC as data acquisition system. The "gammaspectacular" is a small device that integrates a high voltage source and some signal conditioning to make the signals suitable for sound card acquisition.

Theremino and PRA are software tools to extract signal amplitudes from the recorded sound data and create a pulse-height histogram.

However, what I didn't see yet is a setup or software support for gamma coincidence spectroscopy. So I decided to make the software and hardware for this and try out if it is possible.

Update: I just saw this forum post: where it was done before.

http://www.gammaspectacular.com/phpBB3/viewtopic.php?f=5&t=274

Gamma radiation

All atoms have a positively charged core, consisting of protons an neutrons, called the atomic nucleus. Electromagnetic radiation emitted by an atomic nucleus is called gamma radiation.
This emission of gamma radiation can happen if a nucleus is in an excited state (e.g. it vibrates, rotates, or has protons/neutrons in higher orbits), and dissipates its excitation energy by releasing an electromagnetic wave quantum - a photon in the gamma energy range. This process is also known as gamma decay of the nucleus. The energy (inverse to its wavelength) of the released photon matches exactly the energy difference between energy of the nucleus before and after the decay. These excitation energies are very specific for each isotope and can therefore be used to identify gamma radioactive isotopes.
How does an atomic nucleus get into an excited state? In our environment it happens all the time if naturally occurring, unstable nuclei decay into a more stable isotope. There are many elements with unstable isotopes. One example is the isotope 40-K that has a natural abundance of 0.0117%. The number 40 refers to the sum of protons and neutrons inside the atomic nucleus. A sample of natural K (Potassium) contains 93.2581% 39-K (19 protons and 20 neutrons), 0.0117% 40-K (19 protons and 21 neutrons), and 6.7302% 41-K (19 protons and 22 neutrons). Only 40-K decays into 40-Ar (18 protons and 22 neutrons) with a half-life of 1.248 trillion years. (this data can be found here http://www.nndc.bnl.gov/)
during the decay of 40-K into 40-Ar energy is released. This energy goes partly into the kinetic energy of the beta particle that is emitted during the decay process, but partially into internal degrees of freedom of the 40-Ar nucleus.
The excited 40-Ar nucleus is this in an excited state and can emit the aforementioned gamma radiation. In this case it has an energy of 1.461 MeV (mega electron volts) which corresponds to a wavelength of 0.00085 nm (visible red light has a wavelength of approx. 650 nm), i.e. gamma radiation is much(!) higher in energy than visible light.

The gamma radiation can be detected in a Radiation detector, e.g. a scintillation detector or liquid nitrogen cooled high-purity Ge detector.
Potassium is a chemical component of concrete, thus gamma radiation from 40-Ar is all around us. But this is not the only energy component of natural radiation as can be seen in the following picture from an approx. 190 hours measurement inside the SPIRAL facility in France using HPGe detectors (I found it online some time ago, but cannot find the link anymore). It shows a pulse-height histogram of the detector signals, and the pulse-height is given in units of keV (kilo electron volts):

This gamma ray spectrum was measured with a high resolution high-purity Germanium detector (HPGe), so a lot of distinct energies can be seen. The most prominent lines come from the decay of 40-K and 208-Tl

Coincident ramma radiation

Some Isotopes decay into excited states which make multiple gamma emissions. One example for this is 176-Lu.

The excited state in 176-Hf decays into another excited state at lower energy, wich again decays into an excited state wich decays into the nuclear ground state. The half-life of all excited states is very(!) short. In the order of picoseconds or nanoseconds. Much shorter than a sound card acquisition system can resolve. For all practical purposes, these decays happen at the same time. If two detectors are used, this information can be used. the acquisition system only considers signals that appeared simultaneously in both detectors, there is a high chance that these two signals actually came from the sample that we want to investigate, and not from some background radiation source. It is very unlikely that a background 40-K decay happens at the exact same time together with the 176-Lu decay.

A common way of showing such data is a coincidence histogram, or coincidence matrix. It is a two dimensional histogram where the simultaneous detection of a pair of signal amplitudes is counted. Each pixel in the picture corresponds to a bin in the histogram, and the height of each bin is encoded in the color. The amplitude of the detectors is in arbitrary units.

Setup summary

The previous picture shows my first recorded gamma coincidence spectrum of some LYSO crystals using a DIY system at home. The sample source was a buch of LYSO scintillation crystals from Ebay. LYSO is a scintillation material that contains 176-Lu. As detectors I used tow 3x3'' NaI:Tl detector assemblies that I scored on Ebay over the years. For the high-voltage I use potted modules from Hamamatsu that I also got from Ebay. The signal shaping is done with high impedance MOSFET source followers that are directly integrated in the PMT sockets together with the voltage divider. From there, the signals are going directly into the microphone input of my PC. Histograms are produced with a small C++ program that uses the alsa API on Linux. My sound card has a stereo sampling rate of 192 kHz. It is essential that the sound card has true stereo recording, otherwise a coincidence setup is impossible.

Detectors

Gamma radiation can be detected using scintillation crystals, which absorb the gamma photon and convert the energy into a large number of photons in the visible spectral range. The quantity of visible light photons is proportional to the energy of the absorbed gamma photon. Typical photon numbers are a few 10000, which is still a very dim flash of light on a macroscopic scale.
A photomultiplier tube (PMT) can be used to detect these tiny amounts of light and convert it into a measurable electrical signal. Scintillation detectors have an energy resolution that is much worse than the one of HPGe detectors (in the room background picture above). But they are much easier to handle, because HPGe detectors need to be cooled to liquid nitrogen temperature to work. Typical resolution of a HPGe detector is ~0.2%, while a NaI:Tl detector can achieve ~7% - more than one order of magnitude worse.
Here are the detector modules that I use for my coincidence experiments. My PMT bases are open, and one has to be careful with the high voltage! I do not recommend to try this at home! If you make your own setup please make sure you know what you're doing. In any case you do this at your own risk!

On the picture the detectors are not connected to anything, but during a measurement it looks the same, only that some coaxial and HV cables are attached, and a radioactive sample is positioned between them. Because of the coincidence setup no shielding is needed to get clean spectra.

PMT base

The wiring of the PMT socket is shown in the following picture, where the numbers refer to the socket pin numbers (they are written on the socket). The resistor values are all identical. I chose 220k for one socket, and 920k for the other socket. They both work fine, so I guess that anything in that range is fine.

Here is a closeup of the source follower amplifier:

The source follower is very simple and also cheap (I would guess below 20 ct. component price), and sufficient for the purpose. Here is the schematic drawing

Note that there is nothing else connected to the PMT anode. The anode is discharged slowly to ground through the high impedance of the source follower input. The high impedance is needed to get the signals long enough to be captured by the relatively slow sampling rate of the sound card.

The voltage divider is a series of same value resistors... the actual value is not so important as long as your HV power supply can handle the current that is created. There are capacitors on the last three stages of the PMT. I didn't try if they are really needed or if it also works without them.

Software

For single detector setup it is the easiest to take one of the existing software tools for spectroscopy (see links in the introduction). These programs do not support coincidence setups, so I had to write my own.

My software is really simple. It continuously takes the digitized signal from the sound card, puts it through a moving average filter with a width of 10 samples. Then, the difference between the momentary filter signal and the delayed filter signal is calculated (delay is again 10 samples). The resulting signal has sharp dips and peaks. Each time a peak in height is reached, the last dip in height is subtracted and the result is taken as pulse height. If both sound card channels see an amplitude at the same time (same sample) the coincidence spectrum is incremented at the position given by the two amplitudes.
The signal processing is visualized (real data) in the following picture

It is very simple but works well enough. The program to read data from the soundcard, do the signal processing and create histograms can be seen below. It is a bit quick and dirty: global variables and not very modular structure and hard coded parameters (such as smaplingrate). I put some comments to make clear what happens. Feel free to copy, use, modify, improve, redistribute, ...
It can be compiled with this makefile:

CXXFLAGS = -O -std=c++11
LDFLAGS  = -pthread -lasound
all: soundcard_daq 
soundcard_daq: soundcard_daq.cpp

source code highlighting by http://hilite.me/

#include <iostream>
#include <fstream>
#include <vector>
#include <cstdint>
#include <vector>
#include <cstdint>
#include <pthread.h>
#include <cstdlib>
#include <ctime>
#include <alsa/asoundlib.h>

const int sampling_rate = 192000; // set sampling rate here (in Hz)
//const int sampling_rate =  44100; // set sampling rate here (in Hz)

// store the recorded trace for one second
std::vector<int64_t> trace(sampling_rate,0);
std::vector<int64_t> mwatrace(sampling_rate,0);
std::vector<int64_t> filtertrace(sampling_rate,0);
// histograms for channel 1 and 2, and for the coincidence matrix
std::vector<int32_t> histogram1(8192,0);
std::vector<int32_t> histogram2(8192,0);
std::vector<std::vector<int32_t> > matrix(400,std::vector<int32_t>(400,0));

// Signal detection classes 
class Mwa // moving average window
{
public:
 Mwa(int l) : buffer(l,0), idx(0), sum(0)
 {}
 int64_t filter(int64_t value)
 {
  sum -= buffer[idx];
  buffer[idx] = value;
  sum += value;
  ++idx;
  if (idx >= buffer.size()) idx = 0;
  return sum;
 }
private:
 std::vector<int64_t> buffer;
 int idx;
 int sum;
};
class SignalFilter
{
public:
 SignalFilter(int l, int d, bool showtrace) : mwa1(l), mwa2(l), delay(d,0), idx(0), last_filter_value(0), direction(0), top(0), bottom(0), show_trace(showtrace)
 {}
 int64_t filter(int64_t value)
 {
  // we are recording with 32 bit, so first shift the value down 
  value >>= 16;
  static uint64_t sample_nr = 0;
  int64_t delayed_value = delay[idx];
  delay[idx] = value;
  ++idx;
  if (idx >= delay.size()) idx = 0;

  int64_t mwa1f = mwa1.filter(value);
  int64_t mwa2f = mwa2.filter(delayed_value);
  int64_t filter_value = mwa1f - mwa2f;

  // here we save the first second of the signal
  if (sample_nr < sampling_rate && show_trace)
  {
   trace[sample_nr]       = value;
   mwatrace[sample_nr]    = mwa1f;
   filtertrace[sample_nr] = filter_value;
   ++sample_nr;
  }

  // here is the peak detection and determination of the amplitude
  //   (amplitude = 0 means that no peak was detected)
  int64_t amplitude = 0;
  if (direction == 1) // going up before
  { 
   if (last_filter_value > filter_value) // we are suddenly going down
   {
    top = last_filter_value;
    amplitude = top-bottom;
   }
  }
  if (direction == -1) // going down before
  {
   if (last_filter_value < filter_value) // we are suddenly going up
   {
    bottom = last_filter_value;
   }
  }
  // remember if the signal was going up or down in this step
  if (last_filter_value < filter_value) direction =  1; 
  if (last_filter_value > filter_value) direction = -1; 

  last_filter_value = filter_value;
  return amplitude;
 }
private:
 Mwa mwa1, mwa2;
 std::vector<int64_t> delay;
 int idx;

 int64_t last_filter_value;
 int direction;
 int64_t top, bottom;

 bool show_trace;
};

// This function runs in a never ending thread. It receives the sound data
// in chunks of 32 samples per channel, sendst the data to the signal filters
// and determines peas amplitudes. If a peak amplitude was detected, the 
// value is put into the corresponding histogram.
void *pcm_thread(void* ptr)
{
 std::cerr << "start of data acquisition thread" << std::endl;
 snd_pcm_t               *handle;               // a handle to pcm device
 const snd_pcm_uframes_t  period_length = 32;    // 32 frames in one period
 const unsigned int       rate          = sampling_rate;
 const unsigned int       num_channels  = 2;     // stereo
 // the following buffer used to read out one period
 int32_t                  buffer[num_channels*period_length]; 
 SignalFilter             filter_ch1(10,10, true);
 SignalFilter             filter_ch2(10,10, false);
 int rc;

 // open pcm hardware
 rc = snd_pcm_open(&handle, "default", SND_PCM_STREAM_CAPTURE, 0);
 if (rc < 0) 
 {
  std::cerr << "unable to open pcm device: " << snd_strerror(rc) << "\n";
  return 0;
 }

 // get set of default hardware parameters
 snd_pcm_hw_params_t *params;
 snd_pcm_hw_params_alloca(&params);
 snd_pcm_hw_params_any(handle, params);
 
 // change hardware parameters:
 // Interleaved mode: right left right left right left right left ....
 rc = snd_pcm_hw_params_set_access(handle, params, SND_PCM_ACCESS_RW_INTERLEAVED);
 if (rc != 0) std::cerr << "error: snd_pcm_hw_params_set_access" << std::endl;
 // signed 32-bit little-endian format
 rc = snd_pcm_hw_params_set_format(handle, params, SND_PCM_FORMAT_S32_LE);
 if (rc != 0) std::cerr << "error: snd_pcm_hw_params_set_format" << std::endl;
 // two channels (stereo) 
 rc = snd_pcm_hw_params_set_channels(handle, params, 2);
 if (rc != 0) std::cerr << "error: snd_pcm_hw_params_set_channels" << std::endl;
 // set sample rate
 rc = snd_pcm_hw_params_set_rate(handle, params, rate, 0);
 if (rc != 0) std::cerr << "error: snd_pcm_hw_params_set_rate" << std::endl;
 // set period size 
 rc = snd_pcm_hw_params_set_period_size(handle, params, period_length, 0);
 if (rc != 0) std::cerr << "error: snd_pcm_hw_params_set_period_size" << std::endl;

 // write the parameters to the driver 
 rc = snd_pcm_hw_params(handle, params);
 if (rc < 0) 
 {
  std::cerr << "unable to set hw parameters: " << snd_strerror(rc) << "\n";
  return 0;
 }
 
 uint64_t time = 0;
 int samples = 0;
 int counts1 = 0;
 int counts2 = 0;
 int counts12 = 0;

 for(int i = 0;;++i)
 {
  rc = snd_pcm_readi(handle, (char*)buffer, period_length);
  if (rc < 0)
  {
   rc = snd_pcm_recover(handle, rc, 0);
   if (rc < 0)
    exit(1);
  }

  for (int i = 0; i < period_length; ++i)
  {
   ++samples;   
   ++time;
   uint64_t amplitude1 = filter_ch1.filter(buffer[2*i]);
   uint64_t amplitude2 = filter_ch2.filter(buffer[2*i+1]);

   // ratemeter output
   if (samples == sampling_rate)
   {
    std::cout << "ch1-rate | ch2-rate | coincidence-rate : " 
        << counts1 << " | " << counts2 << " | " 
        << counts12 << std::endl;
    samples  = 0;
    counts1  = 0;    
    counts2  = 0;
    counts12 = 0;
   }

   if (amplitude1) // fill histogram1
   { 
    // downscaling of amplitude to better fit in the histogram
    int ampl = amplitude1 / 100; 
    if (ampl > 0 && ampl < histogram1.size())
    {
     ++histogram1[ampl];
    }
    // ratemeter for channel1 (only amplitudes greater 60 are counted)
    if (ampl > 60)
    {
     ++counts1;
    }
   }
   if (amplitude2) // fill histogram2
   { 
    // downscaling of amplitude to better fit in the histogram
    int ampl = amplitude2 / 100;
    if (ampl > 0 && ampl < histogram2.size())
    {
     ++histogram2[ampl];
    }
    // ratemeter for channel1 (only amplitudes greater 60 are counted)
    if (ampl > 60)
    {
     ++counts2;
    }
   }
   if (amplitude1 && amplitude2) // fill coincidence matrix
   {
    // downscaling of amplitude to better fit into the matrix
    int ampl1 = amplitude1 / 1500;
    int ampl2 = amplitude2 / 1500;
    if (ampl1 < matrix.size() && ampl2 < matrix[0].size())
    {
     ++matrix[ampl1][ampl2];
    }
    if (ampl1 > 6 && ampl2 > 6)
    {
     ++counts12;
    }
   }
  }
 }
 std::cerr << "end of thread" << std::endl;
 return 0;
}

int delay_ms(unsigned long millis)
{
 struct timespec ts;
 int err;
 ts.tv_sec = millis / 1000;
 ts.tv_nsec = (millis % 1000) * 1000000L;
 err = nanosleep(&ts, (struct timespec *)NULL);
 return (err);
}

int main()
{
 // Start a never ending thread where the 
 //  recording and signal filtering is done.
 pthread_t recording_thread;
 if (pthread_create(&recording_thread, 0, &pcm_thread, 0))
 {
  std::cerr << "cannot create thread" << std::endl;
 }

 // ouptut updated histograms and coincidence matrix every 10 seconds
 for (int i = 0;;++i)
 {
  delay_ms(10000);
  std::cerr << "updating histograms"<< std::endl;
  std::ofstream histfile("histogram.dat");
  for (int i = 0; i < histogram1.size(); ++i)
   histfile << i << "\t" << ((i>10)?histogram1[i]:0) 
        << "\t" << ((i>10)?histogram2[i]:0)
        << std::endl;
  histfile.close();

  std::ofstream matrixfile("matrix.dat");
  for (int i = 0; i < matrix.size(); ++i)
  {
   for (int j = 0; j < matrix[i].size(); ++j)
   {
    matrixfile << matrix[i][j] <<  " ";
   }
   matrixfile << std::endl;
  }
  matrixfile.close();

  // the signal of channel 1 in the first second of recording
  if (i == 0)
  {
   std::ofstream tracefile("traces.dat");
   for (int i = 0; i < trace.size(); ++i)
    tracefile << i << " " << trace[i] 
            << " " << mwatrace[i] 
            << " " << filtertrace[i] << std::endl;
   tracefile.close();
  }
 }
 return 0;
}