Merge pull request #1302 from cta-observatory/pointing_in_dl3

Computing the run pointing in sky coordinates

lstchain/high_level/hdu_table.py

@@ -14,6 +14,7 @@
 from astropy.coordinates import AltAz, SkyCoord
 from astropy.coordinates.erfa_astrom import ErfaAstromInterpolator, erfa_astrom
 from astropy.io import fits
+from astropy.stats import circmean
 from astropy.table import QTable, Table
 from astropy.time import Time
 
@@ -257,26 +258,49 @@ def get_timing_params(data, epoch=LST_EPOCH):
     return time_pars, time_utc
 
 
-def get_pointing_params(data, source_pos, time_utc):
+def get_pointing_params(data, source_pos, time_utc, exclude_fraction=0.2):
     """
-    Convert the telescope pointing position for the first event from AltAz
-    into ICRS frame of reference.
+    Convert the telescope pointing directions into ICRS frame of reference,
+    and average them for the run (excluding the first events, for which there
+    may be some mispointing due to not-yet-stable tracking).
+
+    exclude_fraction: fraction of the events that will be excluded from the 
+    averaging (the excluded events are the ones at the beginning of the run)
 
     Note: The angular difference from the source is just used for logging here.
+
+    Returns the average pointing in ICRS frame
     """
+
     pointing_alt = data["pointing_alt"]
     pointing_az = data["pointing_az"]
 
-    pnt_icrs = SkyCoord(
-        alt=pointing_alt[0],
-        az=pointing_az[0],
-        frame=AltAz(obstime=time_utc[0], location=LST1_LOCATION),
-    ).transform_to(frame="icrs")
+    max_off_angle = 0.05 * u.deg 
+    # Keep track of how often pointing is beyond max_off_angle w.r.t. the mean
+    # run pointing (computed excluding the beginning of run - see exclude_fraction)
+
+    first_event = int(exclude_fraction * len(data))
+    with erfa_astrom.set(ErfaAstromInterpolator(300 * u.s)):
+        evtwise_pnt_icrs = SkyCoord(
+            alt=pointing_alt,
+            az=pointing_az,
+            frame=AltAz(obstime=time_utc, location=LST1_LOCATION),
+        ).transform_to(frame="icrs")
+
+    mean_ra = circmean(evtwise_pnt_icrs[first_event:].ra)
+    mean_dec = np.mean(evtwise_pnt_icrs[first_event:].dec)
+    pnt_icrs = SkyCoord(ra=mean_ra, dec=mean_dec, frame="icrs")
+
+    off_fraction = ((pnt_icrs.separation(evtwise_pnt_icrs) > max_off_angle).sum() / 
+                    len(evtwise_pnt_icrs))
+
+    log.info(f"Mean pointing RA: {mean_ra.to(u.deg):.3f},  Dec: {mean_dec.to(u.deg):.3f}")
+    log.info(f"Fraction of events with pointing beyond {max_off_angle:.3f} of mean pointing: {off_fraction:.6f}")
 
     source_pointing_diff = source_pos.separation(pnt_icrs)
 
     log.info(
-        f"Source pointing difference with camera pointing is {source_pointing_diff:.3f}"
+        f"Angle between source direction and mean telescope pointing is {source_pointing_diff:.3f}"
     )
 
     return pnt_icrs