Message Boards Message Boards

3 Replies
30 Total Likes
View groups...
Share this post:

Detect radioactivity at home with Mathematica

I thought I would share the following rather simple experiment. It is based on a Geiger counter which interacts via the serial port with Mathematica. It can be used to measure ambient levels of radioactivity and to test different radioactive sources at home. As an additional benefit it generates "perfect/physical" random numbers (as opposed to pseudorandom numbers), some of which I attach to this post. 

The setup is rather simple, unfortunately one of the parts is a bit expensive. I use

1) Sparkfun's Geiger Counter SEN-11345. (about 150 $ and similar in £)
2) A USB 2.0 A/Mini-B cable.

The Geiger tube is on the bottem right on the photo. The parts inside the dashed line are under high voltage when the thing is switched on; avoid touching any of these parts. The Geiger counter has an onboard AMTEL chip (like Arduino) so it communicate easily via the serial port. You might have to install the appropriate VCP driver from this website.

Once the driver is installed you connect the Geiger counter to the USB port of the computer. On the software side you need to install the SerialIO package. Step by step instructions for the installation can be found on this website.

The Sparkfun Geiger counter sends a sequence of 0 and 1 to the serial port, which indicate whether the previous interval is longer or shorter than the current interval between counts. Each number corresponds to one event - the counter measures alpha, beta and gamma radiation; the latter only if the red protective cap is removed. The sequence of zeros and ones can be used as a set of independent random numbers. For most simulations the pseudorandom numbers that computers (and Mathematica) generate, but for some other applications such as encryption physical random numbers are preferable. The code below will produce that string of random digits, but it will also record the times of the measured decay events. The code is for Mac OSX but can easily be adapted for other operating systems.

 (*Clear everything*)
 (*Load SerialIO*)
 << SerialIO`
 (*Connect to the Geiger counter*)
 myGeigerCounter =
 SerialOpen[Quiet[FileNames["tty.usb*", {"/dev"}, Infinity]][[1]]];
SerialSetOptions[myGeigerCounter, "BaudRate" -> 9600];
While[SerialReadyQ[myGeigerCounter] == False, Pause[0.1]];

s = {}; tZero = AbsoluteTime[]; SerialRead[myGeigerCounter];
AppendTo[s, {SerialRead[myGeigerCounter], AbsoluteTime[] - tZero}];,
UpdateInterval -> 0.002],
SynchronousUpdating -> False]

That is basically the entire code for reading the data. The critical point is the UpdateInterval. If the computer is busy, it might measure too slowly and measure a sequence of 0s and 1s instead of a single digit, so it is important to have few background processes running. The UpdateInterval of 0.002 has shown to be sufficiently high for my experiments; Sparkfun's detector can measure up to 100 Hz. Also, if there is no event for 10 seconds the program will read " " and record the time, i.e. last event plus 10 seconds. Therefore we need to clean up the data after the measurement is finished (simply abort the evaluation), see below. If you want to plot the measurements dynamically you can use: 

Dynamic[sc = DeleteCases[s, {"", ___}];
ListPlot[sc[[2 ;;, 2]] - sc[[1 ;; -2, 2]],
  PlotRange -> {All, {0, Max[sc[2 ;;, 2]] - sc[[1 ;; -2, 2]]}},
  ImageSize -> Large]]

This gives the following animation, where every couple of seconds a new point is added.

If you are interested in an estimate of the "counts per minute" this command will do the trick:

Dynamic[10/(sc[[-1, 2]] - sc[[-10, 2]])*60. // N]

So by now we have a working Geiger counter and we can visualise the results. 

To test the counter I ran a measurement of the background radiation for a couple of hours. I import the data:

tstest = Import["~/Desktop/Geigerdata_overnight.csv"];

 Then clean it up

tsBackground = DeleteCases[# & /@ tstest, {"", ___}] // ToExpression

and plot a diagram of the time between two consecutive events

ListPlot[tsBackground[[2 ;; -1, 2]] - tsBackground[[1 ;; -2, 2]],
PlotRange -> All, AxesLabel -> {"Decay Incidence", "Waiting Time"},
ImageSize -> Large]

There are sligthly more than 12000 events in this time series.
To compare that to some radioactive object I dismantled a smoke detector and extracted the radioactive Americanum 241 container. The sourse emits alpha and gamma particles and is supposed to have a rate of about 33000 Bq, i.e. decays per second. (It says 0.9uCi on the detector. Wolfram Alpha knows that that is 33300 Bq.)


I position the Americanum at about 3 cm from the tip of the Geiger counter.

I then measure the counts over a couple of hours and obtain:

This time there are more then 20000 events recorded. Note that the waiting times are much shorter than in the background case. We can now compare the histograms of the waiting time between two consecutive events:

The total number of events are not the same, but it is obvious that the waiting times are much shorter for the Americanum element than for the background (blue).
The waiting time distribution is expected to be exponentially distributed, so we can fit 
distBackground =
  tsBackground[[2 ;; -1, 2]] - tsBackground[[1 ;; -2, 2]],
which gives
distAmericanum =
  tsAmericanum[[2 ;; -1, 2]] - tsAmericanum[[1 ;; -2, 2]],
which gives

The mean waiting time for the background is

and for the Americanum element

Here is another representation
         tsBackground[[2 ;; -1, 2]] - tsBackground[[1 ;; -2, 2]]][[
        1]][[;; -2]],
          tsBackground[[2 ;; -1, 2]] - tsBackground[[1 ;; -2, 2]]][[
         2]] // Log // N}], Mesh -> All,
    AxesLabel -> {"waiting time", "Log number"}, ImageSize -> Large],
        tsAmericanum[[2 ;; -1, 2]] - tsAmericanum[[1 ;; -2, 2]]][[
       1]][[;; -2]],
         tsAmericanum[[2 ;; -1, 2]] - tsAmericanum[[1 ;; -2, 2]]][[
        2]] // Log // N}], PlotStyle -> Red, Mesh -> All]}]

Wolfram Alpha knows a lot of useful stuff about Americanum 241, try:
== AM 241

One little issue with the measurement is that the Americanum element is supposed to have 33000 Bq, i.e. decays per second. Our measurement, corrected for background radiation, shows that there is more or less one additional decay per second. When the measurement was conducted, the red cap at the tip of the Geiger counter was not removed; this filters out the alpha radiation. I read somewhere that only less than 10% of the events come from gamma radiation, even though that seems incorrect given the information on wikipedia. Also the detector was at a distance of about 3 cm from the element. Its opening is about 1cm^2. That might account for another factor of 100, but we would still be much below the at least 33 events per second. I have not really looked into this or done any calculations, but it might be nice to find the missing events.

This does not seem to be a problem with the 0.002 Updateinterval, as visual inspection of the time series shows that there are nearly no incidence of two bits being read at the same time from the serial port. The Geiger conter can measure up to 100Hz. The counts per minute of the background radiation are about 20 which is reasonable. (The measurement was made in Aberdeen (UK) where many houses are made from Granite which is slightly radioactive. Inside the hoses Radon is produced which can sometimes yield very high levels of radiation. This particular measurement was not taken in a Granite house.) 

If anyone cares about good random numbers, they can use

tsAmericanum[[All, 1]]

To extract a string of ones ans zeros, which are "real" random numbers. I attach a file with the sequence of ones and zeros I got from the Americanum experiment. You can also use that to analyse how many times our update interval was too slow so that two digits were transmitted at once. The file can be read with:
data = Import["~/Desktop/RandomAmericanum.csv"] // Flatten

A simple Tally shows that 

(*{{0, 11779}, {1, 11979}, {11, 65}, {10, 1}}*)
Only 66 cases out of about 22000 were double digits.

Much about the code can/should be improved. Using an intermediate Arduino might improve the measurements slightly.


PS: Note that there might be regulations in your country regarding the Americanum element of the smoke detector, e.g. how to dispose of them. Also Sparkfun does not deliver the Geiger Counter to all countries. 
POSTED BY: Marco Thiel
3 Replies

Here's a short update with an implementation in Mathematica 10. Due to the new DeviceOpen function SerialIO is not required anymore.

 (*Connect to the Geiger counter*)
myGeigerCounter = 
DeviceOpen["Serial", {Quiet[FileNames["tty.usb*", {"/dev"}, Infinity]][[1]], "BaudRate" -> 9600}]

s = {}; tZero = AbsoluteTime[]; DeviceRead[myGeigerCounter];
RunScheduledTask[AppendTo[s, {Quiet[ToExpression[FromCharacterCode[DeviceRead[myGeigerCounter]]]], 
AbsoluteTime[] - tZero}], 0.01]

This is obviously shorter and more elegant than before. Note that the update interval is 0.01. In the Mathematica 9/ SerialIO implementation the update interval was 0.002, i.e. a fifth of that. I noticed that such a short update interval causes problems with RunScheduledTask. Cleaning the data is the same.

sc = Cases[s, {_Integer, ___}];

The histogram can be displayed like this:

Histogram[Differences[sc[[All, 2]]], {0, 20, 1}, ScalingFunctions -> "Log", LabelingFunction -> Top, ImageSize -> 500]

enter image description here

I will try to update the post and add some more analysis later.


POSTED BY: Marco Thiel
Dear Arnoud,

that is fantastic! I cannot wait to lay hands upon Mathematica 10. I do have a larger number of project which will greatly benefit from Mathematica 10's new features. I had already planned to "translate" all of these programs, including the Geiger Counter to Mathematica 10. Thanks for your great input on that. I very much like that cloud export feature. That will make many new projects possible. The following post is a bit long, but the main question is:

How fast can Mathematica 10 read serial data when connected to an Arduino?

As I said, unfortunately, I do not (yet) have access to Mathematica 10 and that's why I have to run lots of things using SerialO. One critical thing is the rate with which Mathematica can read the serial port. This is crucial for this particular project but also for my thermal imaging camera. In the latter case I can reduce the scanning time if Mathematica reads quick enough. There is also the limit of the actual sensor though. 

In the radioactivity example, it is indeed quite possible to monitor radioactivity over a rather long time period automatically and with Mathematica 10 it seems to be really quite straight forward to use the cloud and monitor everything from a distance. That is brilliant.

In my post I also tried to read in the actual event times. The Geiger counter does not provide the times, so unless I want to reprogram the chip on the Geiger Counter, I need to read it into Mathematica as quick as I can and time-stamp it. Now that crucially depend on the speed of the SerialRead and was why I tried this in the first place. The waiting times between two consecutive events is exponentially distributed and in the case of the Americanum source Mathematica fitted 
The thing is that sometimes events will occur in very quick succession. The PDF is
PDF[ExponentialDistribution[1.15882], x]
The sampling rate in my example was 0.002 in Mathematica with SerialIO.
s = {}; tZero = AbsoluteTime[]; SerialRead[myGeigerCounter];
AppendTo[s, {SerialRead[myGeigerCounter], AbsoluteTime[] - tZero}];,
UpdateInterval -> 0.002],
SynchronousUpdating -> False]

I noticed that the Raspberry Pi has problems with that because it is often busy with other things and cannot maintain that update interval. My laptop (an old one for experiments) has no problems reading at that speed if it doesn't do anything else. Now here comes the problem. How likely is it that two events occur faster than the 0.002 seconds?
NIntegrate[PDF[ExponentialDistribution[1.15882], x], {x, 0, 0.002}]


If I have 22000 events I would naively expect 


So slightly more than 50 events which are too fast for our data aquisition. In my original post I observed 66 which is ok I guess. I am not sure yet, but the some of the discrepancy might be due to the actual refresh rate being slightly lower than 0.002 seconds; and also because I effectively have a >0.002 binning. 

I think that one can make an Arduino sample analogue data at about 8000 times per second - see this post. This limit is bascially from the AD converter; I am currently trying to link other AD converters to Mathematica to push the limit a bit further. It would be nice to get Mathematica to read this as quick a possible.

I am really curious to see how Mathematica 10 deals with this. Could you by any chance try and see how fast Mathematica 10 can read serial reliably? That would greatly help me, because I could prepare some other projects so that they work the day Mathematica 10 comes out. 

I think that in your example you read the entire device buffer, which of course, makes me see all 0s and 1s as long as I read faster than the buffer fills up. In my example I tried to read only one bit at a time and use Mathematica's time stamp to record the event times. So if a sensor sends data quickly and does not provide times, how well can Mathematica 10 deal with that? (There is a more serious application I have in mind for this...)

Thank you very much for your post. I am very much looking forward to trying it myself on MMA10.


PS: I'll post something on the speed of data aquisition a bit later. 
POSTED BY: Marco Thiel
Very cool post, thanks!

I am working on reproducing your example in the upcoming Mathematica 10 release, using the new Device API functions, which replaces (among other things) the SerialIO package.

Here is how I connected to the device:
serial = DeviceOpen["Serial", {"COM4", "BaudRate" -> 9600}]
In my case I work on a Windows machine and the device presented itself on COM4 after being plugged in. I think I already had the appropriate VCP driver installed (for another project). The first argument to DeviceOpen is "Serial" which, unsurprisingly, opens a serial connection. The baud rate is specified as an option to the "COM4" device.

Next I set up a serial read function:
read[serial_] := Module[{data, count},
data = DeviceReadBuffer[serial];
If[Length[data] > 0, count = Total[ToExpression /@ FromCharacterCode /@ data], count = 0];
{AbsoluteTime[], count}
Here, DeviceReadBuffer will read all available data that is currently in the serial buffer. I use FromCharacterCode to translate the character code '48' to the string "0" and '49' to the string "1". Next I use ToExpression to turn "0" and "1" into the numbers 0 and 1. I use Total to aggregate the number of hits in a given time period. If no data was read (empty list), I set the count to zero as well. I return the absolute time and the count.

Here I set up a scheduled task, which runs every 60 seconds:
data = {};

RemoveScheduledTask /@ ScheduledTasks[];
RunScheduledTask[AppendTo[data, read[serial]], 60];
I start out with data being an empty list. I then clear out any previously scheduled tasks. Then I run the scheduled task every 60 seconds and append the new data point to the list of data points.

Finally, I visualize the data by using DateListPlot:
Dynamic[DateListPlot[data, Joined -> True, Filling -> Axis],
UpdateInterval -> 60]

Here is an example of the data collected in about 60 minutes:

And a longer time period (12+ hours):
POSTED BY: Arnoud Buzing
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract