Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use effective focal length to compute src-dep parameters #1131

Merged
merged 13 commits into from
Jul 23, 2023

Conversation

SeiyaNozaki
Copy link
Collaborator

Related to #1048
This PR allows to using the effective focal length (instead of the equivalent focal length:28 m) to compute source-dep parameters. Due to the version incompatibility, it raises an error when reading subarray table using lstchain v0.10 + DL1 data processed by lstchain v0.9. So, if it fails to read subarray table to get the effective focal length, it just uses the standard LST effective focal length.

tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length
except OSError:
print("subarray table is not readable because of the version inompatibility.")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should use logging system, not prints

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done!

@codecov
Copy link

codecov bot commented Jul 3, 2023

Codecov Report

Patch coverage: 70.52% and project coverage change: -0.09 ⚠️

Comparison is base (78af94b) 73.88% compared to head (836030b) 73.79%.

Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1131      +/-   ##
==========================================
- Coverage   73.88%   73.79%   -0.09%     
==========================================
  Files         125      124       -1     
  Lines       12569    12577       +8     
==========================================
- Hits         9286     9281       -5     
- Misses       3283     3296      +13     
Impacted Files Coverage Δ
lstchain/scripts/tests/test_lstchain_scripts.py 99.65% <ø> (-0.02%) ⬇️
lstchain/tests/test_lstchain.py 97.48% <ø> (ø)
lstchain/scripts/lstchain_dl1_to_dl2.py 76.61% <55.55%> (-2.21%) ⬇️
lstchain/scripts/lstchain_mc_rfperformance.py 56.75% <55.55%> (-2.07%) ⬇️
lstchain/reco/dl1_to_dl2.py 80.35% <74.02%> (-0.37%) ⬇️

... and 1 file with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

@moralejo
Copy link
Collaborator

moralejo commented Jul 5, 2023

@SeiyaNozaki I am not 100% sure this is correct.
I see the effective focal length is used to compute the expected source position (in the call to get_expected_source_pos).

However, if you use the effective focal length for computing the image parameters, then it is like having no aberration in the reco params (e.g. reco direction). Hence the expected source position (to be used as reference for the src-dep parameters computations) should be calculated with the nominal focal length, right?

@moralejo
Copy link
Collaborator

Could you have a look at this, @SeiyaNozaki?

@SeiyaNozaki
Copy link
Collaborator Author

Just let me summarize the current parameterization and how we should correct aberration effect.

If I understand correctly, hillas parameters are computed in camera frame (unit: m) in ctapipe.

When computing parameters related MC true gamma-ray position (including src_x, src_y) inside add_disp_to_parameters_table, equivalent focal length is used for the conversion of MC gamma-ray arrival direction from sky coordinate to camera frame. So those values are not aberration-corrected.

focal_length = subarray.tel[telescope_id].optics.equivalent_focal_length

source_pos_in_camera = sky_to_camera(df.mc_alt.values * u.rad,
df.mc_az.values * u.rad,
focal,
df.mc_alt_tel.values * u.rad,
df.mc_az_tel.values * u.rad,
)

alpha and dist are computed using x, y, expected_src_x, expected_src_y (all in camera frame), and psi. So those expected source position should be computed after aberration correction (by using effective focal length). It means src_x and src_y (for MC) should also computed after aberration correction

if data_type == 'mc_gamma':
expected_src_pos_x_m = data['src_x'].values
expected_src_pos_y_m = data['src_y'].values

Am I correct...?

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

If I understand correctly, hillas parameters are computed in camera frame (unit: m) in ctapipe.

in lstchain

ctapipe supports both, computing in telescope frame (degrees, with configuration if effective or equivalent focal length should be used for the transformation) and camera frame.

Since version 0.13, computing them in telescope frame is the default for the ctapipe-process tool and since 0.15 the effective focal length is used by default.

We actually plan to remove computation of hillas parameters (and related image parameters) in camera frame in one of the upcoming releases, see:

cta-observatory/ctapipe#2061

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

I think the preferred solution here would but to just switch to computation in telescope frame for everything, not piece-by-piece applying effective focal length corrections to different parts of the code.

@SeiyaNozaki
Copy link
Collaborator Author

what do you think @moralejo ? It could be a big change, but we may need to apply this change at some point.

@moralejo
Copy link
Collaborator

If I understand correctly, hillas parameters are computed in camera frame (unit: m) in ctapipe lstchain.

Ok, so they are in meters, and measured "on the camera" (i.e. the effect of aberration is present, e.g. centroids are slightly further out than they would for ideal optics with the "equivalent focal length" of 28 m)

When computing parameters related MC true gamma-ray position (including src_x, src_y) inside add_disp_to_parameters_table, equivalent focal length is used for the conversion of MC gamma-ray arrival direction from sky coordinate to camera frame. So those values are not aberration-corrected.

def add_disp_to_parameters_table(dl1_file, table_path, focal):

I think there is a problem there. You pass to the function the equivalent focal length. This is used to calculate the true source position on the camera (in m). That is the position for ideal optics (28 m focal), but the computed image parameters are affected by aberration. Either we calculate the true source position modified by aberration, OR we correct the aberration in the image parameters (in this case just the centroid position).

The way it is now, the disp_dx, disp_dy, disp_norm will come from aberration-affected (centroid_x, centroid_y) and aberration_free (source_position_x, source_position_y). Probably the effect is not large for the source-INdependent analysis, because the RF would "learn" to estimate the source position (for ideal optics) from the aberration-affected images.

@moralejo
Copy link
Collaborator

focal_length = subarray.tel[telescope_id].optics.equivalent_focal_length

source_pos_in_camera = sky_to_camera(df.mc_alt.values * u.rad,
df.mc_az.values * u.rad,
focal,
df.mc_alt_tel.values * u.rad,
df.mc_az_tel.values * u.rad,
)

alpha and dist are computed using x, y, expected_src_x, expected_src_y (all in camera frame), and psi. So those expected source position should be computed after aberration correction (by using effective focal length). It means src_x and src_y (for MC) should also computed after aberration correction

If you want to use the parameters in camera frame (=aberration-affected) then, for calculating src-dep params, one should use the "true source position" also affected by aberration (=computed using the effective focal length). But this is not what we do now.

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

Ok, so they are in meters, and measured "on the camera" (i.e. the effect of aberration is present, e.g. centroids are slightly further out than they would for ideal optics with the "equivalent focal length" of 28 m)

I prefer a slightly different point of view, which in the end has the same interpretation however. The centroids in the camera are where they are, independent of the frame, i.e. the pixel values don't change.

Where the aberration comes in is the question, to which point on the sky a pixel corresponds to. And the effective focal length is the more correct answer here.

I think there is a problem there. You pass to the function the equivalent focal length. This is used to calculate the true source position on the camera (in m). That is the position for ideal optics (28 m focal), but the computed image parameters are affected by aberration.

Either we calculate the true source position modified by aberration

Which means using the effective focal length for the transformation from AltAz through TelescopeFrame into CameraFrame

OR we correct the aberration in the image parameters (in this case just the centroid position).

Better than correcting the image parameters after the fact would be already computing them so they are correct, which is what computing them in the TelescopeFrame does, if the effective focal length is used to transform the pixel coordinates from CameraFrame to TelescopeFrame.

In r0_to_dl1 and dl1ab, we would need to transform the camera geometry:

geometry_tel_frame = geometry.transform_to(TelescopeFrame())

and then just use that geometry for the computation of all image parameters.

ctapipe_io_lst always fills the effective_focal_length into the frame of CameraGeometry.

@moralejo
Copy link
Collaborator

moralejo commented Jul 12, 2023

I think the preferred solution here would but to just switch to computation in telescope frame for everything, not piece-by-piece applying effective focal length corrections to different parts of the code.

Seems the cleanest solution, yes, if the aberration correction is applied (i.e. effective focal length is always used) when transforming the image parameters from camera to telescope frame.
The computation of true source direction in telescope frame should have nothing to do with the focal lengths, so there is no danger of confusion in that case.

Still, I think it is better to finish this PR with just the needed corrections, right? Or do you have time @SeiyaNozaki to attempt the major change now?

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

Seems the cleanest solution, yes, if the aberration correction is applied (i.e. effective focal length is always used) when transforming the image parameters from camera to telescope frame.

Technically, we don't transform the image parameters. We transform the pixel coordinates into TelescopeFrame using the effective focal length and then compute image parameters in the TelescopeFrame.

@moralejo
Copy link
Collaborator

Better than correcting the image parameters after the fact would be already computing them so they are correct, which is what computing them in the TelescopeFrame does, if the effective focal length is used to transform the pixel coordinates from CameraFrame to TelescopeFrame.

Yes, I agree.

In r0_to_dl1 and dl1ab, we would need to transform the camera geometry:

geometry_tel_frame = geometry.transform_to(TelescopeFrame())

and then just use that geometry for the computation of all image parameters.

ctapipe_io_lst always fills the effective_focal_length into the frame of CameraGeometry.

Seems simple enough... Would it then be ok to make that part of this PR?

BTW, I think all this still allows running dl1_to_dl2 on existing DL1 files (from v0.9), right? (eff. focal length is hard-coded)

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

BTW, I think all this still allows running dl1_to_dl2 on existing DL1 files (from v0.9), right? (eff. focal length is hard-coded)

Starting from the images yes, the parameters will change though.

@moralejo
Copy link
Collaborator

Ah, right, no problem if we re-do the DL1ab part...

Current DL1b:

  • image parameters in camera frame (in m, obviously "aberration-affected" - apologies for the jargon)
  • for MC: true source position in camera frame, src_x, src_y in m, but computed for "ideal f=28 m optics" (hence should not be "mixed" with the image parameters without correcting one or the other!)

If we go for the change, DL1b would be like:

  • image parameters in telescope frame (in deg, and "corrected for aberration")
  • for MC: true source position in deg in telescope frame src_x, src_y would be (src_fov_lon, src_fov_lat). This is fully consistent with the image parameters.

I think we should go for this change before the next reprocessing (from raw) of the full sample, but I hesitate to go for it straightaway. Perhaps it is better to make sure, with as few changes as possible, that the source-dependent analysis is properly done, to check if the discrepancies we see with source-independent may be related to this aberration issue.

@maxnoe
Copy link
Member

maxnoe commented Jul 12, 2023

image parameters in camera frame (in m, obviously "aberration-affected" - apologies for the jargon)

I think its even worse... length and width are transformed by hand into degrees using the equivalent focal length, the rest of the image parameters is left untouched and is in meters or the related units.

@moralejo
Copy link
Collaborator

Let me summarize my understanding of the current code in this PR:

In get_expected_source_pos:

# For gamma MC, expected source position is actual one for each event
if data_type == 'mc_gamma':
expected_src_pos_x_m = data['src_x'].values
expected_src_pos_y_m = data['src_y'].values

Note that the focal length passed as an argument to get_expected_source_pos is not used at all for MC gammas (only for MC protons). The pre-calculated values data['src_x'] (and src_y) are used for gammas - and they were computed with the equivalent focal length (i.e. 28 m), in r0_to_dl1:

source_pos_in_camera = sky_to_camera(df.mc_alt.values * u.rad,
df.mc_az.values * u.rad,
focal,
df.mc_alt_tel.values * u.rad,
df.mc_az_tel.values * u.rad,
)

Now, this source position computed with f=28 m, i.e. assuming perfect optics, is used to calculate the source-dependent parameters (from image parameters in camera frame, i.e. "affected by aberration"):

src_dep_params = calc_source_dependent_parameters(data, expected_src_pos_x_m, expected_src_pos_y_m)

This is what is not correct in my opinion.
I think the solution (for the source-dependent analysis) is to modify get_expected_source_pos to recalculate the "true gamma direction" using the effective focal length (i.e. introducing the aberration). With a call to sky_to_camera using the effective focal length (i.e. the same as shown above in r0_to_dl1, but with the effective focal)

@moralejo
Copy link
Collaborator

That is for source-dependent analysis. For source-independent the target disp values have the same problem, so that should be corrected as well.

@SeiyaNozaki
Copy link
Collaborator Author

Okay, I agree with @moralejo

Strictly speaking, reco_disp_norm_diff is computed as the difference between dist and reco_disp_norm.

train['reco_disp_norm_diff'] = np.abs(train['dist'] - train['reco_disp_norm'])

So I want just to update all disp parameters (including src_x, src_y) for MC simulation events using effective focal length for source-dependent analysis. Does it make sense?

@moralejo
Copy link
Collaborator

Strictly speaking, reco_disp_norm_diff is computed as the difference between dist and reco_disp_norm.

train['reco_disp_norm_diff'] = np.abs(train['dist'] - train['reco_disp_norm'])

So I want just to update all disp parameters (including src_x, src_y) for MC simulation events using effective focal length for source-dependent analysis. Does it make sense?

Sorry, I don't understand what you mean. If we recalculate (src_x, src_y) with effective focal length then everything is in camera frame, and all source-dep parameters will be correct, right?

@moralejo
Copy link
Collaborator

As far as I could see, all calls to apply_models use the equivalent focal length, 28 m:

def apply_models(dl1,

This means that the conversion from reco event direction in camera frame to Alt Az is done ignoring the effect of aberration. This is not correct, but may be partly compensated by the fact that the RF was trained with the target direction (on camera) computed also with the equivalent focal length - so the RF could somehow "learn" to go from images with aberration to reco directions without it.

We should however fix it everywhere, probably it is more clear if done in this PR.

For consistency we should also change the stored values of src_x, src_y in the MC DL1 files. As of now they are also computed with equivalent f. They will anyway not be used in the pipeline if we go for the changes above - so that we can run dl1-to-d2 on existing DL1 files.

@SeiyaNozaki
Copy link
Collaborator Author

ok! So we need to compute src_x, src_y and other disp parameters including disp_norm and disp_sign using effective focal length, then store those updated values in DL2 data. This update always should be done independent to the analysis method (src-indep vs src-dep), right?

I will implement so.

@SeiyaNozaki
Copy link
Collaborator Author

Strictly speaking, reco_disp_norm_diff is computed as the difference between dist and reco_disp_norm.

train['reco_disp_norm_diff'] = np.abs(train['dist'] - train['reco_disp_norm'])

So I want just to update all disp parameters (including src_x, src_y) for MC simulation events using effective focal length for source-dependent analysis. Does it make sense?

Sorry, I don't understand what you mean. If we recalculate (src_x, src_y) with effective focal length then everything is in camera frame, and all source-dep parameters will be correct, right?

Sorry please forget it... I just wanted to say I want to recalculate not only src_x, src_y, but also all disp parameters (disp_norm, disp_sign etc) which is used for RF training.

@@ -794,11 +806,11 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3
# For proton MC, nominal source position is one written in config file
if data_type == 'mc_proton':
expected_src_pos = utils.sky_to_camera(
u.Quantity(data['mc_alt_tel'].values + config['mc_nominal_source_x_deg'], u.deg, copy=False),
u.Quantity(data['mc_az_tel'].values + config['mc_nominal_source_y_deg'], u.deg, copy=False),
u.Quantity(data['mc_alt_tel'].values + np.deg2rad(config['mc_nominal_source_x_deg']), u.rad, copy=False),
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am confused about this.... you add something called "nominal source_{x,y}" to something called az/alt?

This seems wrong, at least if the operation is correct, than the names.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This parameter is already there, it specifies the relative position (0.4 deg on x-axis by default) to compute expected source position for proton MC.

"mc_nominal_source_x_deg": 0.4,
"mc_nominal_source_y_deg": 0.0,

Just I found the bug, though it didn't affect the final result much

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(0.4 deg on x-axis by default)

What x-axis? In which coordinate system?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep, indeed it may be confusing... It assumes in AltAz coordinates (0.4 deg).. This parameter was introduced here:
#247 (comment)

Copy link
Collaborator

@moralejo moralejo Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep, indeed it may be confusing... It assumes in AltAz coordinates (0.4 deg).. This parameter was introduced here: #247 (comment)

u.Quantity(data['mc_az_tel'].values + np.deg2rad(config['mc_nominal_source_y_deg']), u.rad, copy=False),

It is not correct, camera frame's x is "along the altitude axis", but y is definitely not the azimuth axis....

Also, why do you need this? mc_nominal_source_x_deg and mc_nominal_source_y_deg are already the position on the camera (well, in degrees and without aberration) w.r.t. which you want the src-dep reco for protons.
You just have to transform it to meters

expected_src_pos_y_m = np.tan(np.deg2rad(mc_nominal_source_y_deg)) * effective_focal_length
expected_src_pos_x_m = np.tan(np.deg2rad(mc_nominal_source_x_deg)) * effective_focal_length

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just resolved the conflict resulting from the azimuth range change

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

expected_src_pos_y_m = np.tan(np.deg2rad(mc_nominal_source_y_deg)) * effective_focal_length
expected_src_pos_x_m = np.tan(np.deg2rad(mc_nominal_source_x_deg)) * effective_focal_length

We have the ctapipe coordinate systems for this.

The source position is given in AltAz, correct? or where does this mc_nominal_source_{x,y}_deg come from?

Then:

source_pos = SkyCoord(alt=..., az=.., frame=AltAz())
pointing = SkyCoord(alt=..., az=..., frame=AltAz())
camera_frame = CameraFrame(telescope_pointing=pointing, focal_length=effective_focal_length)

source_camera = source_pos.transform_to(camera_frame)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, the "fake position" for the proton analysis is provided in "camera coordinates" x, y, but in degrees and without aberration. So this converts it to true camera coordinates.

Copy link
Member

@maxnoe maxnoe Jul 19, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, then it's:

source_pos = SkyCoord(fov_lon=..., fov_lat=.., frame=TelescopeFrame())
camera_frame = CameraFrame(focal_length=effective_focal_length)
source_camera = source_pos.transform_to(camera_frame)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, then it's:

source_pos = SkyCoord(fov_lon=..., fov_lat=.., frame=TelescopeFrame())
camera_frame = CameraFrame(focal_length=effective_focal_length)
source_camera = source_pos.transform_to(camera_frame)

Yes, and I think fov_lon = -mc_nominal_source_y_deg and fov_lat = mc_nominal_source_x_deg

It does not matter much, however, because the direction of the "nominal" position chosen for protons is not critical I think.

Comment on lines 70 to 71
print("subarray table is not readable because of the version inompatibility.")
print("Use the effective focal lentgh for the standard LST optics")
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typos: inompatibility, lentgh

Also, better "I will use the effective..." (otherwise it seems it is asking the user for some change )

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"The effective focal length for the standard LST optics will be used"

Comment on lines 74 to 75
src_dep_df = pd.concat(get_source_dependent_parameters(
dl1_params, config, effective_focal_length=effective_focal_length), axis=1)
Copy link
Collaborator

@moralejo moralejo Jul 19, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I already commented about this. In

def get_expected_source_pos(data, data_type, config, effective_focal_length=29.30565 * u.m):

for gamma mc the returned nominal position on the camera is (I think) obtained with nominal focal length (pre-calculated):

if data_type == 'mc_gamma':
expected_src_pos_x_m = data['src_x'].values
expected_src_pos_y_m = data['src_y'].values

Wait, it seems the values src_x, src_y finally used in building the models and in applying them are not directly those from dl1 files, but are recalculated with effective_focal_length, so I think it is probably working.

On the other hand, lstchain_add_source_dependent_parameters.py will keep the src_x, src_y that are found in the dl1 file (calculated with "equivalent f" for existing v0.9 files).

I think that what is misleading is the behaviour of get_expected_source_pos:

get_expected_source_pos(data, data_type, config, effective_focal_length=29.30565 * u.m)

I would expect the line above to always recalculate the source position using the focal length that is passed. But for MC gammas it just returns whatever is already stored in "data"

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

src_x and src_y is updated using effective focal length here (I thought this was done before get_expected_source_pos function, but it is not the case... I will update the code)

dl2 = update_disp_with_effective_focal_length(dl2, effective_focal_length = effective_focal_length)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I said, I think you have to update get_expected_source_pos so that it recalculates the position, not just propagates...

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, it is true. Then I will modify codes to always recalculate the source position inside this function.
This function is only used for source-dep analysis. But, it would be better to recalculate the source position somewhere else also for the source-indep analysis addressing this comment, correct?
#1131 (comment)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, it is true. Then I will modify codes to always recalculate the source position inside this function. This function is only used for source-dep analysis. But, it would be better to recalculate the source position somewhere else also for the source-indep analysis addressing this comment, correct? #1131 (comment)

Yes, for source-indep we need that the target for direction RF (true source position) is also in "true camera coordinates", meaning that it contains the effect of aberration (just like the shower images).

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I see.

For lstchain_add_source_dependent_parameters.py, I'd even remove this script since source-dep parameters are also computed during dl1_to_dl2 step if src-dep parameters don't exist.
This script overwrites the existing DL1 data to add source-dependent data frame, so it doesn't make sense for the current data analysis scheme since we analyze the common DL1 data produced by OSA.

@moralejo
Copy link
Collaborator

Looks fine to me. If I understood correctly, the source-independent analysis is also corrected, by using as target in build_models the true gamma direction in camera coordinates computed with effective focal length. Then, after applying the disp RF, the direction is converted to Alt Az using again effective focal length.

@rlopezcoto rlopezcoto merged commit 584f43b into main Jul 23, 2023
7 of 9 checks passed
@rlopezcoto rlopezcoto deleted the use_effective_focal_length_for_srcdep branch July 23, 2023 13:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants