Hillas parameters in telescope frame

[LukasNickel]

Related discussions about this can be found in:

• Calculation of Hillas Parameters in the Nominal System #674
• Calculate Hillas Parameters in TelescopeFrame #1090
• Telescope / NominalFrame for reconstruction / visualization #1315
• Direction reconstruction (also) from Telescope Frame #1408

This changes the calculation of the image parameter for the dl1 output to be performed in the nominal frame.
For this, new containers have been added. This is kind of the brute force way. Since it is only relevant for a few container fields (x,y,length, width...), we could also make their unit configurable. However, that might lead to confusion, because the field names do not reflect that difference. Or let one container have both fields of which only one is filled?
Side note: Some code (e.g. calculation of the concentration) explicitly names things x and y, which I did not change.
So maybe a naming of the coordinate axes, that does not indicate the use of a certain coordinate frame could be a way out?

All of these changes are written out by the stage1-tool. Coordinate transformations happen in the ImageProcessor.
Since new frames are constructed for each event, this might be slower. I do not know how to access all pointings at once without iterating over the events though.

Regarding the change of the data model version:
We probably want to collect test files of older versions in order to make sure, that we do not break their functionality.
Is adding dl1 files to ctapipe-extra the way to go or is there another solution?

---

Not longer relevant (keeping it for context):

The more I think about this, the less certain I am if we even want to calculate features in the nominal frame or if a later conversion makes more sense. Anyway, lets discuss this.

Status:
For now this adds a new method and a new container for the calculation of the parameters in the nominal frame.
Since conversion from the camera geometry to the nominal frame requires the telescope pointings, I figured it would be easiest to perform this transformation outside and call the function with a SkyCoord containing the pixel positions. After writing the code I learnt that the camera geometry also has a transform_to method, so thats another possibility. Would make the API clearer but adds some overhead to the calculation.

There is of course a lot of duplicated code right now since the performed steps are pretty much identical. So if we decide to follow up on this approach, it would be beneficial to put some of that code into an extra function.
Some feature (Intensity/Skewness/...) are identical. If we want to include them in the NominalHillasParametersContainer as well probably depends on whether we want both sets of parameters in stage1 or only one of them.

With the coordinate frames at hand, it is possible to convert from one set of parameters to the other and the results should match.

Questions / Points to discuss

• Is this how we want to do it or do we keep the calculation in the camera frame and convert afterwards?
• The current image processor does not seem to be aware of the telescope pointing (?). This can of course be easily changed but might be an argument towards a calculation in the camera frame and later transformation so long as there will not be more image features, that should be calculated in the nominal frame. (Timing?)
• Do we want to keep both containers or would this replace the calculation in the camera frame? -> Should there be duplicated features like intensity or just the "new" ones? What should be part of the stage-1 tool?

Images, Cog, length and width in the nominal frame:

[HealthyPear]

The more I think about this, the less certain I am if we even want to calculate features in the nominal frame or if a later conversion makes more sense.

From the point of view of ctapipe.reco.HillasReconstructor it would probably be overkill to calculate the parameters in this frame, for 2 reasons:

• the only practical gain to go from CameraFrame to any sky-based frame is to remove the optics from the equation and better compare the images afterward (this is at the core of Direction reconstruction (also) from Telescope Frame #1408 , not only because of the comparison with CTAMARS)
• the only use of the NominalFrame for non-parallel pointing, either from mispointing or from divergent pointing or a mix of the 2 - in the former case we should apply anyway pointing corrections, so the effect would be small, in the latter ctapipe.reco.HillasReconstructor already takes care of this by projecting from the sky to a "fake" array-based CameraFrame, thus nullifying the divergent pointing.

Of course, I see no problem in adding the possibility of having the parameters in yet another frame, but that would complicate #1408 and it's encapsulation in stage2 and I think merging that before would make things clearer also for this PR.

Anyway here you just add the possibility to have the parameters in that frame, so it doesn't harm the rest.

[LukasNickel]

Yeah, the HillasReconstructor does not really care about this.
It is more about what we want to store in dl1 files rather than the dl2 calculation.

[maxnoe]

@LukasNickel

What we discussed was something else. The task for the next ctapipe release should be to have the DL1 in degrees.
I am no expert in the divergent pointing analysis, so I don't know if it makes sense to go to NominalFrame.

For now, let's keep it simple.

• Do not copy the code, just make sure it works when pix_x and pix_y are in degrees.
• Use a different container so the difference is transparent
• Change the ImageProcessor to do this by default

[LukasNickel]

* Do not copy the code, just make sure it works when `pix_x` and `pix_y` are in degrees.

* Use a different container so the difference is transparent

* Change the ImageProcessor to do this by default

Ok, lets do it like this.

[maxnoe]

Since we want to switch to the hillas parameters in degree by default I think it would make sense to actually make the HillasParametersContainer the one with the stuff in degrees and have a CameraHillasParameters container for the one with unit meters.

@kosack @HealthyPear opinions?

[HealthyPear]

Since we want to switch to the hillas parameters in degree by default I think it would make sense to actually make the HillasParametersContainer the one with the stuff in degrees and have a CameraHillasParameters container for the one with unit meters.

@kosack @HealthyPear opinions?

I like the idea!
Should I change this directly into #1408 or we do it here after the merge?

[maxnoe]

• Use telescope frame
• Make it default in the image processor
• Make the normal containers without prefix the ones with stuff in degrees
• Keep the old ones with a prefix Camera....
• Remove any hard coded sites / array pointings /

[kosack]

A few comments:

• I don't see any reason to keep two copies of the containers. Just convert everything spatial to TelescopeFrame, and drop the old CameraFrame parts. (caveat: see second to last point)
• There is a very good reason not to use the NominalFrame (other than it is not used by the HillasReconstructor), and that is that we may update the astrometric (pointing) corrections, which changes the NominalFrame. It is nice to not have to redo DL1 if those change. Also the center position of the nominal frame is not well defined (by default it is the barycenter of the various telescope pointings, but that may be changed/optimized later). If we store in TelescopeFrame, it's easy to add any corrections before reconstruction.
• In principle you don't even need two containers if you want to retain the ability to compute hillas parameters in any frame. The original idea was that hillas.x could be in meters or degrees. The only problem is just that the containers so far force a unit. If we want to support both cases we could make that container have a constructor argument that specifies the unit or something like that. It could be useful to support the CameraFrame as well, as the impact point reconstruction is done in the camera frame - see all the tricks that Michele had to do to have the psi angle in both Telescope and Camera frame without breaking things..., but it's not really necessary
• Another option is to just make a conversion function that converts lengths and angles from one frame to another, and then you can convert HillasParameters (and any other spatial parameters) without having to change the CameraGeometry and re-compute them. That could get rid of the need to store the psi angle in two frames as in Direction reconstruction (also) from Telescope Frame #1408

[LukasNickel]

In principle you don't even need two containers if you want to retain the ability to compute hillas parameters in any frame. The original idea was that hillas.x could be in meters or degrees. The only problem is just that the containers so far force a unit. If we want to support both cases we could make that container have a constructor argument that specifies the unit or something like that.

I feel like this makes it harder to keep track of the units and coordinate frames. The reconstructor could then not rely on getting a HillasParametersContainer alone, but would also need to check if its in the correct frame. Also x/y are heavily connected to the camera frame in my opinion. If its still worth it in order to keep a cleaner looking container structure, I do not really know right now.

It could be useful to support the CameraFrame as well, as the impact point reconstruction is done in the camera frame - see all the tricks that Michele had to do to have the psi angle in both Telescope and Camera frame without breaking things..., but it's not really necessary

Do you mean, that you want to save results in both coordinate frames in the dl1 output?

Another option is to just make a conversion function that converts lengths and angles from one frame to another, and then you can convert HillasParameters (and any other spatial parameters) without having to change the CameraGeometry and re-compute them. That could get rid of the need to store the psi angle in two frames as in #1408

General conversion functions for the dl2 reconstructors could be handy, agreed. If we agree on a list of functions, I can add them here, in a new pr or we see it as part of #1408.

[LukasNickel]

Regarding the failing test:
Travis is consistently calculating 0.06638165 deg as the width, while my local tests run just fine with 0.06022414 deg.
Why is there a difference at all? Shouldn't the toymodel produce the same image and the transformation + hillas calculation be identical?

[maxnoe]

Will be fixed by merging #1589 or removing your cache

[HealthyPear]

I also like the idea of making the HillasContainer a bit more adaptive (I tried to propose this during #1408 in fact).

If there's a way to do this It would be cool to build a HillasContainer with a frame object and a unit, so even if unit is degrees one can differentiate between TelescopeFrame and NominalFrame.

Also because in #1408 I didn't change the Hillass containers, I just write to them in a frame, but I think that for a dummy event the units would be wrong (the values instead would be always either 0 or NaN)

[maxnoe]

@HealthyPear No, please don't! We should keep containers simple without any logic other then storing / validating data. I think they are too complex as it is.

[HealthyPear]

@HealthyPear No, please don't! We should keep containers simple without any logic other then storing / validating data. I think they are too complex as it is.

Then if we want to store Hillas parameters in more than one frame we have to create at least 2 of them, one for meters and one for degrees.
I don't see many other alternatives...

[LukasNickel]

• Added an option to use the camera frame for the ImageProcessor. Its not available from the stage1 tool, is that a thing we want/need? I wasnt sure how that would work with the changed datamodel version.
• Merged the master and resolved some conflicts. Due to the changes in the datamodel in the master (indices lookup) this change now leads to version 2.1 instead of 1.3 of the datamodel

[kosack]

This should be a Traitlet, not a construction option. Otherwise it is not configurable. The default should be set to True.

[kosack]

That will not change the data model - it is only a configuration option, and the current data model supports both units in the output. However, we may in the future deprecate the CameraFrame one and clean up all that code as I said before, but now we want to be able to run the full analysis twice on a large dataset and ensure they results are compatible with both frames.

[maxnoe]

While we are at it, please also remove the is_simulation option. It should not be needed. Code that needs to check that should just check if event.simulation is None

[maxnoe]

Sorry, not really connected. Will just make a PR

for some reason clicking "request changes" didn't work

[kosack]

I did a quick check running in both the telescope (using this version) and camera frame (using the master version) and comparing the differences. It all looks good. There are some slight differences, but this may be due to rounding error. We may need to investigate a bit though, but that can come later. For example, not sure why delta alt is not centered on 0: there is a small bias

[LukasNickel]

Sorry, thought I pushed that yesterday.

@kosack What did you use for that comparison? Diffuse or pointlikes?

[HealthyPear]

Sorry, thought I pushed that yesterday.

@kosack What did you use for that comparison? Diffuse or pointlikes?

It shouldn't matter, but better check with both.

[kosack]

@kosack What did you use for that comparison? Diffuse or pointlikes?

point-like. This was just a quick test to make sure it was not crazy. We should do a proper benchmark later. The errors are all quite small.

[HealthyPear]

@kosack What did you use for that comparison? Diffuse or pointlikes?

point-like. This was just a quick test to make sure it was not crazy. We should do a proper benchmark later.

we should have a small diffuse test file with 10 showers if I recall correctly

[LukasNickel]

It shouldnt, I am just curious especially if there is a bias

[LukasNickel]

we should have a small diffuse test file with 10 showers if I recall correctly

I dont think 10 showers will suffice for this kind of benchmark

[kosack]

@HealthyPear There's already a unit test for this for HillasReconstructor, test_CameraFrame_against_TelescopeFrame so no need here, I just wanted to be sure the mechanism all works with ctapipe-process.