Coma aberration

[aleberti]

Changes to calculate hillas parameters taking into account the coma aberration correction.

[aleberti]

Hi @jsitarek , please have a look at the changes I made to apply the coma aberration correction. For the moment I modified only magic_calib_to_dl1, and running it on some Crab calibrated data to check that the comparison with MARS is now fine.

About lst1_magic_mc_dl0_to_dl1 script, do you think that should be changed as well? I see that simtel files have the following optics definition for MAGIC:

OpticsDescription(name=UNKNOWN, equivalent_focal_length=17.00 m, num_mirrors=-1, mirror_area=235.24 m2)

so it does not use the effective focal length as we do in ctapipe_io_magic (equivalent_focal_length=17*1.0713). So now in principle for simtel_array MCs no correction is applied.

That being said, I will also modify the script to calculate stereo parameters, because now it assumes hillas parameters in the CameraFrame, while now they are saved in the TelescopeFrame.

[YoshikiOhtani]

Hi @aleberti, thank you for the changes. So it saves the length and width parameters with the unit of degree, is it correct? Then we anyway need to modify the lst1_magic_mc_dl0_to_dl1.py because it calculates the Hillas parameters in the CameraFrame.

[aleberti]

Hi @aleberti, thank you for the changes. So it saves the length and width parameters with the unit of degree, is it correct? Then we anyway need to modify the lst1_magic_mc_dl0_to_dl1.py because it calculates the Hillas parameters in the CameraFrame.

yes, correct, see also my comment to Julian. I need to know if I have to use the same effective focal length. Finally, I already changed the stereo reco script in order to use the HillasParametersContainer instead of CameraHillasParametersContainer

[jsitarek]

Hi @aleberti

Thank you for this implementation and sorry for a late feedback.
As far as I understand hillas_reconstructor can work fine both on Hillas parameters in telescope or camera frame:
https://github.com/cta-observatory/ctapipe/blob/d2544512f8cc79f7f655cb17dd7d23fa195a828b/ctapipe/reco/hillas_reconstructor.py#L339
So the new code should provide exactly the same stereo parameters as the old one, have you checked this?

The problem with 17m focal length seems to be already understood by @YoshikiOhtani 's study - without adding an explicit option it was using nominal rather than effective focal length.
If the effective focal length is used, then the stereo reconstruction should work fine no matter if computed from m or deg Hillas container.
If by default ctapipe is not taking into account the coma aberration factor into account in the calculation of the camera frame, I think it is best to keep it like this. We will have different values in mm than MAGIC MARS files, but this is just a scaling factor that can be taken into account in the comparisons. This will only be important if we read in MARS S or Q files. Actually in this sense it is better to have those values in DL1 containers in degrees - since either way a conversion have to be done then when going from MARS files it is less prone to errors.

and definitely +1 from me on modifying lst1_magic_mc_dl0_to_dl1.py, otherwise RFs will not work any more if it starts comparing the values in deg and in m between data and MC.

[YoshikiOhtani]

Hi, ok so let's update lst1_magic_mc_dl0_to_dl1.py to computes the Hillas parameters in a telescope frame, but then which focal length should we use? In my opinion, since LST does not use the effective one and the current ctapipe selects the nominal one by default, it would be better to use the nominal one also for MAGIC. If everyone agrees, I would like to push a commit to ctapipe_io_magic to implement the option "focal_length_choice" also in the MAGICEventSource, which by default gives the "nominal" one. Then when comparing with MARS results one can change the option to "effective". Please let me know what you think, thanks!

EDIT: ctapipe v0.12 that we are using now gives the "nominal" focal length by default, but the most latest version gives the "effective" one, so at some point we will also change the default setting of MAGICEventSource to "effective".

[aleberti]

Ok, so I modified lst1_magic_mc_dl0_to_dl1 in order to use the effective focal lenght when calling EventSource. Also, I modified again lst1_magic_stereo_reco, so that the Hillas container is constructed in the right way depending on how it was saved in the DL1 step. This is because at the moment in LST real data now the Hillas parameters are still calculated in the CameraFrame, with only the conversion to deg of width and length (which is undone in the coincidence script).

[YoshikiOhtani]

OK thanks @aleberti. My original plan was that I will change the coincidence script so that it recomputes the length and width in the degree unit with the effective focal length, otherwise there is the mismatch of the focal length values used in real and MC data. Let me push this commit to this branch later and then it would be great if you could review the change.

[aleberti]

Sure, go ahead!

[YoshikiOhtani]

Hi, so I pushed a commit which recomputes the LST length and width parameters with the effective focal length in case they are originally computed with the nominal one. The nominal and effective values are copied from the LST configuration file uploaded in the lst-sim-config repository. Please have a look at this change @aleberti and @jsitarek if it is fine for you.

[YoshikiOhtani]

Sorry again, but I pushed a new commit because I realized that I should also modify the LST-1 subarray description with the effective focal length which will be used for the stereo reconstruction.

[jsitarek]

Hi, I think there is still something wrong.
I tried running data and MCs with this version.
MAGIC Y files that I reprocess have Hillas parameters in units of deg (they used to be in m). Then after the coincidence script LST1 parameters are in m and MAGIC still in degree, which is very messy, but maybe expected (??) - actually I'm not sure if it is technically possible to do it properly with the QTables, because for different tel_id the same column will be sometimes in one unit and sometimes in another.
Finally the chain fails on stereo reconstruction not finding neither 'x' or 'fov_lon':

Traceback (most recent call last):
File "./lst1_magic_stereo_reco.py", line 350, in
main()
File "./lst1_magic_stereo_reco.py", line 341, in main
stereo_reconstruction(args.input_file, args.output_dir, config, args.magic_only)
File "./lst1_magic_stereo_reco.py", line 236, in stereo_reconstruction
x=u.Quantity(df_tel['x'].iloc[0], u.m),
File "/home/julian.sitarek/anaconda3/envs/magic-lst1/lib/python3.7/site-packages/pandas/core/frame.py", line 3458, in getitem
indexer = self.columns.get_loc(key)
File "/home/julian.sitarek/anaconda3/envs/magic-lst1/lib/python3.7/site-packages/pandas/core/indexes/base.py", line 3363, in get_loc
raise KeyError(key) from err
KeyError: 'x'

[aleberti]

Hi, I think there is still something wrong. I tried running data and MCs with this version. MAGIC Y files that I reprocess have Hillas parameters in units of deg (they used to be in m). Then after the coincidence script LST1 parameters are in m and MAGIC still in degree, which is very messy, but maybe expected (??) - actually I'm not sure if it is technically possible to do it properly with the QTables, because for different tel_id the same column will be sometimes in one unit and sometimes in another. Finally the chain fails on stereo reconstruction not finding neither 'x' or 'fov_lon':

Traceback (most recent call last): File "./lst1_magic_stereo_reco.py", line 350, in main() File "./lst1_magic_stereo_reco.py", line 341, in main stereo_reconstruction(args.input_file, args.output_dir, config, args.magic_only) File "./lst1_magic_stereo_reco.py", line 236, in stereo_reconstruction x=u.Quantity(df_tel['x'].iloc[0], u.m), File "/home/julian.sitarek/anaconda3/envs/magic-lst1/lib/python3.7/site-packages/pandas/core/frame.py", line 3458, in getitem indexer = self.columns.get_loc(key) File "/home/julian.sitarek/anaconda3/envs/magic-lst1/lib/python3.7/site-packages/pandas/core/indexes/base.py", line 3363, in get_loc raise KeyError(key) from err KeyError: 'x'

Ok, then I guess we have to go back to the other idea, that is use effective focal length but do calculations in the CameraFrame, then in principle the stereo reconstruction will take into account the coma aberration correction (if I understand correctly the code). Would this be ok?

[YoshikiOhtani]

Hi @jsitarek and @aleberti,

MAGIC Y files that I reprocess have Hillas parameters in units of deg (they used to be in m). Then after the coincidence script LST1 parameters are in m and MAGIC still in degree, which is very messy, but maybe expected (??)

The commit dbc55f6 should convert the LST length and width to degrees with the effective focal length and so it is strange for me - I guess you just missed the commit?

Finally the chain fails on stereo reconstruction not finding neither 'x' or 'fov_lon'

Yes, now the cog_x and cog_y are saved by different names for MAGIC and LST, "fov_lon" and "fov_lat" for MAGIC, and "x" and "y" for LST-1, which are removed from the data frames since the coincidence script keeps only the common (name) parameters by default. As Julian mentioned in the email, lstchain converts only length and width from meter to degree, not for cog_x and cog_y, so actually there are differences not only the name but also the unit between MAGIC and LST-1. So one solution is to convert them in the coincidence script to degrees with the effective focal length, and rename them to "fov_lon" and "fov_lat". However, I actually agree with what Alessio mentioned:

Ok, then I guess we have to go back to the other idea, that is use effective focal length but do calculations in the CameraFrame, then in principle the stereo reconstruction will take into account the coma aberration correction (if I understand correctly the code). Would this be ok?

Yes, I think this is simpler than the other solutions. I think keeping the Hillas parameters with the unit of meter does not make any problems unless we use the unit for both real and MC data, and actually I don't see any advantage to keep them in degrees. For comparing with MARS results, one just needs to convert the parameters to degrees with the effective focal length. Then the important point for the analysis is to use the effective focal length for the stereo reconstruction, which can be done even if the Hillas parameters are computed in the unit of meter.

Please let me know what you think, especially if you have a strong opinion to keep the parameters in degrees, thanks!

[YoshikiOhtani]

Sorry but since I don't expect a frequent communication due to the shift, I already pushed some commits that I intended. So now all the Hillas parameters are calculated in the CameraFrame as it is done originally, but the stereo reconstruction is performed with the effective focal length. The original purpose of this pull request is to compute the MAGIC Hillas parameters in the TelescopeFrame, but it causes some mess as Julian experienced about the frames and the units, since lstchain calculates the parameters in the CameraFrame. I believe the current change is the most simplest way so far, but please let me know what you think. If you want to keep the Hillas parameters in degrees, I think the option is

• LST real data: convert the Hillas parameters computed in the CameraFrame to degrees with the effective focal length, and rename some parameters in the coincidence script
• MAGIC real data: compute the Hillas parameters in the TelescopeFrame with the effective focal length
• MC data: compute the Hillas parameters in the TelescopeFrame with the effective focal length

[jsitarek]

Hi @YoshikiOhtani
fully agree to simply keep it in units of 'm', especially that LST starts with it
concerning commit dbc55f6 your are right, sorry, once I saw that there is something inconsistent I was playing around and I must have made a mistake checking a wrong file.

[jsitarek]

I'm running now with something that should be compatible with this version, but I really do not understand this part of the code.
I think LST starts in units of 'm' and not 'deg', and here we apply a conversion that will effectively scale down W and L by a factor of 2. And if there is some conversion from m to deg done earlier on and we need to convert it back, why it is only W&L and not the x, y, r, slope, kurtosis, skewness, ... ?
I tried already to remove those two lines, but this has the very bad sensitivity that I've send around, so I must be really missing something here.

[YoshikiOhtani]

Hi Julian, please find the following function that lstchain uses for computing the Hillas parameters:
https://github.com/cta-observatory/cta-lstchain/blob/7b967de9349a196c98d3832b5421b591fd25d926/lstchain/reco/r0_to_dl1.py#L133

Here the input camera geometry is defined as the CameraFrame, so the length, width and etc computed by the hillas_parameters function should have the unit of meter. However, there is a conversion of the unit from meter to degree only for length and width (and their uncertainties) with the _camera_distance_to_angle function afterwards. This is why I implemented the lines to get them back to meter only for length and width, not for the other parameters (cog_x, cog_y, etc...).

[YoshikiOhtani]

Hi @jsitarek, thank you for your comments, good that we all agreed with this solution. I just pushed a minor commit for refactoring the code but it does not harm anything. So if you are fine please give an approval, then we will merge this branch to the master.

[jsitarek]

Hi Yoshiki, Thanks for the explanations. Now it adds up to a consistent behaviour. Then the current code 8n the PR should be perfectly fine, you undo the conversion using the old focal length and then update the focal length to effective value. I am curious why in lstchain it is done like this. Julian pt., 17 cze 2022, 15:42 użytkownik Yoshiki Ohtani ***@***.***> napisał:

…

***@***.**** commented on this pull request. ------------------------------ In magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py <#72 (comment)> : > @@ -143,19 +147,24 @@ def load_lst_data_file(input_file): inplace=True, ) + # Read the subarray description: + subarray = SubarrayDescription.from_hdf(input_file) + # Change the units of some parameters: optics = pd.read_hdf(input_file, key='configuration/instrument/telescope/optics') focal_length = optics['equivalent_focal_length'][0] event_data['length'] = focal_length * np.tan(np.deg2rad(event_data['length'])) # [deg] -> [m] Hi Julian, please find the following function that lstchain uses for computing the Hillas parameters: https://github.com/cta-observatory/cta-lstchain/blob/7b967de9349a196c98d3832b5421b591fd25d926/lstchain/reco/r0_to_dl1.py#L133 Here the input camera geometry is defined as the CameraFrame, so the length, width and etc computed by the hillas_parameters function should have the unit of meter. However, there is a conversion of the unit from meter to degree only for length and width parameters (and their uncertainties) with the _camera_distance_to_angle function afterwards. This is why I implemented the lines to get them back to meter only for length and width, not for the other parameters (cog_x, cog_y, etc...). — Reply to this email directly, view it on GitHub <#72 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH36DYPV2INQ5GH4236QSPDVPR6ELANCNFSM5XY7FUJQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>