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

Computing the run pointing in sky coordinates #1302

Open
wants to merge 12 commits into
base: main
Choose a base branch
from
Open

Computing the run pointing in sky coordinates #1302

Changes from 9 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
befdc81
Computing the run pointing in sky coordinates
moralejo Oct 8, 2024
fa440af
Update conf.py (remove source of docs error?)
moralejo Oct 8, 2024
65f6e94
Check pointing stability
moralejo Oct 8, 2024
82ffcbc
remove f in string with no vars
moralejo Oct 8, 2024
43080af
Merge branch 'main' of github.com:cta-observatory/cta-lstchain into p…
moralejo Oct 9, 2024
c90f8c6
Merge branch 'main' of github.com:cta-observatory/cta-lstchain into p…
moralejo Oct 9, 2024
4609a58
Fix comment
moralejo Oct 10, 2024
8406eaf
warn => warning
moralejo Oct 10, 2024
6fc28d3
Update lstchain/high_level/hdu_table.py
moralejo Oct 10, 2024
b81354f
Changed log message related to unstable pointing
moralejo Oct 14, 2024
14399c2
Corrected comment
moralejo Oct 18, 2024
9df60c6
Fixed comment
moralejo Oct 18, 2024
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
45 changes: 36 additions & 9 deletions lstchain/high_level/hdu_table.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -257,26 +258,52 @@
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):
moralejo marked this conversation as resolved.
Show resolved Hide resolved
"""
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_median_angle = 0.05 * u.deg # Is this a reasonable limit?
# If median of the angle between the pointing directions and their mean in the run
moralejo marked this conversation as resolved.
Show resolved Hide resolved
# is larger than max_median_angle a warning will be displayed about unstable pointing.

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

median_angle_wrt_mean_pointing = np.median(pnt_icrs.separation(evtwise_pnt_icrs))
median_angle_wrt_mean_pointing_2 = np.median(pnt_icrs.separation(evtwise_pnt_icrs[first_event:]))

log.info("Median of angle between pointing direction and mean pointing:")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This compares the median over all events to the median over the events without the ones excluded by "fraction".

Is this really what we want? Normal runs are long (~20 min) and the unstable pointing we saw was for the first ~35 seconds in the bad cases.

So the median won't be affected by the the first 35 seconds of unstable pointing and then stable pointing.

We looked at the deviation of the individual pointings from the "good pointing" in chunks of 2 seconds.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we want to detect unstable outlier values at the beginning, why don't we just use a sigma clipping routine from scipy or astropy?

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 compares the median over all events to the median over the events without the ones excluded by "fraction".

Is this really what we want? Normal runs are long (~20 min) and the unstable pointing we saw was for the first ~35 seconds in the bad cases.

So the median won't be affected by the the first 35 seconds of unstable pointing and then stable pointing.

We looked at the deviation of the individual pointings from the "good pointing" in chunks of 2 seconds.

Right, this will not be sensitive to the instabilities if they are so small... Perhaps it makes more sense to simply write out in the log what fraction of the run has a pointing beyond, say, 0.05 deg of the target. But I do not know where it would make sense to set a threshold for a warning.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we want to detect unstable outlier values at the beginning, why don't we just use a sigma clipping routine from scipy or astropy?

For simplicity, we know the problem occurs only at the beginning of the run. At least the problem we are addressing in this PR. And sigma clipping requires some settings, not obvious to me which ones would be adequate here.

log.info(f" All events: {median_angle_wrt_mean_pointing:.3f}")
if median_angle_wrt_mean_pointing > max_median_angle:
log.warning(" ==> The telescope pointing seems to be unstable during the Run!")

Check warning on line 300 in lstchain/high_level/hdu_table.py

View check run for this annotation

Codecov / codecov/patch

lstchain/high_level/hdu_table.py#L300 

Added line #L300 was not covered by tests
log.info(f" Excluding the first {exclude_fraction:.1%} of events: {median_angle_wrt_mean_pointing_2:.3f}")

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
Expand Down