Direction Reconstruction Problems on Diffuse Gammas

[kbruegge]

Hello Everyone.

I'm having trouble reconstructing the direction on diffuse gammas. Here is the code I use. I basically took it from the example scripts. The code below is executed for each event.

    reco = HillasReconstructor()

    features = {}
    params = {}
    pointing_azimuth = {}
    pointing_altitude = {}

    for telescope_id, dl1 in event.dl1.tel.items():
        camera = event.inst.subarray.tels[telescope_id].camera

        mask = tailcuts_clean(camera, dl1.image[0], boundary_thresh=5, picture_thresh=10)
        telescope_type_name = event.inst.subarray.tels[telescope_id].optics.tel_type

        image = dl1.image[0].copy()
        image[~mask] = 0

        h = hillas_parameters(
            camera,
            image,
            container=True
        )

        moments = MomentParameters(size=h.intensity, cen_x=h.x, cen_y=h.y, length=h.length, width=h.width, psi=h.psi )
        params[telescope_id] = moments
        pointing_azimuth[telescope_id] = event.mc.tel[telescope_id].azimuth_raw * u.rad
        pointing_altitude[telescope_id] = ((np.pi/2) - event.mc.tel[telescope_id].altitude_raw )* u.rad # this is weird to say the least. 

    reconstruction = reco.predict(params, event.inst, pointing_azimuth, pointing_altitude)
    d = {
        'mc_alt': event.mc.alt.rad,
        'mc_az': event.mc.az.rad,
        'mc_core_x': event.mc.core_x.value,
        'mc_core_y': event.mc.core_y.value,
        'mc_height_first_interaction': event.mc.h_first_int.value,
        'alt_prediction': reconstruction.alt.si.value,
        'az_prediction': reconstruction.az.si.value,
        'core_x_prediction': reconstruction.core_x.value,
        'core_y_prediction': reconstruction.core_y.value,
        'h_max_prediction': reconstruction.h_max.value,
    }
    return d 

This seems to work fine on point-source gammas. However for diffuse gammas it fails. See below a plot which shows the distribution of the distances (in degrees) between the mc_az, mc_alt and az_prediction, alt_prediction (to be exact the histogram shows astropy.coordinates.angle_utilities.angular_separation)

Any ideas? Do I have a stupid bug?

[TarekHC]

Hi @mackaiver ,

Can you there 2 plots? 2D histogram of mc_alt vs mc_az and alt_prediction vs az_prediction

They may hint to what is really happening...

[kbruegge]

Yes. Good Idea. I'll do so this afternoon. Currently in a meeting.

[kbruegge]

Here you can see that the altitude prediction does not seem to work for me.

[kbruegge]

Here are two more plots to show the same distirbutions broken down to 1D

[TarekHC]

Another test that comes to my mind is actually testing if the azimuthal reconstruction is correct (creating the angular_separation of mc_az and az_prediction).

Regarding on why the predicted altitude is so off, I'm afraid I'm not the one to ask yet...

[kbruegge]

Another test that comes to my mind is actually testing if the azimuthal reconstruction is correct (creating the angular_separation of mc_az and az_prediction).

Yes. Right. I think this is basically what you see here:

So the azimuth reconstruction is not working correctly in my opinion

[vuillaut]

Hi @mackaiver. This issue seems very much related to issue #697. I see it's closed now and you showed a scatter plot that looked fine but have you made a histogram of the impact point errors?

[kbruegge]

I see it's closed now and you showed a scatter plot that looked fine but have you made a histogram of the impact point errors?

Hmm. Not sure I understand. You think its related to wrong reconstrunction of the impact positions?

[TarekHC]

So the azimuth reconstruction is not working correctly in my opinion

Just to confirm this. Can you recalculate those same plots for those events with mc_az close to the center of the FoV (e.g. 0.5 deg)? If they are clearly not well reconstructed then there is a technical issue somewhere with the diffuse gamma reconstruction. If they are well reconstructed, then its probably the reconstruction algorithm.

[kbruegge]

Can you recalculate those same plots for those events with mc_az close to the center of the FoV (e.g. 0.5 deg)?

Here you go. In red you see diffuse gammas as above.

[kbruegge]

have you made a histogram of the impact point errors?

@vuillaut I also made histograms of impact positions. Is that what you mean?

[vuillaut]

@mackaiver yes! but it's always difficult to see the goodness of the results from a diffusion matrix. Would have the same kind of plots done for alt/az but for impact point? Because it would be weird to have a good impact point reconstruction and not a direction one (it's essentially the same thing aside from ground projection).

[kbruegge]

@vuillaut

Would have the same kind of plots done for alt/az but for impact point?

I can do some 1D histograms as well. First the fixed plot form above:

Now here you see the squared distance bettween true and reconstructed postion.

And here the same for each axis:

Colros are the same as above. Red is diffuse gammas

[vuillaut]

So it seems the reconstruction is as good for diffuse as for point-like, meaning that the Hillas parameters (especially ellipse directions) for diffuse are well calculated. There must be some bug related to angles (and only the azimuthal ?). Maybe @TarekHC will have more ideas.

[TarekHC]

Hi @mackaiver and @vuillaut ,

In these plots, are you only cutting in mc_az?

The first thing I would do is replot these distributions in different camera offset bins, cutting on the angular distance between the center of the FoV and the mc_alt and mc_az. What should be happening is that the angular reconstruction gets worse as you get farther away from the center of the FoV. If we don't see that, there must be something quite wrong there.

Another test that you might do is apply quality cuts. Some reasonable ones could be to require a minimum energy (e.g. 50 GeV?) or a minimum multiplicity (4?) as these plots might be currently dominated by ugly 2 telescope triggers with horrible reconstruction (which would NOT be sufficient to reconstruct so horribly the alt...).

[vuillaut]

@TarekHC

these plots might be currently dominated by ugly 2 telescope triggers with horrible reconstruction (which would NOT be sufficient to reconstruct so horribly the alt...)

But wouldn't that be the case for the impact point reconstruction as well? Moreover the average multiplicity is >> 2.

[kbruegge]

What should be happening is that the angular reconstruction gets worse as you get farther away from the center of the FoV.

I checked. Its bad all across the FoV. (Don't have the plots at the moment)

Another test that you might do is apply quality cuts. Some reasonable ones could be to require a minimum energy (e.g. 50 GeV?) or a minimum multiplicity (4?) as these plots might be currently

Heres the resolution vs energy. Which is basically multiplicity.

If we don't see that, there must be something quite wrong there.

There is most definitively something wrong here. It might be the code I wrote or some bug in the reconstruction. Either way maybe @kosack or @watsonjj can have a look at the script I posted above. I'm I using the reconstruction algorithm in the wrong way?

[kosack]

I'd strongly suspect coordinate transform problems - it's a quite important issue we have in ctapipe right now, and nobody is working on it. The definition of all coordinate frames is a mess. Having one rotation wrong somewhere could easily lead to the situation where point-sources come out fine, and diffuse or offset sources do not.

I think all of the current sensitivity curves tino produced were for on-axis sims, which means some of the coord transforms were likely not debugged fully.

[maxnoe]

Corsikas azimuth is 7 degrees off of astropy for la Palma, and I don't know how much for paranal. (Magnetic vs geographic North).

Is that take into account?

[vuillaut]

Hi. I see the issue is still open so I'll add some tests I have done with the file proton_20deg_0deg_run24___cta-prod3_desert-2150m-Paranal-merged.simtel.gz.

Instead of looking at distributions I like to look at single events to see where issues are. Here is one very powerful proton event reconstructed:

Hillas directions are well reconstructed (so, as suspeced, the problem is not here). However, once put together, here is what we get on an array level:

The notebook I used: diffuse_events.ipynb.zip

[kbruegge]

However, once put together, here is what we get on an array level:

Hey thomas,

are you sure you're on latest master? I see in you're notebook you do

pointing_altitude[tel_id] = ((np.pi/2) - event.mc.tel[tel_id].altitude_raw )* u.rad # this is weird to say the least. 

This is not necessary anymore since #758. If I use

pointing_altitude[tel_id] =event.mc.tel[tel_id].altitude_raw * u.rad

I get this:

[kbruegge]

BTW heres my current angular resolution comparing diffuse gammas to pointlike ones. Not sure if this is okay or not.

[vuillaut]

Hi @mackaiver Thank you for pointing out this change!

So interestingly now, there is a converging intersection for the impact point. However, it is still far from the true impact point. That would explain the not so good angular resolution you get for diffuse.

Let me point out a few things:

• angular resolution worsen with energy, which is unexpected.
• angular resolution looks good at low energy (at least closer to what we expect).
• these two points together could be explained by a difference between LST and other telescopes, probably in the coordinates system, or camera rotation.
• this is somewhat confirmed by the LST arrows converging to a different point than the rest of the array in plot 2

I'll do some tests with the updated script.

[moralejo]

On 3 August 2018 at 14:43, Thomas Vuillaume notifications@github.com wrote:

Hi @mackaiver https://github.com/mackaiver Thank you for pointing out this change!

So interestingly now, there is a converging intersection for the impact point. However, it is still far from the true impact point. That would explain the not so good angular resolution you get for diffuse.

A naïve suggestion (likely wrong) : from the example shown by @mackaiver https://github.com/mackaiver, is it possible there is a simple error in the coordinates (90 degree rotation)? Corsika's x-axis points north (which one would usually take as the +y direction). The error in impact point looks pretty much like a 90-deg rotation (on the other hand the mismatch between LSTs and the other arrays might just be due to distant showers and images close to the edge of some LSTs FoV).

--

Abelardo Moralejo Olaizola Institut de Física d'Altes Energies

Tel : +34 931641662 Fax: +34 935811938 Avís - Aviso - Legal Notice - (LOPD) - http://legal.ifae.es

-- Avís - Aviso - Legal Notice - (LOPD) - http://legal.ifae.es http://legal.ifae.es/

[kpfrang]

Hi all, I've been working on using lookup tables for the weighting of the telescopes in the direction reconstruction like @moralejo did before with MARS right here in slide 12.

For each image, we calculate the distance of closest approach (DCA) as shown in the plot below:

We fill lookup tables as a function of size and ratio of width / length for each camera type. Each row in the following plot represents one camera and the columns show different offset bins. A darker color in this plot means that a telescope in this bin in average has a positive impact on the direction reconstruction.

As expected, you see some worsening of the DCA for higher off-axis angles, however the values definitely look reasonable compared to the point source (left column). From this we probably can say that the parametrization of the images look reasonable and that the problem with the diffuse reconstruction is inside the HillasReconstructor.

[thomasgas]

This should be solved after #830? Did you have changes to try it @mackaiver?

[maxnoe]

@kbruegge that was not a problem in your phdthesis anymore, right?