#!/usr/bin/env python3 """No-fit xsherpa comparison of M104 R0-15 MOS1, MOS2, and pn0 spectra.""" from __future__ import annotations import csv import json from pathlib import Path from typing import Any import numpy as np from astropy.io import fits ROOT = Path(__file__).resolve().parent OUT = ROOT / ( "joint_spectrum_fitting_2T_basedon_region_v22/" "r015_epic_xsherpa_inspect_pn0_pattern0_20260713" ) ENERGY_RANGE = (0.4, 3.2) FEL_ZOOM = (0.6, 1.2) SB_FLUX_YLIM = (9.0e-8, 3.0e-5) SAMPLING_RATE = 5 SNR = 6 CENTER_RA_DEG = 189.99763333333333 CENTER_DEC_DEG = -11.623055555555556 _INPUTS = ( { "label": "MOS1", "nominal_region": "R0-15", "region_qualification": "R0-15 compatibility symlink to masked camera full-FOV extraction", "pha": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005-obj.pi", "bkg": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005-back.pi", "arf": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005.arf", "rmf": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005.rmf", "image": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005-obj-im-500-2000.fits", "mask": ROOT / "data/0900170101/rpc/4.background/bkg_mos1_00500-02000_exclude_extent_source_regionR0-15_mask/mos1U005-mask-im.fits", }, { "label": "MOS2", "nominal_region": "R0-15", "region_qualification": "R0-15 compatibility symlink to masked camera full-FOV extraction", "pha": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005-obj.pi", "bkg": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005-back.pi", "arf": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005.arf", "rmf": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005.rmf", "image": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005-obj-im-500-2000.fits", "mask": ROOT / "data/0900170101/rpc/4.background/bkg_mos2_00500-02000_exclude_extent_source_regionR0-15_mask/mos2U005-mask-im.fits", }, { "label": "PN0 (PATTERN=0)", "nominal_region": "R0-15", "region_qualification": "effective R0-15 within PN camera footprint", "pha": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002-obj.pi", "bkg": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002-back.pi", "arf": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002.arf", "rmf": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002.rmf", "image": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002-obj-im-500-2000.fits", "mask": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002-mask-im.fits", }, ) BANDS = ( (0.4, 3.2), (0.6, 1.2), (0.7, 1.05), (0.7, 0.875), (0.875, 1.05), (1.15, 1.45), (2.0, 3.2), ) def _hdu_header(path: Path) -> fits.Header: with fits.open(path, memmap=False) as hdul: return hdul["SPECTRUM"].header.copy() def validate_inputs() -> list[dict[str, Any]]: """Validate the explicit PHA/QPB/ARF/RMF manifest and return metadata.""" rows: list[dict[str, Any]] = [] for spec in _INPUTS: for key in ("pha", "bkg", "arf", "rmf", "image", "mask"): path = Path(spec[key]) if not path.is_file() or path.stat().st_size == 0: raise FileNotFoundError(f"Missing or empty {key}: {path}") header = _hdu_header(Path(spec["pha"])) image_header = fits.getheader(spec["image"]) mask_header = fits.getheader(spec["mask"]) expr = str(header.get("SLCTEXPR", "")) instrument = str(header.get("INSTRUME", "")) if spec["label"].startswith("PN0"): compact = expr.replace(" ", "").upper() if "PATTERN<=0" not in compact or "FLAG==0" not in compact: raise ValueError(f"PN0 selection is not PATTERN<=0 and FLAG==0: {expr}") elif "PATTERN<=12" not in expr.replace(" ", "").upper(): raise ValueError(f"MOS selection is not PATTERN<=12: {expr}") path_keys = {"pha", "bkg", "arf", "rmf", "image", "mask"} region_dir = Path(spec["pha"]).parent rows.append( { **{key: str(Path(value).absolute()) if key in path_keys else value for key, value in spec.items()}, "resolved_paths": { key: str(Path(spec[key]).resolve()) for key in path_keys }, "region_directory_is_symlink": region_dir.is_symlink(), "resolved_region_directory": str(region_dir.resolve()), "instrument": instrument, "exposure_s": float(header["EXPOSURE"]), "backscal": float(header["BACKSCAL"]), "skyarea_arcmin2": float(header["BACKSCAL"]) * (1.0 / 20.0 / 60.0) ** 2, "detector_channels": int(header["DETCHANS"]), "image_exposure_s": float(image_header.get("EXPOSURE", float("nan"))), "mask_has_celestial_wcs": "CTYPE1" in mask_header and "CTYPE2" in mask_header, "selection_expression": expr, "selection_expression_status": ( "Pattern/FLAG provenance only; the long spatial SLCTEXPR contains malformed/truncated continuation text" ), } ) instruments = [row["instrument"] for row in rows] if instruments != ["EMOS1", "EMOS2", "EPN"]: raise ValueError(f"Unexpected detector order: {instruments}") return rows def _pha_arrays(pha_path: str, bkg_path: str, rmf_path: str) -> dict[str, np.ndarray | float]: with fits.open(pha_path, memmap=False) as hdul: src_h = hdul["SPECTRUM"].header src = np.asarray(hdul["SPECTRUM"].data["COUNTS"], dtype=float) src_exp = float(src_h["EXPOSURE"]) src_backscal = float(src_h["BACKSCAL"]) src_areascal = float(src_h.get("AREASCAL", 1.0)) with fits.open(bkg_path, memmap=False) as hdul: bkg_h = hdul["SPECTRUM"].header bkg_exp = float(bkg_h["EXPOSURE"]) bkg_backscal = float(bkg_h["BACKSCAL"]) bkg_areascal = float(bkg_h.get("AREASCAL", 1.0)) bkg_data = hdul["SPECTRUM"].data # type: ignore[attr-defined] if "COUNTS" in bkg_data.names: bkg = np.asarray(bkg_data["COUNTS"], dtype=float) error_scale = 1.0 elif "RATE" in bkg_data.names: bkg = np.asarray(bkg_data["RATE"], dtype=float) * bkg_exp error_scale = bkg_exp else: raise ValueError(f"QPB PHA has neither COUNTS nor RATE: {bkg_path}") if "STAT_ERR" in hdul["SPECTRUM"].data.names: bkg_err = np.asarray(bkg_data["STAT_ERR"], dtype=float) * error_scale else: bkg_err = np.sqrt(np.clip(bkg, 0.0, None)) with fits.open(rmf_path, memmap=False) as hdul: ebounds = hdul["EBOUNDS"].data e_min = np.asarray(ebounds["E_MIN"], dtype=float) e_max = np.asarray(ebounds["E_MAX"], dtype=float) if not (len(src) == len(bkg) == len(e_min) == len(e_max)): raise ValueError( f"Channel mismatch: src={len(src)}, bkg={len(bkg)}, rmf={len(e_min)}" ) alpha = ( src_exp * src_backscal * src_areascal / (bkg_exp * bkg_backscal * bkg_areascal) ) return { "src": src, "bkg": bkg, "bkg_err": bkg_err, "e_min": e_min, "e_max": e_max, "alpha": float(alpha), "exposure": src_exp, "skyarea": src_backscal * (1.0 / 20.0 / 60.0) ** 2, } def _mean_arf(arf_path: str, lo: float, hi: float) -> float: with fits.open(arf_path, memmap=False) as hdul: tab = hdul["SPECRESP"].data elo = np.asarray(tab["ENERG_LO"], dtype=float) ehi = np.asarray(tab["ENERG_HI"], dtype=float) area = np.asarray(tab["SPECRESP"], dtype=float) overlap = np.clip(np.minimum(ehi, hi) - np.maximum(elo, lo), 0.0, None) if overlap.sum() <= 0: return float("nan") return float(np.sum(area * overlap) / np.sum(overlap)) def _band_channel_mask( e_min: np.ndarray, e_max: np.ndarray, lo: float, hi: float ) -> np.ndarray: """Assign each RMF channel to at most one half-open energy interval.""" e_mid = 0.5 * (np.asarray(e_min, dtype=float) + np.asarray(e_max, dtype=float)) return (e_mid >= lo) & (e_mid < hi) def compute_band_statistics(manifest: list[dict[str, Any]]) -> list[dict[str, Any]]: """Compute counts-space QPB-subtracted statistics for explicit bands.""" rows: list[dict[str, Any]] = [] for spec in manifest: arrays = _pha_arrays(spec["pha"], spec["bkg"], spec["rmf"]) src = np.asarray(arrays["src"]) bkg = np.asarray(arrays["bkg"]) bkg_err = np.asarray(arrays["bkg_err"]) e_min = np.asarray(arrays["e_min"]) e_max = np.asarray(arrays["e_max"]) alpha = float(arrays["alpha"]) exposure = float(arrays["exposure"]) skyarea = float(arrays["skyarea"]) for lo, hi in BANDS: select = _band_channel_mask(e_min, e_max, lo, hi) source_counts = float(np.sum(src[select])) qpb_counts = float(alpha * np.sum(bkg[select])) net_counts = source_counts - qpb_counts variance = source_counts + float(alpha**2 * np.sum(bkg_err[select] ** 2)) snr = net_counts / np.sqrt(variance) if variance > 0 else float("nan") arf = _mean_arf(spec["arf"], lo, hi) band_width = hi - lo net_rate = net_counts / exposure net_sb = net_rate / skyarea net_sb_flux_density = net_sb / arf / band_width if arf > 0 else float("nan") rows.append( { "detector": spec["label"], "instrument": spec["instrument"], "band": f"{lo:g}-{hi:g}", "energy_lo_keV": lo, "energy_hi_keV": hi, "source_counts": source_counts, "scaled_qpb_counts": qpb_counts, "net_counts": net_counts, "variance": variance, "snr": snr, "qpb_fraction": qpb_counts / source_counts if source_counts > 0 else float("nan"), "exposure_s": exposure, "skyarea_arcmin2": skyarea, "mean_arf_cm2": arf, "net_rate_cts_s": net_rate, "net_rate_per_arcmin2": net_sb, "approx_net_sb_flux_density": net_sb_flux_density, "background_scale_alpha": alpha, } ) return rows def _write_statistics(rows: list[dict[str, Any]]) -> Path: OUT.mkdir(parents=True, exist_ok=True) path = OUT / "r015_epic_mos1_mos2_pn0_band_statistics.csv" with path.open("w", newline="") as handle: writer = csv.DictWriter(handle, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) return path def _gain_summary(rows: list[dict[str, Any]]) -> dict[str, Any]: summary: dict[str, Any] = {} for band in ("0.4-3.2", "0.7-1.05", "0.7-0.875", "0.875-1.05"): selected = {row["detector"]: row for row in rows if row["band"] == band} mos = [selected["MOS1"], selected["MOS2"]] pn = selected["PN0 (PATTERN=0)"] mos_net = sum(row["net_counts"] for row in mos) mos_var = sum(row["variance"] for row in mos) mos_snr = mos_net / np.sqrt(mos_var) summary[band] = { "mos1": selected["MOS1"], "mos2": selected["MOS2"], "pn0": pn, "mos12_combined_net_counts": mos_net, "mos12_combined_snr": mos_snr, "pn0_to_mos12_net_count_ratio": pn["net_counts"] / mos_net, "pn0_to_mos12_snr_ratio": pn["snr"] / mos_snr, "mos12_plus_pn0_snr": (mos_net + pn["net_counts"]) / np.sqrt(mos_var + pn["variance"]), "mos12_plus_pn0_snr_gain": ( (mos_net + pn["net_counts"]) / np.sqrt(mos_var + pn["variance"]) ) / mos_snr, } low = summary["0.7-0.875"] high = summary["0.875-1.05"] for label, key in (("MOS1", "mos1"), ("MOS2", "mos2"), ("PN0 (PATTERN=0)", "pn0")): summary.setdefault("fe_l_high_to_low", {})[label] = ( high[key]["approx_net_sb_flux_density"] / low[key]["approx_net_sb_flux_density"] ) return summary def generate_inspect(manifest: list[dict[str, Any]]) -> dict[str, dict[str, int]]: """Run plot-only xsherpa SpectrumInspector and return grouped-bin counts.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import sherpa.astro.ui as ui from xsherpa.inspect import SpectrumInspector try: plt.style.use(["science", "no-latex"]) except OSError: pass OUT.mkdir(parents=True, exist_ok=True) ui.clean() inspector = SpectrumInspector( output_dir=str(OUT), colors=["#4C78A8", "#F58518", "#D62728"] ) spectra = [ {key: row[key] for key in ("pha", "bkg", "arf", "rmf")} for row in manifest ] inspector.load_spectra( spectra, labels=[row["label"] for row in manifest], energy_range=ENERGY_RANGE, ) if len(inspector.spectra) != 3: raise RuntimeError(f"xsherpa loaded {len(inspector.spectra)}/3 spectra") variants = {"1col": (3.35, 2.7), "2col": (7.0, 4.5)} plot_sets = ( ("r015_epic_inspect", ENERGY_RANGE), ("r015_epic_felzoom_inspect", FEL_ZOOM), ) for stem, xlim in plot_sets: for unit in ("rate", "sb", "sb_flux"): for width, figsize in variants.items(): for ext in ("png", "pdf"): with np.errstate(divide="ignore", invalid="ignore"): inspector.plot_spectra( output_filename=f"{stem}_{unit}_{width}_ylog.{ext}", flux_unit=unit, ylog=True, xlim=xlim, ylim=SB_FLUX_YLIM if unit == "sb_flux" else None, figsize=figsize, dpi=300, sampling_rate=SAMPLING_RATE, snr=SNR, ) grouped: dict[str, dict[str, int]] = {} for row, spec in zip(manifest, inspector.spectra): plot = ui.get_data_plot(spec["spec_id"]) xlo = np.asarray(plot.xlo, dtype=float) xhi = np.asarray(plot.xhi, dtype=float) grouped[row["label"]] = { "total_0.4_3.2_keV": int(xlo.size), "overlap_0.6_1.2_keV": int( np.sum((xhi > FEL_ZOOM[0]) & (xlo < FEL_ZOOM[1])) ), } return grouped def make_footprint_overlay(manifest: list[dict[str, Any]]) -> list[str]: """Show existing camera-specific count images, masks, and nominal R0-15 circle.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import astropy.units as u from astropy.coordinates import SkyCoord from astropy.visualization.wcsaxes import SphericalCircle from astropy.wcs import WCS from scipy.ndimage import gaussian_filter try: plt.style.use(["science", "no-latex"]) except OSError: pass outputs: list[str] = [] variants = {"1col": (3.35, 8.0), "2col": (7.0, 2.8)} center = SkyCoord(CENTER_RA_DEG * u.deg, CENTER_DEC_DEG * u.deg, frame="fk5") for width, figsize in variants.items(): fig = plt.figure(figsize=figsize) for i, row in enumerate(manifest, start=1): image = np.asarray(fits.getdata(row["image"]), dtype=float) mask = np.asarray(fits.getdata(row["mask"]), dtype=float) wcs = WCS(fits.getheader(row["image"])).celestial ax = fig.add_subplot(3, 1, i, projection=wcs) if width == "1col" else fig.add_subplot(1, 3, i, projection=wcs) smoothed = gaussian_filter(image, 1.5) positive = smoothed[smoothed > 0] vmax = float(np.percentile(positive, 99.5)) if positive.size else 1.0 ax.imshow(smoothed, origin="lower", cmap="magma", vmin=0.0, vmax=max(vmax, 1e-6)) included = np.ma.masked_where(mask <= 0, mask) ax.imshow(included, origin="lower", cmap="Blues", alpha=0.16, vmin=0, vmax=1) ax.contour(mask, levels=[0.5], colors="#4DD0E1", linewidths=0.45, alpha=0.9) circle = SphericalCircle( (center.ra, center.dec), 15.0 * u.arcmin, transform=ax.get_transform("fk5"), edgecolor="white", facecolor="none", linestyle="--", linewidth=0.8, ) ax.add_patch(circle) ax.scatter( [center.ra.deg], [center.dec.deg], transform=ax.get_transform("fk5"), marker="+", s=18, linewidths=0.8, color="lime", ) ax.set_title( f"{row['label']} | PHA area={row['skyarea_arcmin2']:.1f} arcmin$^2$", fontsize=7, ) ax.coords[0].set_axislabel("RA", fontsize=6) ax.coords[1].set_axislabel("Dec", fontsize=6) ax.tick_params(labelsize=5) fig.suptitle( "M104 R0-15 camera coverage | 0.5-2.0 keV counts | " "cyan=broad detector mask grid (not PHA BACKSCAL extraction area), white=15'", fontsize=8, ) fig.tight_layout() for ext in ("png", "pdf"): path = OUT / f"r015_epic_camera_footprints_region_mask_overlay_{width}.{ext}" fig.savefig(path, dpi=300, bbox_inches="tight") outputs.append(str(path)) plt.close(fig) return outputs def main() -> None: manifest = validate_inputs() rows = compute_band_statistics(manifest) stats_path = _write_statistics(rows) grouped_bins = generate_inspect(manifest) overlays = make_footprint_overlay(manifest) summary = { "purpose": "No-fit first-look R0-15 EPIC camera comparison", "qualification": ( "The MOS R0-15 directories are compatibility symlinks to the current 4.background masked full-FOV products; " "PN0 uses the analogous legacy masked camera-footprint selection with PATTERN<=0 and an approximately 15-arcmin detector circle. " "The sky regions overlap the nominal R0-15 field but are not pixel-identical because each camera has different chip gaps and masks. " "Use sb/sb_flux for shape comparison; rate also contains camera area differences." ), "energy_range_keV": list(ENERGY_RANGE), "fe_l_zoom_keV": list(FEL_ZOOM), "grouping": {"method": "group_subtracted_oversampling_SNR", "sampling_rate": SAMPLING_RATE, "snr": SNR}, "fit_performed": False, "manifest": manifest, "grouped_bins": grouped_bins, "band_statistics_csv": str(stats_path.resolve()), "comparison": _gain_summary(rows), "overlay_outputs": [str(Path(path).resolve()) for path in overlays], } summary_path = OUT / "r015_epic_mos1_mos2_pn0_summary.json" summary_path.write_text(json.dumps(summary, indent=2, ensure_ascii=False, allow_nan=False) + "\n") print(json.dumps({"output_dir": str(OUT.resolve()), "grouped_bins": grouped_bins, "summary": str(summary_path.resolve())}, ensure_ascii=False)) if __name__ == "__main__": main()