Merge pull request #1027 from cta-observatory/disp_sign_fix

Fix definition of "true disp_sign" parameter

lstchain/reco/disp.py

@@ -14,7 +14,7 @@
     ]
 
 
-def disp(cog_x, cog_y, src_x, src_y):
+def disp(cog_x, cog_y, src_x, src_y, hillas_psi):
     """
     Compute the disp parameters
 
@@ -24,6 +24,7 @@ def disp(cog_x, cog_y, src_x, src_y):
     cog_y: `numpy.ndarray` or float
     src_x: `numpy.ndarray` or float
     src_y: `numpy.ndarray` or float
+    hillas_psi: `numpy.ndarray` or float
 
     Returns
     -------
@@ -36,7 +37,16 @@ def disp(cog_x, cog_y, src_x, src_y):
     """
     disp_dx = src_x - cog_x
     disp_dy = src_y - cog_y
-    disp_norm = np.sqrt(disp_dx**2 + disp_dy**2)
+
+    disp_norm = disp_dx * np.cos(hillas_psi) + disp_dy * np.sin(hillas_psi)
+    disp_sign = np.sign(disp_norm)
+    disp_norm = np.abs(disp_norm)
+
+    # disp_sign : indicates in which direction, "positive" or "negative", we must move along the 
+    # reconstructed image axis (with direction defined by the versor cos(hillas_psi), sin(hillas_psi)) 
+    # we must move from cog_x, cog_y to get closest to the true direction (src_x, src_y)
+
+
     if hasattr(disp_dx, '__len__'):
         disp_angle = np.arctan(disp_dy / disp_dx)
         disp_angle[disp_dx == 0] = np.pi / 2. * np.sign(disp_dy[disp_dx == 0])
@@ -46,8 +56,6 @@ def disp(cog_x, cog_y, src_x, src_y):
         else:
             disp_angle = np.arctan(disp_dy/disp_dx)
 
-    disp_sign = np.sign(disp_dx)
-
     return disp_dx, disp_dy, disp_norm, disp_angle, disp_sign
 
 
@@ -68,14 +76,15 @@ def miss(disp_dx, disp_dy, hillas_psi):
     return np.abs(np.sin(hillas_psi) * disp_dx - np.cos(hillas_psi)*disp_dy)
 
 
-def disp_parameters(cog_x, cog_y, mc_alt, mc_az, mc_alt_tel, mc_az_tel, focal):
+def disp_parameters(cog_x, cog_y, hillas_psi, mc_alt, mc_az, mc_alt_tel, mc_az_tel, focal):
     """
     Compute disp parameters.
 
     Parameters
     ----------
     cog_x: `numpy.ndarray` or float
     cog_y: `numpy.ndarray` or float
+    hillas_psi: `numpy.ndarray` or float
     mc_alt: `numpy.ndarray` or float
     mc_az: `numpy.ndarray` or float
     mc_alt_tel: `numpy.ndarray` or float
@@ -87,7 +96,7 @@ def disp_parameters(cog_x, cog_y, mc_alt, mc_az, mc_alt_tel, mc_az_tel, focal):
     (disp_dx, disp_dy, disp_norm, disp_angle, disp_sign) : `numpy.ndarray` or float
     """
     source_pos_in_camera = utils.sky_to_camera(mc_alt, mc_az, focal, mc_alt_tel, mc_az_tel)
-    return disp(cog_x, cog_y, source_pos_in_camera.x, source_pos_in_camera.y)
+    return disp(cog_x, cog_y, source_pos_in_camera.x, source_pos_in_camera.y, hillas_psi)
 
 
 
@@ -110,20 +119,21 @@ def disp_parameters_event(hillas_parameters, source_pos_x, source_pos_y):
     """
     disp_container = lstcontainers.DispContainer()
 
-    d = disp(hillas_parameters.x.to(u.m).value,
-             hillas_parameters.y.to(u.m).value,
-             source_pos_x.to(u.m).value,
-             source_pos_y.to(u.m).value,
+    d = disp(hillas_parameters.x.to_value(u.m),
+             hillas_parameters.y.to_value(u.m),
+             source_pos_x.to_value(u.m),
+             source_pos_y.to_value(u.m),
+             hillas_parameters.psi.to_value(u.rad)
              )
 
     disp_container.dx = d[0] * u.m
     disp_container.dy = d[1] * u.m
     disp_container.norm = d[2] * u.m
     disp_container.angle = d[3] * u.rad
     disp_container.sign = d[4]
-    disp_container.miss = miss(disp_container.dx.value,
-                               disp_container.dy.value,
-                               hillas_parameters.psi.to(u.rad).value) * u.m
+    disp_container.miss = miss(disp_container.dx.to_value(u.m),
+                               disp_container.dy.to_value(u.m),
+                               hillas_parameters.psi.to_value(u.rad))
     return disp_container
 
 

lstchain/reco/r0_to_dl1.py

@@ -771,7 +771,8 @@ def add_disp_to_parameters_table(dl1_file, table_path, focal):
     disp_parameters = disp.disp(df.x.values * u.m,
                                 df.y.values * u.m,
                                 source_pos_in_camera.x,
-                                source_pos_in_camera.y)
+                                source_pos_in_camera.y,
+                                df.psi.values * u.rad)
 
     with tables.open_file(dl1_file, mode="a") as file:
         tab = file.root[table_path]

lstchain/scripts/lstchain_dl1ab.py

@@ -317,7 +317,8 @@ def main():
                         dl1_container['x'].to_value(u.m),
                         dl1_container['y'].to_value(u.m),
                         params['src_x'][ii],
-                        params['src_y'][ii]
+                        params['src_y'][ii],
+                        dl1_container['psi'].to_value(u.rad)
                     )
 
                     dl1_container['disp_dx'] = disp_dx