I have two separate algorithms (call them "A1" and "A2") which reconstruct the $(x, y)$-position of an event in a particle detector. I can test both of these algorithms on simulated events from a very accurate Monte Carlo of the experiment. Note that A1 and A2 work on different observables in the detector and their errors are not correlated (though the observables themselves are correlated.) There is no systematic bias with either algorithm and so -- in aggregate -- the reconstruction error in $x$ and $y$ is approximately zero over all events.

Say A1 reconstructs a given MC event at some position $(x_1, y_1)$, and A2 reconstructs the same MC event at some different position $(x_2, y_2)$. I run A1 and A2 on every MC event, and end up with two distributions: one which gives the distance (scalar distance, not vector) from the true MC position of the event to $(x_1, y_1)$ for all events, and one which gives the distance from the true MC position of the event to $(x_2, y_2)$ for all events. Since this distance is necessarily non-negative, these distributions have some positive mean value, and some RMS.

This is all fine: these two distributions each have a mean (which characterizes the accuracy of the algorithm) and an RMS (which characterizes the precision of the algorithm), and the two distributions are more or less Gaussian. Next, I want to use these algorithms A1 and A2 on real data from the detector, and use the properties of these MC distributions to put a limit on the uncertainty of my reconstructed positions.

My question is this: knowing the RMSes and means of the A1 and A2 distance distributions, when I use A1 and A2 on real data, I should be able to find a "best fit" point and put some uncertainty on it. One one hand, I feel like I should just average the measurements and sum the errors in quadrature, but this feels incorrect for some reason (perhaps because the mean of A2 tends to be much higher than the mean of A1).

Is this the correct way to go about analyzing this data? Am I overthinking this?

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Manishearth Jun 3 '16 at 13:10

As explained in the comments, for a given MC event, the errors $\delta x$ and $\delta y$ for the two reconstruction methods A1 and A2 are uncorrelated. Therefore, no extra information is gained by using both reconstructions. We should just use the better one.

Now, what does "better" mean. As stated in the OP, for both reconstructions, the distributions of $\delta x$ and $\delta y$ are Gaussian with zero mean. However, the standard deviations are different. Assuming that for either reconstruction method $\sigma_{\delta x} = \sigma_{\delta y} \equiv \sigma$, the best strategy is to simply pick the method with the smaller $\sigma$. Unlike $\delta x$ and $\delta y$, the error "distance" $\delta r \equiv \sqrt{\delta x^2 + \delta y^2}$ has a distribution which is strictly positive and has nonzero mean.$^{[a]}$ As pointed out in the OP, A1 and A2 lead to two different distributions of $\delta r$, each with its own (positive) mean and variance. Note, however, that the mean and variance of the Rayleigh distribution are not independent (see article linked below), so really you just want to use the reconstruction method with the smaller mean $\delta r$, which is the same as picking the one with the smaller Gaussian $\sigma$ as described above.

Our statement above that we should just use the reconstruction method which has the smaller $\sigma$ critically depends on the fact that A1 and A2 use the same MC data. However, OP mentioned that in the real experiment A1 and A2 use different observables in the detector. In this case, we do benefit by combining the data from both methods, in essence because more data means more information. The question then is how to combine the results. This is actually a really interesting question.

For a given detector event, we get four numbers: $$x_1, \, y_1, \, x_2, \, y_2$$ where $x_i$ is the $x$ coordinate reconstructed by method Ai, and similarly for $y$. Now, based only on the results from A1, our knowledge of the real $x$ position is

$$P_{X_1}(x|x_1) = \frac{1}{\sqrt{2\pi \sigma_1^2}} \exp \left[-\frac{(x-x_1)^2}{2 \sigma_1^2}\right] \, .$$ The expressions for $P_{Y_1}$, $P_{X_2}$, and $P_{Y_2}$ are the obvious similar things. Now the cool part: the distribution of the $x$ position of the event, given both measurements A1 and A2 is $$P_{X_{1+2}}(x|x_1,x_2) = \frac{1}{\mathcal{N}} P_{X_1}(x|x_1) P_{X_2}(x|x_2)$$ where $\mathcal{N}$ is a normalization factor.

Note that if either of A1 or A2 is much better than the other, you could just ignore the worse one without losing much precision.

$[a]$: For what it's worth, the distribution of $\delta r$ is called the Rayleigh distribution.

  • $\begingroup$ Forgive my ignorance, but, given that $\sigma_1$, $\sigma_2$, $x_1$, $y_1$, $x_2$, and $y_2$ are all constants, I want to just substitute them into the above formula and find the values of $x$ and $y$ which give the maximum probability, correct? My next question is: how do I combine the sigmas to get an error bound for the resulting $x$ and $y$? Just take the root of the sum of the squares? $\endgroup$ – Andrew Watson Jun 1 '16 at 3:48
  • $\begingroup$ @AndrewWatson First, yes, the formula I gave for $P_{X_{1+2}}$ can be written out explicitly with all the $x_1$, $x_2$, etc. values. I didn't do it because I trust that you can ;) Second, once you write out that formula, you can certainly find the maximum probability point by the usual means (i.e. take derivative and set to zero). Finally, you ask about combining sigmas. Be careful! When $x_1 \neq x_2$ the distribution for $x$ is no longer a Gaussian so you can't necessarily just add the variances. $\endgroup$ – DanielSank Jun 1 '16 at 3:57
  • $\begingroup$ @AndrewWatson Remember that because you measure each event twice you have a distribution which depends on two spacial points. If you want to know the variance for a given measured pair of $x_1$ and $x_2$, you can just numerically do the appropriate integral. I'm not sure if there's a closed form representation. There might be! $\endgroup$ – DanielSank Jun 1 '16 at 4:00
  • $\begingroup$ Great! Thanks! And yes I realized that after I posted my comment above. I don't really want much more than a single number for the error (I don't want a 2D distribution) for each event, so I can just ignore that for now. $\endgroup$ – Andrew Watson Jun 1 '16 at 4:02
  • $\begingroup$ @AndrewWatson Happy to help. This was a very fun question! I would ask if you might go back and tidy up the original post. For example, on this site we don't put "EDIT" tags in questions; we just edit the question :) If you could please do this and integrate the useful information from the comments into the main post, that would be much appreciated. We can then ask a mod to delete the (rather long) comment chain. $\endgroup$ – DanielSank Jun 1 '16 at 4:04

To me, using the distance to measure 2 dimensional data is dubious, especially if the events are near to the reconstructed position (supposed 0 distance). Say if the events form a circle (the circle is perfect gaussian). The distribution will not be anymore gaussian since the blue area (in the drawing) will be lesser than the yellow:

enter image description here

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.