MPE Pandel likelihood

Introduction

In order to try to understand the trends in the performance of the MPE, the MPE likelihood was investigated on a lower level. Plots of the MPE likelihood versus the time residual of the first hit were made, for various combinations of the perpendicular distance and the total charge in the DOM.

UPDATE Dima found the solution (better and still fast approximation of the integral of the convoluted Pandel function) in the appendix of a paper by Mathieu and George. This is now implemented both in rime and in ipdf/lilliput (select "GaussConvolutedFastApproximation" for the "PEProb" option of the I3GulliverIPDFPandel service). The approximation is quite good, as can be seen in the plots below.

The plots below are made for different versions of Pandel or different ways to compute/approximate the integral of the PDF:

Numerical Gauss-convoluted: the supposedly most correct computation is obtained by using the Gauss convoluted Pandel (new implementation), where the integral of the PDF is computed by numerically integrating the PDF.
Plain Pandel: this does not take into account jitter. The integral has an easy analytical expression.
Gauss semi-convoluted Pandel: this takes jitter into account in the PDF by a convolution of the Pandel PDF with a Gaussian (in the plots on this page we use jitter=sigma=4ns). The PDF itself is exactly computable, mathematically, but numerically it is approximated in different ways for larger distances and time residuals. The integral does not have an analytical expression (as far as we know). A crude and fast approximation is to use the easy analytical expression of the plain Pandel integral instead. This strategy is called "semi-convoluted" because the PDF takes jitter into account while the integral does not.
Gauss convoluted Pandel with Mathieu/George/Dima's fast approximation: Like above, except that the integral is computed following a prescription by Mathieu and George in the appendix of astro-ph/0506136. they rewrite the integral of the Gauss-convoluted Pandel as a sum of incomplete Gamma functions plus a residual integral; and Dima found that the residual integral can be computed numerically with satisfactory precision using only 4 evaluations of the integrand.
Box-convoluted Pandel: takes jitter (here we set this to 4ns) into account by convoluting the Pandel function with a rectangular jitter-PDF of width 4ns, that is, the jitter is assumed to be flatly distributed in the interval -2ns to +2ns (zero otherwise). Both the PDF and the integral have simple analytical expressions that are relatively cheap to compute.

Discussion

Definition of proper MPE:
MPE_n(t_res) = n pdf(t_res) I(t_res)^n-1
I(t_res) = _{t_res}∫^∞pdf(t)dt
(Exercise for the student: prove that if pdf(t_res) is normalized to 1, then MPE_n(t_res) is also normalized to 1 for all n≥1.)

For MPE with "semi-convoluted" Pandel:
MPE_n(t_res) = n cpandel(t_res) I(t_res)^n-1
I(t_res) = _{t_res}∫^∞plainpandel(t)dt

Normalization

We are most interested in the MPE behavior with semi-convoluted Pandel. From the plots it is clear that the approximation (of using the plain Pandel integral) is good enough at large distances, but that at short distances (less than 20m) the normalization of the MPE likelihood is obviously off. It was kind of expected that the normalization would be off (since we use a different pdf in I(t_res)), but it was underestimated how much it would be off. Note that a constant normalization error would not be a serious problem, since this just shows up as an additive constant in the log. The real problem is that the normalization error depends on the distance from the DOM.

This distance-dependent normalization error may well explain why the Paraboloid sigma with this version of MPE gives a too optimistic resolution estimate: the minimizer finds the track which passes as close by as many DOMs as possible with a slightly negative t_res, and for the tracks on the paraboloid grid the DOMs are probably on average further away, so that the MPE normalization errors are smaller, which makes the the MPE minimum seem steeper.

(Apparent) discontinuity at t_res=0

The integral I(t_res) of the plain Pandel function is equal to 1.0 for negative t_res, and decreases at positive t_res values. For DOMs with very large NPE, the decrease from 1.0 to <1.0 results in a sudden drop in MPE. Moreover, at small distances, the Pandel function is strongly peaked, so the integral I(_res) also decreases very quickly. These two effects makes MPE at t_res=0 so steep that it looks like it has a discontinuity. Even for NPE=2 the sudden decrease looks like a discontinuity.

Discontinuities are bad for the reconstruction, the likelihood landscape should be as smooth as possible. For large charges, the exact MPE would also force nearby tracks to have an even more negative time residual for the first hit, but for modest NPE values the (apparent) discontinuity will hamper the reconstruction process.

This might explain why the MPE reconstruction performs well at high energies, but has slightly worse performance than SPE at low energies.

Box-convoluted Pandel

The above described problems are not very critical, but they degrade the MPE performance and should be addressed. They are also in a way embarassing, because convoluted Pandel was designed precisely in order to eliminate problems at t_res=0; and by using "semi"-convoluted Pandel we have reintroduced them.

We have tried out semi-convoluted Pandel as an attempt to get a fast MPE reconstruction: computing the integral of the Gauss convoluted Pandel is CPU intensive. An alternative approach is to use the box-convoluted Pandel, as defined above. One could think of this PDF as taking for each t_res the average of the plain Pandel PDF over an interval [t_res-jitter/2,t_res+jitter/2].

The curves in the plots (right column) look a little odd at first, but when you think about it, it's clear that they should look like this. At short distances the distributions still have a steep edges (at t_res=jitter/2, because that's where t_res=jitter/2=2ns falls out of the "averaging window"), but not as severe as semi-convoluted Pandel, because MPE with box-convoluted Pandel is actually correctly normalized. The drawback is that like plain Pandel, there is a start of the distribution at t_res=-jitter/2=-2ns, though it is rounded, not sharp, so it's still better than plain Pandel.

Outlook

The problem seems to be solved using Mathieu/George/Dima's approximation for the integral of the convoluted Pandel function. Yay!

There are roughly three strategies we can pursue:

Leave it as it is, use MPE with semi-convoluted Pandel for mass processing and MPE with numerical Pandel for high level processing.
Test box-convoluted Pandel in reconstruction. This was tried before, but while making these plots I found that there was a serious bug in its implementation, so those tests need to be redone.
The discrepancies between the precise numerical Gauss convoluted Pandel and semi-convoluted Pandel only occur at small time residuals and small distances. So we could try to make lookup tables for that region of phase space, and use semi-convoluted Pandel for the rest of the cases.
The use of look-up tables was already pioneered in ipdf (by Robert Franke from DESY), but with a different logic: Robert's tables stretch all the way to large time residuals and large distances. Not all table entries yield precise enough values, especially at short distances and times, and for these cases the integral was computed numerically.
Instead we could try a scheme were we make fine grained tables which work for distances down to e.g. 1m and up to e.g. 40m. For tracks with a shorter distance than 1m we evaluate MPE with a distance of 1m. For distances larger than 40m we use semi-convoluted Pandel.

The consensus in Madison is that the last option seems the most reasonable. I am pursuing this right now. Unfortunately I found that the existing code for Pandel table generation and lookup code in ipdf needed quite some work, but I'm making progress. Once I'm done I'll post some new plots.

Plots

The MPE likelihood obtained with numerically Gauss convoluted Pandel is plotted in every figure as the reference, in shades of green (and dash-dotted). The other MPE curves are plotted in shades of blue.
The (SPE) Pandel PDF distributions are also given: Gauss-convoluted in red, plain (unconvoluted) in magenta, box-convoluted in orange.
Note that the y-axis scale is the same for all plots, but the x-axis varies with distance.
Click on the plots to get an (arbitrarily resizeable) encapsulated postscript version.
These plots are not friendly for color-blind people, my apologies.

Distance	plain, unconvoluted Pandel	NEW Gauss-convoluted Pandel with "FastApproximation" integral	Gauss "semi-convoluted" Pandel	"box"-convoluted Pandel
5m
10m
20m
40m
80m