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
|
import numpy as np
from histogram import HistogramDD
from uncertainties import ufloat, umath, unumpy
from math import sqrt
from itertools import izip, compress
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.
"""
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)
# 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.tools import enable_debug_on_crash
import time
enable_debug_on_crash()
detector = detectors.find('lbne')
sim = Simulation(detector, water)
event = sim.simulate(1, constant_particle_gun('e-',(0,0,0),(1,0,0),100.0)).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)
|
