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

Changes from 3 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
78 changes: 65 additions & 13 deletions lstchain/reco/dl1_to_dl2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@
'train_energy',
'train_reco',
'train_sep',
'update_disp_with_effective_focal_length'
]


Expand Down Expand Up @@ -109,16 +110,16 @@
features = config['disp_regression_features']
model = RandomForestRegressor

logger.info("Given features: ", features)
logger.info("Number of events for training: ", train.shape[0])
logger.info("Training model {} for disp vector regression".format(model))

Check warning on line 115 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L113-L115 

Added lines #L113 - L115 were not covered by tests

reg = model(**disp_regression_args)
x = train[features]
y = np.transpose([train[f] for f in predict_features])
reg.fit(x, y)

logger.info("Model {} trained!".format(model))

Check warning on line 122 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L122 

Added line #L122 was not covered by tests

return reg

Expand Down Expand Up @@ -215,24 +216,24 @@
disp_features = config['disp_regression_features']
model = RandomForestRegressor

logger.info("Given energy_features: ", energy_features)
logger.info("Number of events for training: ", train.shape[0])
logger.info("Training Random Forest Regressor for Energy Reconstruction...")

Check warning on line 221 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L219-L221 

Added lines #L219 - L221 were not covered by tests

reg_energy = model(**energy_regression_args)
reg_energy.fit(train[energy_features],
train['log_mc_energy'])

logger.info("Random Forest trained!")
logger.info("Given disp_features: ", disp_features)
logger.info("Training Random Forest Regressor for disp_norm Reconstruction...")

Check warning on line 229 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L227-L229 

Added lines #L227 - L229 were not covered by tests

reg_disp = RandomForestRegressor(**disp_regression_args)
reg_disp.fit(train[disp_features],
train['disp_norm'])

logger.info("Random Forest trained!")
logger.info("Done!")

Check warning on line 236 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L235-L236 

Added lines #L235 - L236 were not covered by tests
return reg_energy, reg_disp


Expand Down Expand Up @@ -346,22 +347,26 @@
df_gamma = pd.read_hdf(filegammas, key=dl1_params_lstcam_key)
df_proton = pd.read_hdf(fileprotons, key=dl1_params_lstcam_key)

# Update parameters related to target direction on camera frame for gamma MC
# taking into account of the abrration effect using effective 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:
logger.warning("subarray table is not readable because of the version inompatibility.")
logger.warning("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

Check warning on line 359 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L356-L359 

Added lines #L356 - L359 were not covered by tests

df_gamma = update_disp_with_effective_focal_length(df_gamma, effective_focal_length = effective_focal_length)

if config['source_dependent']:
# if source-dependent parameters are already in dl1 data, just read those data
# if not, source-dependent parameters are added here
if dl1_params_src_dep_lstcam_key in get_dataset_keys(filegammas):
src_dep_df_gamma = get_srcdep_params(filegammas)

else:
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:
logger.warning("subarray table is not readable because of the version inompatibility.")
logger.warning("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
)
Expand All @@ -378,10 +383,10 @@
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:
logger.warning("subarray table is not readable because of the version inompatibility.")
logger.warning("Use the effective focal lentgh for the standard LST optics")
effective_focal_length = OPTICS.effective_focal_length

Check warning on line 389 in lstchain/reco/dl1_to_dl2.py

View check run for this annotation

Codecov / codecov/patch

lstchain/reco/dl1_to_dl2.py#L386-L389 

Added lines #L386 - L389 were not covered by tests

src_dep_df_proton = get_source_dependent_parameters(
df_proton, config, effective_focal_length=effective_focal_length
Expand Down Expand Up @@ -607,6 +612,13 @@
+ config['disp_classification_features'],
)

# Update parameters related to target direction on camera frame for MC data
# taking into account of the abrration effect using effective focal length
is_simu = 'disp_norm' in dl2.columns
if is_simu:
dl2 = update_disp_with_effective_focal_length(dl2, effective_focal_length = effective_focal_length)


# Reconstruction of Energy and disp_norm distance
if isinstance(reg_energy, (str, bytes, Path)):
reg_energy = joblib.load(reg_energy)
Expand Down Expand Up @@ -794,11 +806,11 @@
# 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.

cta-lstchain/lstchain/data/lstchain_src_dep_config.json

Lines 62 to 63 in c9e4c4f

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

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

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.

u.Quantity(data['mc_az_tel'].values + np.deg2rad(config['mc_nominal_source_y_deg']), u.rad, copy=False),
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)
u.Quantity(data['mc_alt_tel'].values, u.rad, copy=False),
u.Quantity(data['mc_az_tel'].values, u.rad, copy=False)
)

expected_src_pos_x_m = expected_src_pos.x.to_value(u.m)
Expand Down Expand Up @@ -838,3 +850,43 @@
)

return expected_src_pos_x_m, expected_src_pos_y_m


def update_disp_with_effective_focal_length(data, effective_focal_length=29.30565 * u.m):
"""Update disp parameters using effective focal length

Parameters
----------
data: Pandas DataFrame
config: dictionnary containing configuration
"""

source_pos_in_camera = utils.sky_to_camera(
u.Quantity(data['mc_alt'].values, u.rad, copy=False),
u.Quantity(data['mc_az'].values, u.rad, copy=False),
effective_focal_length,
u.Quantity(data['mc_alt_tel'].values, u.rad, copy=False),
u.Quantity(data['mc_az_tel'].values, u.rad, copy=False)
)

expected_src_pos_x_m = source_pos_in_camera.x.to_value(u.m)
expected_src_pos_y_m = source_pos_in_camera.y.to_value(u.m)

data['src_x'] = expected_src_pos_x_m
data['src_y'] = expected_src_pos_y_m

disp_dx, disp_dy, disp_norm, disp_angle, disp_sign = disp.disp(
data['x'].values,
data['y'].values,
expected_src_pos_x_m,
expected_src_pos_y_m,
data['psi'].values
)

data['disp_dx'] = disp_dx
data['disp_dy'] = disp_dy
data['disp_norm'] = disp_norm
data['disp_angle'] = disp_angle
data['disp_sign'] = disp_sign

return data
Loading