diff options
Diffstat (limited to 'likelihood.py')
| -rw-r--r-- | likelihood.py | 54 |
1 files changed, 36 insertions, 18 deletions
diff --git a/likelihood.py b/likelihood.py index c797afe..51d7ef9 100644 --- a/likelihood.py +++ b/likelihood.py @@ -2,7 +2,7 @@ import numpy as np from histogram import HistogramDD from uncertainties import ufloat, umath, unumpy from math import sqrt -from itertools import izip, compress +from itertools import islice, izip, compress, repeat from chroma.tools import profile_if_possible class Likelihood(object): @@ -47,16 +47,19 @@ class Likelihood(object): """ Return the negative log likelihood that the detector event set in the constructor or by set_event() was the result of a particle generated - by `vertex_generator`. If `nreps` set to > 1, each set of photon vertices - will be propagated `nreps` times. + by `vertex_generator`. If `nreps` set to > 1, each set of photon + vertices will be propagated `nreps` times. """ - hitcount, pdf_prob, pdf_prob_uncert = self.sim.eval_pdf(self.event.channels, - nevals, vertex_generator, - 2.0e-9, self.trange, - 1, self.qrange, - nreps=nreps, - time_only=self.time_only) + vertex_generator = islice(vertex_generator, nevals) + + hitcount, pdf_prob, pdf_prob_uncert = \ + self.sim.eval_pdf(self.event.channels, + vertex_generator, + 2.0e-9, self.trange, + 1, self.qrange, + nreps=nreps, + time_only=self.time_only) # Normalize probabilities and put a floor to keep the log finite hit_prob = hitcount.astype(np.float32) / (nreps * nevals) @@ -79,7 +82,8 @@ class Likelihood(object): hit_channel_prob_uncert = ( (nevals * nreps * hit_prob * (1.0 - hit_prob)) / hit_prob**2 ).sum()**0.5 log_likelihood = ufloat((hit_channel_prob, 0.0)) - # Then include the probability densities of the observed charges and times + # Then include the probability densities of the observed + # charges and times. hit_pdf_ufloat = unumpy.uarray((pdf_prob[self.event.channels.hit], pdf_prob_uncert[self.event.channels.hit])) log_likelihood += unumpy.log(hit_pdf_ufloat).sum() @@ -91,21 +95,35 @@ if __name__ == '__main__': from chroma.sim import Simulation from chroma.optics import water from chroma.generator import constant_particle_gun - from chroma.tools import enable_debug_on_crash + from chroma import tools import time - enable_debug_on_crash() + tools.enable_debug_on_crash() detector = detectors.find('lbne') - sim = Simulation(detector, water) + sim = Simulation(detector, seed=0) - event = sim.simulate(1, constant_particle_gun('e-',(0,0,0),(1,0,0),100.0)).next() + event = sim.simulate(islice(constant_particle_gun('e-',(0,0,0),(1,0,0),100.0), 1)).next() print 'nhit = %i' % np.count_nonzero(event.channels.hit) likelihood = Likelihood(sim, event) - for x in np.linspace(-10.0, 10.0, 2*10*5+1): - start = time.time() - l = likelihood.eval(constant_particle_gun('e-',(0,0,0),(1,0,0),100.0),100, nreps=200) - print 'x = %5.1f, %s (%1.1f sec)' % (x, l, time.time() - start) + x = np.linspace(-10.0, 10.0, 100) + l = [] + + for pos in izip(x, repeat(0), repeat(0)): + t0 = time.time() + ev_vertex_iter = constant_particle_gun('e-',pos,(1,0,0),100.0) + l.append(likelihood.eval(ev_vertex_iter, 1000)) + elapsed = time.time() - t0 + + print '(%.1f, %.1f, %.1f), %s (%1.1f sec)' % \ + (pos[0], pos[1], pos[2], tools.ufloat_to_str(l[-1]), elapsed) + + import matplotlib.pyplot as plt + + plt.errorbar(x, [v.nominal_value for v in l], [v.std_dev() for v in l]) + plt.xlabel('X Position (m)') + plt.ylabel('Negative Log Likelihood') + plt.show() |