Geiger-Muller Detector and Counting Statistics in Python with Chi Square minimization


Posted by 4 years ago
426.1k points

Using a Geiger-Muller counter enables us to count the number of radiation events from some source. A few months ago I tried this experiment with Cesium-137 and interpolated the data to see if it fits a Gaussian/Poisson curve. After performing the experiment the data were as follows,

264 1 294 15 324 13
265 1 295 18 325 16
266 4 296 22 326 13
267 2 297 23 327 6
268 2 298 24 328 11
269 2 299 16 329 6
270 2 300 25 330 11
271 2 301 15 331 6
272 5 302 27 332 14
273 3 303 28 333 5
274 2 304 24 334 5
275 6 305 22 335 3
276 7 306 22 336 9
277 3 307 21 337 4
278 7 308 21 338 4
279 5 309 23 339 6
280 9 310 36 340 4
281 4 311 21 341 2
282 8 312 12 342 4
283 11 313 19 343 3
284 9 314 15 344 1
285 15 315 23 345 3
286 19 316 19 346 3
287 11 317 13 347 2
288 12 318 17 348 1
289 10 319 13 349 1
290 13 320 19 350 1
291 18 321 19 351 1
292 20 322 18 352 2
293 14 323 17

Which was the number of radiation events per time. The data were interpolated using Python. Here are some important code snippets,

from lmfit import minimize, Parameters, fit_report  
from numpy import * 
import seaborn as sns  
from matplotlib import rc  
import matplotlib as mpl  
import scipy as sc  
from scipy.interpolate import interp1d  

Using Scipy may have worked on it's own, without the help of lmfit, but for me this was made easier using lmfit. I defined my functions as follows,

def gaussian(x,amp,mu,sigma):
    return (amp/(sqrt(2*pi)*sigma)) * exp(-(x-mu)**2 /(2*sigma**2))

def poisson(x,amp,mu):
    return amp*exp(-mu)*(mu**x)/sc.special.gamma(x+1)

def chisqg(params,x,y,err):
    amp = params['amp'].value
    mu = params['mu'].value
    sigma = params['sigma'].value
    model = (amp/(sqrt(2*pi)*sigma)) * exp(-(x-mu)**2 /(2*sigma**2))
    return (y-model)/err

def chisq(params,x,y,err):
    amp = params['amp'].value
    mu = params['mu'].value
    model = amp*exp(-mu)*(mu**x)/sc.special.gamma(x+1)
    return (y-model)/err

and additionally,

params = Parameters()
params.add('amp',value=1000)
params.add('mu',value=2)
out = minimize(chisq,params,args=(x,y,err,))

f = interp1d(x,poisson(x,amp,mu),kind='cubic')
f2 = interp1d(x,gaussian(x,ampg,mug,sigma),kind='cubic')

these snippets along with the help of libraries already mentioned gives the following result,

pgl.png