diff options
| author | Stan Seibert <stan@mtrr.org> | 2011-12-21 11:36:25 -0500 |
|---|---|---|
| committer | tlatorre <tlatorre@uchicago.edu> | 2021-05-09 08:42:38 -0700 |
| commit | 127dda0bae74829eed4b274797db26e6741ffdb5 (patch) | |
| tree | 38396dafa0de363eb7887b74fcb4c348e58cdd2b | |
| parent | e150de0d229401b045bf17cd5a1b658f69b2230c (diff) | |
| download | chroma-127dda0bae74829eed4b274797db26e6741ffdb5.tar.gz chroma-127dda0bae74829eed4b274797db26e6741ffdb5.tar.bz2 chroma-127dda0bae74829eed4b274797db26e6741ffdb5.zip | |
Apply a floor to the probability of each channel being in the observed
state (hit or not), rather than a floor on the hit probability.
A channel that is impossible to hit should have zero probability in
the likelihood and not be hit in the actual data. Before this change,
that channel would be forced to have a hit probability of 0.5 /
ntotal, which is wrong. We only need to ensure the probability of the
observed state of the channel is not zero so that the log() function
is defined.
| -rw-r--r-- | chroma/likelihood.py | 13 |
1 files changed, 9 insertions, 4 deletions
diff --git a/chroma/likelihood.py b/chroma/likelihood.py index 506131d..5b746dd 100644 --- a/chroma/likelihood.py +++ b/chroma/likelihood.py @@ -64,9 +64,8 @@ class Likelihood(object): time_only=self.time_only, min_bin_content=320) - # Normalize probabilities and put a floor to keep the log finite + # Normalize probabilities 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) @@ -88,14 +87,20 @@ class Likelihood(object): by `vertex_generator`. If `nreps` set to > 1, each set of photon vertices will be propagated `nreps` times. """ + ntotal = nevals * nreps * ndaq hit_prob, pdf_prob, pdf_prob_uncert = \ self.eval_channel_vbin(vertex_generator, nevals, nreps, ndaq) # 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() + + # Flip probability for channels in event that were not hit + hit_prob[~self.event.channels.hit] = 1.0 - hit_prob[~self.event.channels.hit] + # Apply a floor so that we don't take the log of zero + hit_prob = np.maximum(hit_prob, 0.5 / ntotal) + + hit_channel_prob = np.log(hit_prob).sum() log_likelihood = ufloat((hit_channel_prob, 0.0)) # Then include the probability densities of the observed |