### 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


params = Parameters()