Introduction
People have built DIY gamma spectrometers and information can be found in the internet:
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.
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:
1
2
3
4
| CXXFLAGS = -O -std=c++11
LDFLAGS = -pthread -lasound
all: soundcard_daq
soundcard_daq: soundcard_daq.cpp
|
source code highlighting by
http://hilite.me/
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
| #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(¶ms);
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;
}
|
