HillasReconstructor and PSF

[jjlk]

Hi everyone, @thomasgas

The geometrical reconstruction with the method implemented in the HillasReconstructor class gives interesting results when one computes the point spread function (PSF) for event with multiplicity 2. You can see below:

• Left, PSF with a 68 % containment radius as a function of true energy for gamma-rays seen by the LST subarray for different event multiplicities
• Right, corresponding number of events

I'm using ctapipe 0.6.0 and for both plots, the events come from ON-axis gamma-ray simulation at the Southern site, for 20 deg zenith and North azimulth (production called, Prod3_Paranal_Baseline_NSB1x_gamma_North_20deg_DL0). I applied a cleaning base on tail-cuts (3/6) and in order to reconstruct the shower I selected the images with the following cuts:

• intensity >= 50 p.e.
• #pixels >=3
• 0.1 <= ellipticity <= 0.6
• nominal distance (distance between the center of the camera and the barycentre of the charges) <= 80 % of the camera radius

I didn't find an explaination for the 200 GeV 'bump' for events with multiplicity 2 by looking at the image or the shower properties. Note that the effect seems also present for multiplicty 3 but with a lower intensity.

@thomasgas can reproduce this behaviour independently.

Is there an expert that would be willing to have a look at this problem? I'll be happy to help by providing high-level diagnostics but I guess that a closer look at the reconstruction implementation will be necessary.

Thanks!

[thomasgas]

Hi Julien! Thanks for reporting this here!
My first guess at this is that when the impact point on the ground is more or less along the line connecting the two triggered telescopes, then the reconstruction is worst: maybe we see this just at lower energies due to the small images that we have on the cameras, but I'm just doing some speculation.
Experts are welcome!

[vuillaut]

maybe we see this just at lower energies due to the small images that we have on the cameras

Well, it's not at lowest energies.
There are more events below 100GeV and the reconstruction is better.

[thomasgas]

There are more events below 100GeV and the reconstruction is better.

Ah yes you are right

[thomasgas]

hi @jjlk, do you still have this issue with 2 LSTs?
Can @kpfrang weighting method using a LUT for the distance of closest approach be useful?

[kpfrang]

Hi all,
I'm not sure in the moment, when I'll have time to do it. But if you expect that the LUTs might help, I can run through some files to see how it looks.

[moralejo]

Hi,

I can reproduce the bump issue with the MARS analysis of Prod3b. I used 4 LSTs in the northern site, pointing south (most favourable azimuth B-field-wise, to compare with north-wise pointing from Paranal). The same image-wise cuts are applied. Caveat: I think the intensity cut (and the cleaning, too) probably have a different meaning in MARS and ctapipe. I believe ctapipe applies some correction to the conversion factors which depends on the number of integrated samples, to make the outcome closer to the actual number of detected p.e. in a pixel. MARS, uses instead the conversion factors from simtelarray as provided.

The color codes are the same as in the plots above: red, blue and green are respectively Ntels==2, ==3 and ==4. The different angular resolution values (worse for the MARS plots) may be related to the different reconstruction methods, and partly also to larger B-field in the northern site.

This is with no g/h separation cuts:

This is with hadronness<0.3, a soft cut which has gamma efficiency > 0.7 above 30 GeV

The plot below shows that with no g/h separation cut, and just requiring a minimum angle between the images in 2-tel events, does not solve the issue. I guess because the problem is not mainly images with their axes genuinely too parallel (due to large impact), but rather events in which at least one axis is poorly reconstructed (be it due to spurious islands, roundish image from low-impact or just too few usable pixels):

@HealthyPear @kosack can you confirm that the bump is mostly gone as soon as you apply some gamma/hadron separation?

[vuillaut]

Hi @moralejo
Good to know that this bump is reproducible with another analysis chain!

but rather events in which at least one axis is poorly reconstructed (be it due to spurious islands, roundish image from low-impact or just too few usable pixels)

Do you have an idea why this would happen at a particular energy though?

[moralejo]

Hi @moralejo
Good to know that this bump is reproducible with another analysis chain!

but rather events in which at least one axis is poorly reconstructed (be it due to spurious islands, roundish image from low-impact or just too few usable pixels)

Do you have an idea why this would happen at a particular energy though?

It is actually a rather broad peak, but indeed there must be some reason for such a funny shape. That it gets better at ~TeV and above may perhaps be understandable... larger showers (fewer fluctuations), and likely no low impact involved (hence no roundish image) because only 2 LSTs saw the image, and at such energies if one LST has low impact at least other two would have seen the shower. What is quite hard to understand is that the curve goes down at low energies.

[vuillaut]

Are these point source on-axis events?

[moralejo]

Yes, they are.

[HealthyPear]

Hi, everyone - I'll try to answer and follow your replies, but this weekend I am travelling back home, so I cannot launch or check stuff.

@moralejo : for the moment I couldn't improve much over the code left by Julien, so I don't know how the 68% containment for N=2 after optimization - I saw the complete one (meaning all telescopes) and even if the bump can be seen it has a lower amplitude.

@vuillaut : sometime ago I launched a full-array La Palma N production to see what changes when adding the MSTs; in the end I see 2 bumps, one at ~100 GeV and the other ~3 TeV. As abelardo mentions they are not very narrow, but they are peaked approximately on the central energy value corresponding to the best sensitivity of each camera. An answer to your question would then be simply that the bump is alway there, it's just that the subset of telescopes involved are more sensitive around those energies, and so events at neighbour energies spill into the region of what becomes a bump.

I asked Abelardo (and I will also check ASAP) what happens if I select 2 specfic telescopes (I expect no bump - should be like MAGIC): this is because I think that, in combination with the energy/sensitivity factor there is also a second one related to a selection effect. When we have an array and we select only a subset of pairs, we are intrinsically selecting only the worst events (either very far, or very faint, or a combination of both).

Also, it appears that in the current ctapipe reconstruction, this happens pair-wise via a weight applied on each pair of images (as it should be - better pairs of images get more trust in the reconstruction).
Problem is that apparently when we select N=2 there is always 1 pair, thus no weight, which means that the only reconstruction available is the one that will be used, regardless of the quality (after image cuts survival, of course).

An update for Lugano: I have implemented the MARS 1st pass tailcut in a local version of ctapipe (I will send the pull request next week) + protopipe (this to run things fast) - also we found a bug in the use of the tailcut. I don't think that this will magically solve the problem, but it will surely make things better.

[moralejo]

Below are a few of the distributions from which the sigma_68 are calculated (for the case with no g/h cuts, no image-angle cut). The value delta_angle is the angular distance from true to reconstructed direction. There is nothing especially weird about the distributions (nicely unimodal). And the change from 25 to 160 GeV (~ 0.34 to 0.45 deg in sigma_68, i.e. ~ -0.47 to -0.35 in the x-axis) is rather un-spectacular.

[moralejo]

@HealthyPear: Below is the 2-LST array result (La Palma Prod3b's tel id's 5 and 6, which correspond to LST-1 and LST-2 in the actual array). Pointing south, everything else same as above. The "bump" in angular resolution is still there, and is still due to the least gamma-like events.

Note that MARS-based CTA analysis uses simple image axes crossing for 2-telescope events (and a weighted multi-telescope minimization for >2 telescopes events). In the MAGIC analysis, instead, (for any number of telescopes in an event) a Random-forest-based Disp method is used telescope-wise, then an average direction is calculated. Note also that we never plot angular resolution before gamma selection cuts (which make the bump disappear).

[maxnoe]

I am closing this as this seems to be an inherent property of arrays with different telescopes with this reconstructor and not a bug.