1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
import numpy as np
from histogram import HistogramDD
from uncertainties import ufloat, umath, unumpy
from math import sqrt
from itertools import islice, izip, compress, repeat
from chroma.tools import profile_if_possible
class Likelihood(object):
"Class to evaluate likelihoods for detector events."
def __init__(self, sim, event=None, tbins=100, trange=(-0.5e-9, 999.5e-9),
qbins=10, qrange=(-0.5, 9.5), time_only=True):
"""
Args:
- sim: chroma.sim.Simulation
The simulation object used to simulate events and build pdfs.
- event: chroma.event.Event, *optional*
The detector event being reconstructed. If None, you must call
set_event() before eval().
- tbins: int, *optional*
Number of time bins in PDF
- trange: tuple of 2 floats, *optional*
Min and max time to include in PDF
- qbins: int, *optional*
Number of charge bins in PDF
- qrange: tuple of 2 floats, *optional*
Min and max charge to include in PDF
- time_only: bool, *optional*
Only use time observable (instead of time and charge) in computing
likelihood. Defaults to True.
"""
self.sim = sim
self.tbins = tbins
self.trange = trange
self.qbins = qbins
self.qrange = qrange
self.time_only = time_only
if event is not None:
self.set_event(event)
def set_event(self, event):
"Set the detector event being reconstructed."
self.event = event
@profile_if_possible
def eval(self, vertex_generator, nevals, nreps=1):
"""
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.
"""
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)
hit_prob = np.maximum(hit_prob, 0.5 / (nreps*nevals))
# Set all zero or nan probabilities to limiting PDF value
bad_value = (pdf_prob <= 0.0) | np.isnan(pdf_prob)
if self.time_only:
pdf_floor = 1.0 / (self.trange[1] - self.trange[0])
else:
pdf_floor = 1.0 / (self.trange[1] - self.trange[0]) / (self.qrange[1] - self.qrange[0])
pdf_prob[bad_value] = pdf_floor
pdf_prob_uncert[bad_value] = pdf_floor
print 'channels with no data:', (bad_value & self.event.channels.hit).astype(int).sum()
# NLL calculation: note that negation is at the end
# Start with the probabilties of hitting (or not) the channels
hit_channel_prob = np.log(hit_prob[self.event.channels.hit]).sum() + np.log(1.0-hit_prob[~self.event.channels.hit]).sum()
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.
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()
return -log_likelihood
if __name__ == '__main__':
from chroma import detectors
from chroma.sim import Simulation
from chroma.optics import water
from chroma.generator import constant_particle_gun
from chroma import tools
import time
tools.enable_debug_on_crash()
detector = detectors.find('lbne')
sim = Simulation(detector, seed=0)
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)
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()
|
