Coverage for src / lstautorta / crab_cleanup.py: 0%
167 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-08-10 11:56 +0000
1#!/usr/bin/env python3
3"""
4Scan DL1 HDF5 files, estimate telescope pointing vs Crab in AltAz, and decide KEEP/DELETE per obs_id.
5Also collects associated logs and exports deletion/keep lists + CSV summary.
7Log file is timestamped to avoid overwriting between runs.
9Example : "python3 crab_cleanup.py --base-day-dir /fefs/onsite/pipeline/rta/data/2025/12/ --out-dir ./crab_cleanup_out --max-sep-deg 5.0 --tel-group tel_001 --step 200"
10"""
12import argparse
13import logging
14import re
15from datetime import UTC, datetime
16from pathlib import Path
18import h5py
19import numpy as np
20import pandas as pd
21from astropy import units as u
22from astropy.coordinates import AltAz, EarthLocation, SkyCoord
23from astropy.time import Time
25# -------------------------
26# Patterns / helpers
27# -------------------------
29DL1_RE = re.compile(r"^(dl1_|dl1_v).+\.h5$")
30OBS_RE = re.compile(r"obs_id_(\d+)")
33def is_dl1_file(p: Path) -> bool:
34 return p.is_file() and DL1_RE.match(p.name) is not None
37def extract_obs_id_from_name(filename: str):
38 m = OBS_RE.search(filename)
39 return int(m.group(1)) if m else None
42def scan_dl1_files(day_dir: Path):
43 """
44 Recursively scan all DL1 *.h5 under day_dir.
45 Returns dict: obs_id -> list[Path]
46 """
47 by_obs = {}
48 for p in day_dir.rglob("*.h5"):
49 if not is_dl1_file(p):
50 continue
51 obs = extract_obs_id_from_name(p.name)
52 if obs is None:
53 continue
54 by_obs.setdefault(obs, []).append(p)
56 for obs in by_obs:
57 by_obs[obs] = sorted(by_obs[obs])
58 return by_obs
61# -------------------------
62# DL1 reading + separation
63# -------------------------
66def open_dl1_get_pointing_and_time(dl1_path: Path, tel_group: str, step: int):
67 """
68 Extract a subsample of (time, alt_tel, az_tel) from DL1.
70 Assumptions:
71 - dl1/event/subarray/trigger/time = unix seconds (UTC)
72 - dl1/event/telescope/parameters/<tel_group>/alt_tel, az_tel = radians
73 """
74 with h5py.File(dl1_path, "r") as f:
75 trig = f["dl1/event/subarray/trigger"]
76 pars = f[f"dl1/event/telescope/parameters/{tel_group}"]
78 n = len(trig)
79 if n == 0:
80 return np.array([]), np.array([]), np.array([])
82 idx = np.arange(0, n, step, dtype=int)
84 times_unix = trig["time"][idx].astype(float)
85 alt_rad = pars["alt_tel"][idx].astype(float)
86 az_rad = pars["az_tel"][idx].astype(float)
88 return times_unix, alt_rad, az_rad
91def crab_separation_stats(
92 dl1_path: Path, site: EarthLocation, crab: SkyCoord, tel_group: str, step: int, max_sep_deg: float
93):
94 """
95 Compute separation stats (deg) between telescope pointing and Crab in AltAz.
96 Default decision: median separation < max_sep_deg -> KEEP else DELETE.
97 """
98 times_unix, alt_tel, az_tel = open_dl1_get_pointing_and_time(dl1_path, tel_group=tel_group, step=step)
100 if len(times_unix) == 0:
101 return {
102 "sep_min_deg": np.nan,
103 "sep_med_deg": np.nan,
104 "sep_mean_deg": np.nan,
105 "sep_max_deg": np.nan,
106 "decision": "ERROR",
107 "error": "no events",
108 }
110 t = Time(times_unix, format="unix", scale="utc")
112 crab_altaz = crab.transform_to(AltAz(obstime=t, location=site))
113 tel_altaz = SkyCoord(az=az_tel * u.rad, alt=alt_tel * u.rad, frame=AltAz(obstime=t, location=site))
115 sep_deg = tel_altaz.separation(crab_altaz).to_value(u.deg)
117 stats = {
118 "sep_min_deg": float(np.min(sep_deg)),
119 "sep_med_deg": float(np.median(sep_deg)),
120 "sep_mean_deg": float(np.mean(sep_deg)),
121 "sep_max_deg": float(np.max(sep_deg)),
122 "error": "",
123 }
125 keep = stats["sep_med_deg"] < max_sep_deg
126 stats["decision"] = "KEEP" if keep else "DELETE"
127 return stats
130# -------------------------
131# Logs discovery
132# -------------------------
135def find_obs_root_from_dl1(dl1_path: Path, obs_id: int) -> Path | None:
136 """
137 Try to locate the directory path that corresponds to ".../<obs_id>/" within dl1_path.
138 """
139 parts = dl1_path.parts
140 obs_str = str(obs_id)
141 if obs_str not in parts:
142 return None
143 i = parts.index(obs_str)
144 return Path(*parts[: i + 1])
147def find_logs_for_dl1(dl1_path: Path, log_subdir: Path):
148 """
149 Look for logs in <obs_root>/<log_subdir>/ matching:
150 <dl1_prefix>.o-* and <dl1_prefix>.e-*
151 """
152 obs = extract_obs_id_from_name(dl1_path.name)
153 if obs is None:
154 return []
156 obs_root = find_obs_root_from_dl1(dl1_path, obs)
157 if obs_root is None:
158 return []
160 log_dir = obs_root / log_subdir
161 if not log_dir.exists():
162 return []
164 prefix = dl1_path.stem
165 logs = []
166 logs.extend(sorted(log_dir.glob(f"{prefix}.o-*")))
167 logs.extend(sorted(log_dir.glob(f"{prefix}.e-*")))
168 return [p for p in logs]
171# -------------------------
172# Main
173# -------------------------
176def setup_logger(out_dir: Path, prefix: str = "crab_pointing_cleanup"):
177 out_dir.mkdir(parents=True, exist_ok=True)
179 # Timestamp (UTC) in filename to avoid overwriting
180 ts = datetime.now(UTC).strftime("%Y%m%d_%H%M%S_UTC")
181 log_path = out_dir / f"{prefix}_{ts}.log"
183 logger = logging.getLogger("crab_cleanup")
184 logger.setLevel(logging.INFO)
185 logger.handlers.clear()
187 fmt = logging.Formatter("%(asctime)s | %(levelname)s | %(message)s")
189 fh = logging.FileHandler(log_path)
190 fh.setLevel(logging.INFO)
191 fh.setFormatter(fmt)
193 sh = logging.StreamHandler()
194 sh.setLevel(logging.INFO)
195 sh.setFormatter(fmt)
197 logger.addHandler(fh)
198 logger.addHandler(sh)
200 logger.info("Log file: %s", log_path)
201 logger.info("Script start UTC: %s", datetime.now(UTC).isoformat())
202 return logger, log_path
205def main():
206 parser = argparse.ArgumentParser(description="KEEP/DELETE DL1 based on telescope pointing separation to Crab.")
207 parser.add_argument(
208 "--base-day-dir",
209 type=Path,
210 required=True,
211 help="Day directory to scan (e.g. /.../2026/01/20 or /.../2026/01/)",
212 )
213 parser.add_argument(
214 "--max-sep-deg", type=float, default=5.0, help="Max median separation (deg) to KEEP (default: 5.0)"
215 )
216 parser.add_argument(
217 "--tel-group", type=str, default="tel_001", help="HDF5 telescope group name (default: tel_001)"
218 )
219 parser.add_argument("--step", type=int, default=200, help="Subsampling step (default: 200)")
220 parser.add_argument(
221 "--log-subdir",
222 type=Path,
223 default=Path("logs/dl1_alt_az"),
224 help="Logs path relative to obs_id directory (default: logs/dl1_alt_az)",
225 )
226 parser.add_argument(
227 "--out-dir",
228 type=Path,
229 default=Path(),
230 help="Output directory for txt/csv/log files (default: current directory)",
231 )
232 parser.add_argument("--orm-height-m", type=float, default=2200.0, help="Site height in meters (default: 2200)")
233 args = parser.parse_args()
235 logger, _ = setup_logger(args.out_dir)
237 base_day_dir = args.base_day_dir
238 if not base_day_dir.exists():
239 raise FileNotFoundError(f"base-day-dir does not exist: {base_day_dir}")
241 # ORM / La Palma (LST-1)
242 site = EarthLocation(lat=28.7617 * u.deg, lon=-17.89 * u.deg, height=args.orm_height_m * u.m)
244 # Crab (ICRS J2000)
245 crab = SkyCoord(ra=83.6331 * u.deg, dec=22.0145 * u.deg, frame="icrs")
247 logger.info("Scanning DL1 under: %s", base_day_dir)
248 dl1_by_obs = scan_dl1_files(base_day_dir)
249 logger.info("Unique obs_id: %d", len(dl1_by_obs))
250 logger.info("Example obs_id: %s", list(dl1_by_obs.keys())[:5])
252 rows = []
253 for obs_id, files in sorted(dl1_by_obs.items()):
254 rep = files[0]
255 try:
256 stats = crab_separation_stats(
257 rep, site=site, crab=crab, tel_group=args.tel_group, step=args.step, max_sep_deg=args.max_sep_deg
258 )
259 row = {"obs_id": obs_id, "rep_file": str(rep), **stats}
260 rows.append(row)
261 logger.info(
262 "obs_id=%s rep=%s decision=%s sep_med=%.3f deg",
263 obs_id,
264 rep.name,
265 row["decision"],
266 row["sep_med_deg"] if np.isfinite(row["sep_med_deg"]) else float("nan"),
267 )
268 except Exception as e:
269 rows.append(
270 {
271 "obs_id": obs_id,
272 "rep_file": str(rep),
273 "sep_min_deg": np.nan,
274 "sep_med_deg": np.nan,
275 "sep_mean_deg": np.nan,
276 "sep_max_deg": np.nan,
277 "decision": "ERROR",
278 "error": repr(e),
279 }
280 )
281 logger.exception("obs_id=%s rep=%s ERROR: %r", obs_id, rep.name, e)
283 df = pd.DataFrame(rows).sort_values(["decision", "obs_id"])
284 logger.info("Decision counts:\n%s", df["decision"].value_counts(dropna=False).to_string())
286 decision_by_obs = {int(r.obs_id): r.decision for _, r in df.iterrows()}
288 to_delete_dl1 = []
289 to_delete_logs = []
290 to_keep_dl1 = []
292 for obs_id, files in dl1_by_obs.items():
293 dec = decision_by_obs.get(obs_id, "ERROR")
294 if dec == "KEEP":
295 to_keep_dl1.extend(files)
296 elif dec == "DELETE":
297 to_delete_dl1.extend(files)
298 for f in files:
299 to_delete_logs.extend(find_logs_for_dl1(f, log_subdir=args.log_subdir))
301 to_delete_dl1 = sorted(set(to_delete_dl1))
302 to_delete_logs = sorted(set(to_delete_logs))
303 to_keep_dl1 = sorted(set(to_keep_dl1))
305 logger.info("DL1 to delete: %d", len(to_delete_dl1))
306 logger.info("Logs to delete: %d", len(to_delete_logs))
307 logger.info("DL1 to keep: %d", len(to_keep_dl1))
309 # Outputs
310 out_dir = args.out_dir
311 out_del_dl1 = out_dir / "to_delete_dl1.txt"
312 out_del_logs = out_dir / "to_delete_logs.txt"
313 out_keep = out_dir / "to_keep_dl1.txt"
314 out_csv = out_dir / "obsid_crab_pointing_summary.csv"
316 out_del_dl1.write_text("\n".join(str(p) for p in to_delete_dl1) + ("\n" if to_delete_dl1 else ""))
317 out_del_logs.write_text("\n".join(str(p) for p in to_delete_logs) + ("\n" if to_delete_logs else ""))
318 out_keep.write_text("\n".join(str(p) for p in to_keep_dl1) + ("\n" if to_keep_dl1 else ""))
319 df.to_csv(out_csv, index=False)
321 logger.info("Wrote: %s", out_del_dl1)
322 logger.info("Wrote: %s", out_del_logs)
323 logger.info("Wrote: %s", out_keep)
324 logger.info("Wrote: %s", out_csv)
326 keep_obs = df.loc[df["decision"] == "KEEP", "obs_id"].to_list()
327 logger.info("KEEP obs_id list: %s", keep_obs)
329 # Shell helpers
330 logger.info("### Dry-run command:")
331 logger.info('cat "%s" "%s" | while read f; do echo rm -f "$f"; done', out_del_dl1, out_del_logs)
333 logger.info("### Delete command:")
334 logger.info('cat "%s" "%s" | while read f; do rm -f "$f"; done', out_del_dl1, out_del_logs)
337if __name__ == "__main__":
338 main()