Uncertainty on stereo parameters (h_max)

[Elisa-Visentin]

Please describe the use case that requires this feature. This could help in event reconstruction through Random Forests (and in optimizing reconstruction algorithms)

Describe the solution you'd like Maybe a new key in ReconstructedGeometryContainer

Additional context While comparing different Random Forest implementations, some of the results are 'anomalous' (it looks like the higher the number of telescopes detecting an event, the worse is the stereo reconstruction): knowing the uncertainty on this parameter could help in improving training and event reconstruction.

I hope this won't bother you too much... (I can't implement it by myself, sorry)

[maxnoe]

it looks like the higher the number of telescopes detecting an event, the worse is the stereo reconstruction

That is strange in some respects, generally the reconstruction should improve with more telescopes, but:

• More telescopes triggered usually means higher energy
• Higher energies have a larger likelihood of being truncated

In ctapipe itself, we can apply selections on the telescope events in several steps, so that we don't include truncated images in the stereo reconstructions (both for applying the HillasReconstructor and for the machine learning models).

While providing an uncertainty would be in any case helpful, I doubt that it directly addresses the problem you mention

[maxnoe]

To comment on the actual issue here: the ReconstructedGeometryContainer already has fields for uncertainties of all reconstructed values, see: https://github.com/cta-observatory/ctapipe/blob/7a131ac1cbe2fa92b6e8ab14813f3c74f104a151/ctapipe/containers.py#L747

However, this might not be filled by all reconstructors and indeed, the HillasReconstructor does not compute an uncertainty for the h_max prediction.

In general, this is a very hard problem, since it is not straight forward how to compute and propogate the uncertainties from different steps.

For the specific case of the HillasReconstructor and h_max for example:

• h_max is calculated as the least squares intersection point of the "line of sight" of each telescope of the center of gravity hillas parameters. How do we compute uncertainties of the cog estimation to the uncertainty of the least squares intersection in 3d? Is the cog even the correct thing? Shouldn't it be the brightest pixel instead of the cog?

[Elisa-Visentin]

As for your last question, the cog is already an average over pixel coordinates weighted by their charge (if I understood the hillas code), while by using the brightest pixel some images could lead to absurd results that could disrupt the geometric reconstruction (e.g., the brightest pixel in a corner of the ellipse, far away from the cog).

As for the uncertainty, a coarse solution could be to propagate 'brutally' errors on the cog (pixels position could be error-free, the charge is a number of photo-electrons, so maybe Poissonian?) as a 'maximum error' (e.g., compute h_max both for the evaluated cog and for the cog+cog_error (all the telescopes) and take the difference as the error; hopefully telescope pointing directions should have a negligible error, and, if we don't move them and we produce consistent simulations, the uncertainty on their position should not matter). Maybe somebody has a better idea to evaluate this uncertainty properly.

This is a residual effect of my lab classes ('every measured and evaluated parameter must have a unit and an error'), sorry...

[maxnoe]

the cog is already an average over pixel coordinates weighted by their charge

Yes, but showers have usually a skewness and in that case the mean is not equal to the mode.

[Elisa-Visentin]

In case of asymmetric distributions (e.g., a Landau or something like this) a weighted mean (or a truncated mean, or a combination of both) will be 'safer' as the mode approximation than using the brightest pixel: a single pixel, in particular in real data, can return a really high charge value (PMTs are 'nice', Silicon detectors can be more annoying) even if it is located in the 'queue' of the distribution, and it could not be removed by 'masking'. But, instead of a single pixel, the brightest pixel cluster (summing the charge on how many pixels?) could be used to 'level out' anomalous peaks