This is an instructor notebook (written by BWO, 1/26/2016) that mocks up a radioactive decay of some substance including a background, and writes out two files: one containing count rates as a function of time (count_rates.txt) and the other containing a much smaller estimate of the amount of radioactive substance as a function of time (sample_amount.txt). The count rate includes the expected radioactive decay plus an arbitrary positive noise background, and the sample amount has an arbitrary measurement error built in (no offset, just noise).

The point of this notebook is that the actual assignment has students calculating radioactive decay by hand, and then reading some data, looking at it, and trying to figure out what is going on in the data (i.e., why there is noise, and if the amount of stuff you'd expect based on the measured count rate actually lines up with the measured amount of stuff).



In [ ]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import numpy.random as rdm
import math

initial_amount = 1.0e+6  # initial number of atoms
half_life = 1712.0       # half-life in seconds (or whatever units)

stop_time = 20000.0      # stop time for data output; make this several half-lives

# create a value that is sensible for the noise rate (i.e., not dominating but not negligible)
noise_norm = initial_amount/half_life/10.0

# interval between measurements (same units as half_life)
dt = 5.0

# subsampling of outputs
subsample = 180

time = np.arange(0.0,stop_time,dt)



In [ ]:

    
# expected amount of stuff as f(t)(no noise)  (float)
amount_of_stuff = initial_amount * 0.5**(time/half_life)   

# expected count rate per dt, as f(t) (no noise)  (float)
count_rate = initial_amount / half_life * math.log(2.0) * 0.5**(time/half_life) * dt



In [ ]:

    
# calculate the background count rate per dt.  (float; scale of 4.0 is arbitrary)
background_counts = rdm.normal(noise_norm,scale=noise_norm/4.0,size=time.size)*dt 

# make sure background counts are always positive
background_counts[background_counts<0] = np.abs(background_counts[background_counts<0.0])

# add noise to the amount of stuff; basically just mocking up poisson sampling error
# (This really doesn't show up, it's just amusing to do)
noisy_amount = amount_of_stuff + rdm.normal(0.0,scale=np.sqrt(amount_of_stuff),size=time.size) 

# add noise to the count rate to make it look real(ish), and convert to ints
noisy_counts = (count_rate + background_counts).astype(dtype=int)

# convert non-noisy counts to integers as well (since counts have to be ints in our output)
count_rate = count_rate.astype(dtype=int)

# some sanity checks
plt.plot(time,background_counts,'b.')#,x,noise,'r.')
plt.plot(time,noisy_counts,'r-')
plt.plot(time,count_rate,'g-',linewidth=3)

print(noisy_counts.min(), noisy_counts.max(), noisy_counts.mean())



In [ ]:

    
# What does the background look like?  Just checking...  (histogram of counts per interval)

n,bins,patches=plt.hist(background_counts,100)



In [ ]:

    
# convert amount of stuff to ints

noisy_amount = noisy_amount.astype(dtype=int)
amount_of_stuff = amount_of_stuff.astype(dtype=int)

plt.plot(time,noisy_amount,'r-',time,amount_of_stuff,'b-')
print("noisy amount:", noisy_amount.min(),noisy_amount.max(),noisy_amount.mean())



In [ ]:

    
# subsample time and noisy_amount, because we don't want to make the students' lives too easy and
# give them everything with the same number of elements on the x-axis!

# note: subsample defined up at the top.
time_cut = time[::subsample].astype(dtype=int)
noisy_amt_cut = noisy_amount[::subsample]  # already an int



In [ ]:

    
# do a quick test to make sure that when you naively remove counts from 
# the total amount of stuff it ends up being negative (est_amount)
# also do a simple integration of the noiseless solution (calc_amount)
# to double-check we're doing things right.

mytime = np.arange(0,stop_time,dt)
est_amount = np.zeros_like(mytime)

calc_amount = np.zeros_like(mytime)

print(mytime.size, est_amount.size, noisy_counts.size)

est_amount[0] = noisy_amount[0]

calc_amount[0] = initial_amount 


for i in range(1,mytime.size-2):
    est_amount[i] = est_amount[i-1] - noisy_counts[i-1]  # noisy_counts already has a dt in it
    calc_amount[i] = calc_amount[i-1] - initial_amount / half_life * math.log(2.0) * 0.5**(mytime[i]/half_life) * dt
    
# wow, that sure looks bad, doesn't it?
plt.plot(mytime,est_amount,'b.',time_cut,noisy_amt_cut,'r-',mytime,calc_amount,'k--',linewidth=3)



In [ ]:

    
# write out count rates vs. time
combined_array = np.array([time,noisy_counts])
np.savetxt("count_rates.txt", combined_array, fmt='%d', 
           delimiter=' ', newline='\n', header='', footer='', comments='# ')

# write out amount of sample vs. time
combined_array_2 = np.array([time_cut,noisy_amt_cut])
np.savetxt("sample_amount.txt", combined_array_2, fmt='%d', 
           delimiter=' ', newline='\n', header='', footer='', comments='# ')



In [ ]:

    
# just checking to make sure this works.

x = np.loadtxt("sample_amount.txt", dtype=int)[0]
y = np.loadtxt("sample_amount.txt", dtype=int)[1]

plt.plot(x,y,'ro')



In [ ]: