crappy pixel noise simulation

crap_mc.py

#!/usr/bin/env python

import numpy as np

'''
 Try to simulate the nearAvg5 distribution with pixel
 noise modelled as a gamma distribution.
'''
def simulate_navg(threshold=40., nsamples=100000, mean=25., shape=2., gridsize=5):
  if not gridsize%2:
    raise ValueError('gridsize should be odd')

  # get a NxN array of zeros with a 1 in the middle
  mid = int(gridsize)/2
  center = np.zeros((gridsize, gridsize))
  center[mid,mid] = 1

  # iterator to generate NxN random samples
  trials = (np.int_(np.random.gamma(shape, float(mean)/shape, (gridsize, gridsize)))
    for i in xrange(nsamples))

  # trigger on events where the middle is greater than threshold,
  triggered = (t for t in trials if t[mid, mid] > threshold)

  # subtract the center for triggered events
  subtracted = (t - t[mid,mid]*center for t in triggered)

  # return the means
  return [t.mean() for t in subtracted]


'''
 run as a program. sample the distribution and plot it as a TH1F.
'''
if __name__ == '__main__':
  import ROOT as r

  h = r.TH1F('n5', 'n5', 100, 0, 50)
  map(h.Fill, simulate_navg())

  h.Draw()
  r.gROOT.FindObject('c1').SetLogy()

  print "Press enter to continue..."
  raw_input()