The Timepix hybrid silicon pixel detector is an impressive bit of kit. Developed at CERN by the Medipix Collaboration, and used in the Large Hadron Collider (LHC), it can detect and visualise ionising radiation in real time. In this challenge, we will use code to count and identify different types of radiation emitted by a radioactive notebook at the Royal Institution of Great Britain.
Figure 1: A CERN@school MX-10 detector (left) mounted above one of William Crookes' notebooks at the Royal Institution of Great Britain. You can also see a Geiger counter to the right.
Ionising radiation, such as alphas (helium nuclei), betas (electrons and positrons) and gammas (photons), registers as different patterns of pixels in each frame of data. The pixel information is recorded in a text file with the x and y position of the hit on the sensor and the count value C of the hit, which is related to the energy of the radiation.
We have provided data from a frame recorded in one second from the above experimental setup. We can read this data into Python as follows:
In [1]:
# Import the data analysis tools we need from the Pandas toolkit.
import pandas as pd
## Our pixel data stored as a Panda DataFrame.
df000 = pd.read_table('./data/crookes/data000.txt', header=None, names=["x", "y", "C"])
# Show the first five rows of pixel data.
df000.head()
Out[1]:
The DataFrame object provides a convenient wrapper for the pixel data read from the raw text file the detector produced. For example, we can find out how many pixels there are in total in the frame:
In [2]:
print("There are %d pixels in the data frame." % (len(df000)))
Of course, what we're really interested in is what the pixels look like as a whole so we can see the patterns and start to identify what has been detected. We've written some code to do this for you:
In [3]:
# Import a function that extracts a pixel dictionary from the data frame.
from cernatschool.helpers import get_pixel_dictionary_from_dataframe
# Import a function that makes an image from a pixel dictionary.
from cernatschool.visualisation import make_frame_image
# (This line just allows us to show images as we go...)
%matplotlib inline
## The pixel dictionary - contains the pixel information in the format {X:C}.
pixel_dict = get_pixel_dictionary_from_dataframe(df000)
# Make the frame image from the pixel dictionary.
make_frame_image(pixel_dict)
And there you go! A 256 x 256 pixel image of the radiation detected in one second of recording data from above Crookes' notebook. The colours indicate how much energy was deposited in each pixel using the scale on the right of the image.
In [4]:
# Exercise 1: Manual particle counting
#
# How many particles can you count in the above image? Write your answer below.
#
#
Now, you should be able to make out three sorts of clusters - groups of pixels adjacent pixels - in the image, corresponding to the three types of radiation we'd expect to see from this measurement.
Firstly, there are the alphas. Alpha radiation consists of helium nuclei - two protons and two neutrons with an overall charge of 2+. These deposit large amounts of energy in the detector. To cut a long story short, this results in many pixels being activated over a large area in the image, resulting in "blobs" like these:
In [5]:
# Find the clusters - groups of adjacent pixels - in the data.
from cernatschool.kluster import KlusterFinder
# Get one of the alpha clusters.
alpha_cluster_1 = KlusterFinder(pixel_dict, 256, 256).getListOfKlusters()[4]
# Import a function for making an image of just one cluster.
from cernatschool.visualisation import make_kluster_image
# Make the image.
make_kluster_image(alpha_cluster_1)
Notice too how all of the energy is centred in the middle of the cluster.
Secondly, we have beta radiation - electrons and positrons - that scatters through the silicon of the detector leaving curly tracks of pixels like in the following cluster:
In [6]:
# An example of a beta cluster from the data.
beta_cluster_1 = KlusterFinder(pixel_dict, 256, 256).getListOfKlusters()[11]
# Make the image!
make_kluster_image(beta_cluster_1)
Finally, we have gamma radiation - photons that set off single pixels, or two, three of four when they hit the boundaries between pixels:
In [7]:
# An example of a tripixel gamma candidate.
gamma_cluster_1 = KlusterFinder(pixel_dict, 256, 256).getListOfKlusters()[25]
# Make the image!
make_kluster_image(gamma_cluster_1)
We can use the KlusterFinder Python object to find the clusters of pixels in the data, and the make_kluster_image function to make an image of the cluster that lets us look at it in more detail. Why not try looking at more clusters yourself?
In [8]:
# Exercise 2: Making your own cluster images
#
# Use code like the following to make your own cluster images from the data!
#
# First, let's make a Python list of the clusters:
klusters = KlusterFinder(pixel_dict, 256, 256).getListOfKlusters()
# We can now make an image by specifying the index in the list.
# Try changing the number 12 to another number to see the different
# clusters!
make_kluster_image(klusters[12])
It's fun to make the pictures of the clusters, but we can actually use code to do some important - but, dare we say, tedious - work for us. For example, how did you get on with the particle counting in Exercise 1? Our code makes tasks like this a lot easier...
In [9]:
print("Number of clusters (particle) in the data frame: %d" % (len(klusters)))
Did you get it right?
Now, it turns out KlusterFinder is a bit cleverer than just counting the particles. It also does a bunch of calculations on the clusters that tell us something about the clusters themselves. For example, it calculates:
We can make a DataFrame of these properties of the clusters and present it as a table:
In [10]:
# Import a function that makes a DataFrame of kluster information.
from cernatschool.helpers import get_dataframe_of_klusters
## The cluster DataFrame.
kdf = get_dataframe_of_klusters(klusters)
# Show the whole table, representing all 34 clusters.
kdf
Out[10]:
That's quite a lot of information to take in - but you can use code to help. Want to find only the particles above a certain size? No problem.
In [11]:
## Let's look for particles with a minimum size of 50 pixels.
min_size = 50
# Apply a filter to the table using the `loc` method and a conditional statement.
kdf.loc[kdf['size'] >= min_size]
Out[11]:
Which particles are on the edge of the frame?
In [12]:
kdf.loc[(kdf['isedgekluster']==True)]
Out[12]:
Which we can check with:
In [13]:
make_kluster_image(klusters[9])
We can also combine our search criteria by using boolean operators in the loc method like so:
In [14]:
## The minimum particle size.
min_size = 25
## The minimum density.
min_density = 0.5
# Find clusters that aren't on the edge, have a minimum size of 25 pixels and
# a minimum density of 0.8 using multiple conditions linked by AND & operators.
kdf.loc[(kdf['isedgekluster']==False ) & ((kdf['size'] >= min_size) & (kdf['density_uw'] >= min_density))]
Out[14]:
In [15]:
# Exercise 3: Particle spotting
#
# Use the make_kluster_image function to examine these particles more closely.
#
# What do you notice?
Hmm... so we seem to be able to sort the particles into groups by requiring that they meet certain criteria. We can make use of this to automatically create lists of clusters that represent these groups, and then count them:
In [16]:
## A list of big particles!
big_clusters = []
# Loop over the clusters.
for kluster in klusters:
# Apply the cut on cluster size.
if kluster.getNumberOfPixels() > 50:
# If it passes, add it to the list!
big_clusters.append(kluster)
print("%d big clusters found!" % (len(big_clusters)))
In [17]:
# Exercise 4: Checking our working
#
# Use the DataFrame of clusters to check that there are indeed four "big clusters"
# as defined here.
The other methods to get the cluster properties from a cluster are:
getNumberOfPixels() - the cluster size;getRadiusUW() - the cluster radius;getDensityUW() - the cluster density;getInnerPixelFraction() - the fraction of pixels that are completely surrounded by other pixels;isEdgeCluster() - returns True if the cluster is on the edge of the frame.So we can use multiple if statements and these values to sort the clusters in our for loop.
In [18]:
# Exercise 5: GRAND CHALLENGE - COUNT THE ALPHAS, BETAS, AND GAMMAS
#
# Can you come up with a for loop containing a series of if statements
# that will sort the clusters into four lists?
#
# We've made a start for you here - good luck!
## A list of single alpha particles.
alphas = []
## A list of beta particles.
betas = []
## A list of gammas.
gammas = []
## A list of "other" particles.
others = []
# Add to the loop here:
for kluster in klusters:
# Dummy if statement - you can do better!
if kluster.getNumberOfPixels() >= 1:
alphas.append(kluster)
print("Number of alphas = % 3d" % (len(alphas)))
print("Number of betas = % 3d" % (len(betas)))
print("Number of gammas = % 3d" % (len(gammas)))
print("Number of 'others' = % 3d" % (len(others)))
print("------------------------")
print("TOTAL = % 3d" % (len(alphas) + len(betas) + len(gammas) + len(others)))