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
130
131
|
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, ndaq=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.
"""
ntotal = nevals * nreps * ndaq
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,
ndaq=ndaq,
time_only=self.time_only)
# Normalize probabilities and put a floor to keep the log finite
hit_prob = hitcount.astype(np.float32) / ntotal
hit_prob = np.maximum(hit_prob, 0.5 / ntotal)
# 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 = ( (ntotal * 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()
|
