diff options
Diffstat (limited to 'likelihood.py')
| -rw-r--r-- | likelihood.py | 10 |
1 files changed, 6 insertions, 4 deletions
diff --git a/likelihood.py b/likelihood.py index 51d7ef9..08b090f 100644 --- a/likelihood.py +++ b/likelihood.py @@ -43,13 +43,14 @@ class Likelihood(object): self.event = event @profile_if_possible - def eval(self, vertex_generator, nevals, nreps=1): + 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) @@ -59,11 +60,12 @@ class Likelihood(object): 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) / (nreps * nevals) - hit_prob = np.maximum(hit_prob, 0.5 / (nreps*nevals)) + 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) @@ -79,7 +81,7 @@ class Likelihood(object): # 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 + 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 |