fix bug in Likelihood class; pdfs for each pmt were the same histogram

1 files changed, 9 insertions, 7 deletions

diff --git a/likelihood.py b/likelihood.py
index f0323fe..c6b6d99 100644
--- a/likelihood.py
+++ b/likelihood.py
@@ -19,7 +19,7 @@ class Likelihood(object):
         if event is not None:
             self.set_event(event)
 
-        self.pdfs = dict(zip(sim.detector.pmtids,[HistogramDD((100,10), range=[(0.0,1000.0),(-0.5,9.5)])]*len(sim.detector.pmtids)))
+        self.pdfs = dict(zip(sim.detector.pmtids,[HistogramDD((100,10), range=[(0.0,500.0),(-0.5,9.5)]) for i in sim.detector.pmtids]))
 
     def set_event(self, event):
         "Set the detector event being reconstructed."
@@ -35,8 +35,8 @@ class Likelihood(object):
             pdf.reset()
 
         for ev in self.sim.simulate(nevals, vertex_generator):
-            for i, t, q in compress(izip(range(ev.hits.hit.size),ev.hits.t,ev.hits.q),ev.hits.hit):
-                self.pdfs[i].fill((ev.hits.t[i],ev.hits.q[i]))
+            for i, t, q in compress(izip(range(ev.channels.hit.size),ev.channels.t,ev.channels.q),ev.channels.hit):
+                self.pdfs[i].fill((ev.channels.t[i],ev.channels.q[i]))
 
         for pdf in self.pdfs.values():
             if pdf.nentries > 0:
@@ -51,7 +51,7 @@ class Likelihood(object):
                 if self.pdfs[i].nentries > 0:
                     probability = ufloat([0.5/self.pdfs[i].nentries]*2)
                 else:
-                    probability = ufloat([1.0/len(self.pdfs[i].bincenters)]*2)
+                    probability = ufloat([1.0/self.pdfs[i].hist.size]*2)
 
             log_likelihood += umath.log(probability)
 
@@ -63,12 +63,14 @@ if __name__ == '__main__':
     from chroma.optics import water
     from chroma.generator import constant_particle_gun
 
-    detector = detectors.find('minilbne')
+    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.hits.hit)
+
     likelihood = Likelihood(sim, event)
 
-    for x in np.linspace(-1.0, 1.0, 20):
-        print 'x = %f, %s' % (x, likelihood.eval(constant_particle_gun('e-',(x,0,0),(1,0,0),100.0),100))
+    for x in np.linspace(-10.0, 10.0, 10*5+1):
+        print 'x = %5.1f, %s' % (x, likelihood.eval(constant_particle_gun('e-',(x,0,0),(1,0,0),100.0),100))