Merge pull request #1132 from cta-observatory/wrap_az_tel

Normalize `az_tel` to the range of [0, 360) and compute sin_az_tel in dl1_to_dl2 step

lstchain/reco/dl1_to_dl2.py

@@ -13,7 +13,7 @@
 import joblib
 import numpy as np
 import pandas as pd
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, Angle
 from astropy.time import Time
 from pathlib import Path
 from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
@@ -392,6 +392,16 @@ def build_models(filegammas, fileprotons,
                                                   + config['disp_classification_features'],
                                     )
 
+    # Normalize all azimuth angles to the range [0, 360) degrees
+    df_gamma.az_tel = Angle(df_gamma.az_tel, u.rad).wrap_at(360 * u.deg).rad
+    df_proton.az_tel = Angle(df_proton.az_tel, u.rad).wrap_at(360 * u.deg).rad
+
+    # Dealing with `sin_az_tel` missing data because of the former version of lstchain
+    if 'sin_az_tel' not in df_gamma.columns:
+        df_gamma['sin_az_tel'] = np.sin(df_gamma.az_tel)
+    if 'sin_az_tel' not in df_proton.columns:
+        df_proton['sin_az_tel'] = np.sin(df_proton.az_tel)
+
     # Training MC gammas in reduced viewcone
     src_r_min = config['train_gamma_src_r_deg'][0]
     src_r_max = config['train_gamma_src_r_deg'][1]

lstchain/scripts/lstchain_dl1_to_dl2.py

@@ -12,6 +12,8 @@
 
 import numpy as np
 import pandas as pd
+import astropy.units as u
+from astropy.coordinates import Angle
 from ctapipe.instrument import SubarrayDescription
 from tables import open_file
 
@@ -95,6 +97,13 @@ def apply_to_file(filename, models_dict, output_dir, config):
             data.alt_tel = - np.pi / 2.
             data.az_tel = - np.pi / 2.
 
+    # Normalize all azimuth angles to the range [0, 360) degrees 
+    data.az_tel = Angle(data.az_tel, u.rad).wrap_at(360 * u.deg).rad
+
+    # Dealing with `sin_az_tel` missing data because of the former version of lstchain
+    if 'sin_az_tel' not in data.columns:
+        data['sin_az_tel'] = np.sin(data.az_tel)
+
     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