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.

Sign up for GitHub

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

Jump to bottom

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

Merged
merged 13 commits into from
Jul 23, 2023
Merged

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

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
e31408b
use effective focal length to compute src-dep parameters
SeiyaNozaki Jul 3, 2023
7c146d4
small bug fix, and update the value for pytest
SeiyaNozaki Jul 3, 2023
8543db8
Use logging instead of print
SeiyaNozaki Jul 3, 2023
b9b5b55
update disp parameters using effective focal length to take into acco…
SeiyaNozaki Jul 13, 2023
709219d
small bug fix
SeiyaNozaki Jul 14, 2023
7871b34
small update
SeiyaNozaki Jul 14, 2023
a34fcff
update the calculation of assumed source position for proton MC
SeiyaNozaki Jul 14, 2023
500d7b2
Merge branch 'main' into use_effective_focal_length_for_srcdep
moralejo Jul 19, 2023
fd6b0ca
correct typo
SeiyaNozaki Jul 19, 2023
0945601
recalculate MC gamma source position inside get_expected_source_pos f…
SeiyaNozaki Jul 19, 2023
09dcf48
remove lstchain_add_source_dependent_parameters.py
SeiyaNozaki Jul 19, 2023
4b95835
update for the calculation of proton MC expected source position
SeiyaNozaki Jul 19, 2023
836030b
remove lstchain_add_source_dependent_parameters from index.rst
SeiyaNozaki Jul 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 35 additions & 17 deletions lstchain/reco/dl1_to_dl2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@

from ctapipe.image.hillas import camera_to_shower_coordinates
from ctapipe.instrument import SubarrayDescription
from ctapipe_io_lst import OPTICS

__all__ = [
'apply_models',
Expand Down Expand Up @@ -349,10 +350,18 @@ def build_models(filegammas, fileprotons,
src_dep_df_gamma = get_srcdep_params(filegammas)

else:
subarray_info = SubarrayDescription.from_hdf(filegammas)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length
src_dep_df_gamma = get_source_dependent_parameters(df_gamma, config, focal_length=focal_length)
try:
subarray_info = SubarrayDescription.from_hdf(filegammas)
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!

print("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

src_dep_df_gamma = get_source_dependent_parameters(
df_gamma, config, effective_focal_length=effective_focal_length
)

df_gamma = pd.concat([df_gamma, src_dep_df_gamma['on']], axis=1)

Expand All @@ -362,10 +371,18 @@ def build_models(filegammas, fileprotons,
src_dep_df_proton = get_srcdep_params(fileprotons)

else:
subarray_info = SubarrayDescription.from_hdf(fileprotons)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length
src_dep_df_proton = get_source_dependent_parameters(df_proton, config, focal_length=focal_length)
try:
subarray_info = SubarrayDescription.from_hdf(fileprotons)
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.")
print("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

src_dep_df_proton = get_source_dependent_parameters(
df_proton, config, effective_focal_length=effective_focal_length
)

df_proton = pd.concat([df_proton, src_dep_df_proton['on']], axis=1)

Expand Down Expand Up @@ -536,7 +553,7 @@ def apply_models(dl1,
reg_disp_vector=None,
reg_disp_norm=None,
cls_disp_sign=None,
focal_length=28 * u.m,
effective_focal_length=29.30565 * u.m,
custom_config=None,
):
"""
Expand All @@ -562,7 +579,7 @@ def apply_models(dl1,
cls_disp_sign: string | Path | bytes | sklearn.ensemble.RandomForestClassifier
Path to the random forest filename or file or pre-loaded RandomForestClassifier object
for disp sign reconstruction
focal_length: `astropy.unit`
effective_focal_length: `astropy.unit`
custom_config: dictionary
Modified configuration to update the standard one

Expand Down Expand Up @@ -654,7 +671,7 @@ def apply_models(dl1,
dl2.y.values * u.m,
dl2.reco_disp_dx.values * u.m,
dl2.reco_disp_dy.values * u.m,
focal_length,
effective_focal_length,
alt_tel * u.rad,
az_tel * u.rad)

Expand Down Expand Up @@ -688,7 +705,7 @@ def apply_models(dl1,
return dl2


def get_source_dependent_parameters(data, config, focal_length=28 * u.m):
def get_source_dependent_parameters(data, config, effective_focal_length=29.30565 * u.m):
"""Get parameters dict for source-dependent analysis.

Parameters
Expand All @@ -704,8 +721,9 @@ def get_source_dependent_parameters(data, config, focal_length=28 * u.m):
else:
data_type = 'real_data'

expected_src_pos_x_m, expected_src_pos_y_m = get_expected_source_pos(data, data_type, config,
focal_length=focal_length)
expected_src_pos_x_m, expected_src_pos_y_m = get_expected_source_pos(
data, data_type, config, effective_focal_length=effective_focal_length
)

src_dep_params = calc_source_dependent_parameters(data, expected_src_pos_x_m, expected_src_pos_y_m)
src_dep_params_dict = {'on': src_dep_params}
Expand Down Expand Up @@ -755,7 +773,7 @@ def calc_source_dependent_parameters(data, expected_src_pos_x_m, expected_src_po
return src_dep_params


def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m):
def get_expected_source_pos(data, data_type, config, effective_focal_length=29.30565 * u.m):
"""Get expected source position for source-dependent analysis .

Parameters
Expand All @@ -775,7 +793,7 @@ def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m):
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),
focal_length,
effective_focal_length,
u.Quantity(data['mc_alt_tel'].values, u.deg, copy=False),
u.Quantity(data['mc_az_tel'].values, u.deg, copy=False)
)
Expand Down Expand Up @@ -806,7 +824,7 @@ def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m):
obstime = Time(time, scale='utc', format='unix')
pointing_alt = u.Quantity(data['alt_tel'], u.rad, copy=False)
pointing_az = u.Quantity(data['az_tel'], u.rad, copy=False)
source_pos = utils.radec_to_camera(source_coord, obstime, pointing_alt, pointing_az, focal_length)
source_pos = utils.radec_to_camera(source_coord, obstime, pointing_alt, pointing_az, effective_focal_length)

expected_src_pos_x_m = source_pos.x.to_value(u.m)
expected_src_pos_y_m = source_pos.y.to_value(u.m)
Expand Down
17 changes: 12 additions & 5 deletions lstchain/scripts/lstchain_add_source_dependent_parameters.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@

import pandas as pd
from ctapipe.instrument import SubarrayDescription
from ctapipe_io_lst import OPTICS

from lstchain.io import (
get_standard_config,
Expand Down Expand Up @@ -61,11 +62,17 @@ def main():
pass

dl1_params = pd.read_hdf(dl1_filename, key=dl1_params_lstcam_key)
subarray_info = SubarrayDescription.from_hdf(dl1_filename)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length

src_dep_df = pd.concat(get_source_dependent_parameters(dl1_params, config, focal_length=focal_length), axis=1)
try:
subarray_info = SubarrayDescription.from_hdf(dl1_filename)
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.")
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"

effective_focal_length = OPTICS.effective_focal_length

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 •
edited
Loading

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

cta-lstchain/lstchain/reco/dl1_to_dl2.py

Line 801 in 500d7b2

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):

cta-lstchain/lstchain/reco/dl1_to_dl2.py

Lines 812 to 814 in 500d7b2

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)

cta-lstchain/lstchain/reco/dl1_to_dl2.py

Line 629 in 500d7b2

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.


metadata = global_metadata()
write_dataframe(src_dep_df, dl1_filename, dl1_params_src_dep_lstcam_key, config=config, meta=metadata)
Expand Down
22 changes: 14 additions & 8 deletions lstchain/scripts/lstchain_dl1_to_dl2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
import numpy as np
import pandas as pd
from ctapipe.instrument import SubarrayDescription
from ctapipe_io_lst import OPTICS
from tables import open_file

from lstchain.io import (
Expand Down Expand Up @@ -95,9 +96,14 @@ def apply_to_file(filename, models_dict, output_dir, config):
data.alt_tel = - np.pi / 2.
data.az_tel = - np.pi / 2.

subarray_info = SubarrayDescription.from_hdf(filename)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length
try:
subarray_info = SubarrayDescription.from_hdf(filename)
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.")
print("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

# Apply the models to the data

Expand All @@ -116,15 +122,15 @@ def apply_to_file(filename, models_dict, output_dir, config):
models_dict['cls_gh'],
models_dict['reg_energy'],
reg_disp_vector=models_dict['disp_vector'],
focal_length=focal_length,
effective_focal_length=effective_focal_length,
custom_config=config)
elif config['disp_method'] == 'disp_norm_sign':
dl2 = dl1_to_dl2.apply_models(data,
models_dict['cls_gh'],
models_dict['reg_energy'],
reg_disp_norm=models_dict['disp_norm'],
cls_disp_sign=models_dict['disp_sign'],
focal_length=focal_length,
effective_focal_length=effective_focal_length,
custom_config=config)

# Source-dependent analysis
Expand All @@ -136,7 +142,7 @@ def apply_to_file(filename, models_dict, output_dir, config):
# if not, source-dependent parameters are added now
else:
data_srcdep = pd.concat(dl1_to_dl2.get_source_dependent_parameters(
data, config, focal_length=focal_length), axis=1)
data, config, effective_focal_length=effective_focal_length), axis=1)

dl2_srcdep_dict = {}
srcindep_keys = data.keys()
Expand All @@ -157,15 +163,15 @@ def apply_to_file(filename, models_dict, output_dir, config):
models_dict['cls_gh'],
models_dict['reg_energy'],
reg_disp_vector=models_dict['disp_vector'],
focal_length=focal_length,
effective_focal_length=effective_focal_length,
custom_config=config)
elif config['disp_method'] == 'disp_norm_sign':
dl2_df = dl1_to_dl2.apply_models(data_with_srcdep_param,
models_dict['cls_gh'],
models_dict['reg_energy'],
reg_disp_norm=models_dict['disp_norm'],
cls_disp_sign=models_dict['disp_sign'],
focal_length=focal_length,
effective_focal_length=effective_focal_length,
custom_config=config)

dl2_srcdep = dl2_df.drop(srcindep_keys, axis=1)
Expand Down
16 changes: 11 additions & 5 deletions lstchain/scripts/lstchain_mc_rfperformance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
import matplotlib.pyplot as plt
import pandas as pd
from ctapipe.instrument import SubarrayDescription
from ctapipe_io_lst import OPTICS

from lstchain.io import (
read_configuration_file,
Expand Down Expand Up @@ -87,10 +88,15 @@ def main():

config = replace_config(standard_config, custom_config)

subarray_info = SubarrayDescription.from_hdf(args.gammatest)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length

try:
subarray_info = SubarrayDescription.from_hdf(args.gammatest)
tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1
effective_focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length
except OSError:
print("subarray table is not readable because of the version inompatibility.")
print("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

reg_energy, reg_disp_norm, cls_disp_sign, cls_gh = dl1_to_dl2.build_models(
args.gammafile,
args.protonfile,
Expand All @@ -110,7 +116,7 @@ def main():
data = pd.concat([gammas, proton], ignore_index=True)

dl2 = dl1_to_dl2.apply_models(data, cls_gh, reg_energy, reg_disp_norm=reg_disp_norm,
cls_disp_sign=cls_disp_sign, focal_length=focal_length,
cls_disp_sign=cls_disp_sign, effective_focal_length=effective_focal_length,
custom_config=config)

####PLOT SOME RESULTS#####
Expand Down
6 changes: 3 additions & 3 deletions lstchain/tests/test_lstchain.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,7 @@ def test_get_source_dependent_parameters_mc(simulated_dl1_file):
assert (src_dep_df_gamma['on']['expected_src_y'] == dl1_params['src_y']).all()

np.testing.assert_allclose(
src_dep_df_proton['on']['expected_src_x'], 0.195, atol=1e-2
src_dep_df_proton['on']['expected_src_x'], 0.205, atol=1e-2
)
np.testing.assert_allclose(
src_dep_df_proton['on']['expected_src_y'], 0., atol=1e-2
Expand Down Expand Up @@ -224,13 +224,13 @@ def test_get_source_dependent_parameters_observed(observed_dl1_files):
assert (src_dep_df_on['on']['expected_src_y'] == 0).all()

np.testing.assert_allclose(
src_dep_df_wobble['on']['expected_src_x'], -0.195, atol=1e-2
src_dep_df_wobble['on']['expected_src_x'], -0.205, atol=1e-2
)
np.testing.assert_allclose(
src_dep_df_wobble['on']['expected_src_y'], 0., atol=1e-2
)
np.testing.assert_allclose(
src_dep_df_wobble['off_180']['expected_src_x'], 0.195, atol=1e-2
src_dep_df_wobble['off_180']['expected_src_x'], 0.205, atol=1e-2
)
np.testing.assert_allclose(
src_dep_df_wobble['off_180']['expected_src_y'], 0., atol=1e-2
Expand Down