#!/usr/bin/env python3 """Compare the previous mixed R0-15 EPIC contract with traditional 3B full-region spectra.""" from __future__ import annotations import csv import hashlib import json from pathlib import Path from typing import Any import numpy as np from astropy.io import fits import inspect_r015_epic_mos12_pn0 as previous ROOT = Path(__file__).resolve().parent OUT = ( ROOT / "joint_spectrum_fitting_2T_basedon_region_v22" / "r015_epic_3b_fullregion_xsherpa_comparison_20260713" ) ENERGY_RANGE = previous.ENERGY_RANGE FEL_ZOOM = previous.FEL_ZOOM SB_FLUX_YLIM = previous.SB_FLUX_YLIM SAMPLING_RATE = previous.SAMPLING_RATE SNR = previous.SNR THREE_B_SPECS = ( { "label": "MOS1", "instrument": "EMOS1", "pha": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_mos1_00500-02000_exclude_extent_source_mask/mos1U005-obj.pi", }, { "label": "MOS2", "instrument": "EMOS2", "pha": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_mos2_00500-02000_exclude_extent_source_mask/mos2U005-obj.pi", }, { "label": "PN0 (PATTERN=0)", "instrument": "EPN", "pha": ROOT / "data/0900170101/rpc/3B.background_20250124/bkg_pn0_00500-02000_mask/pnU002-obj.pi", }, ) def _companions(pha: Path) -> dict[str, Path]: stem = pha.name.removesuffix("-obj.pi") return { "pha": pha, "bkg": pha.parent / f"{stem}-back.pi", "arf": pha.parent / f"{stem}.arf", "rmf": pha.parent / f"{stem}.rmf", } def _sha256(path: str | Path) -> str: digest = hashlib.sha256() with Path(path).open("rb") as handle: for block in iter(lambda: handle.read(1024 * 1024), b""): digest.update(block) return digest.hexdigest() def validate_3b_inputs() -> list[dict[str, Any]]: """Validate the traditional 3B masked full-region three-camera contract.""" rows: list[dict[str, Any]] = [] for spec in THREE_B_SPECS: companions = _companions(Path(spec["pha"])) for role, path in companions.items(): if not path.is_file() or path.stat().st_size <= 0: raise FileNotFoundError(f"Missing or empty {role}: {path}") header = fits.getheader(companions["pha"], "SPECTRUM") expression = str(header.get("SLCTEXPR", "")) if str(header.get("INSTRUME")) != spec["instrument"]: raise ValueError(f"Instrument mismatch: {companions['pha']}") if spec["instrument"] == "EPN": if "PATTERN <= 0" not in expression or "FLAG == 0" not in expression: raise ValueError("Traditional PN input is not verified pn0") elif "PATTERN<=12" not in expression or "FLAG == 0" not in expression: raise ValueError(f"Unexpected MOS selection: {companions['pha']}") rows.append( { **spec, **{role: str(path.absolute()) for role, path in companions.items()}, "lineage": "3B.background_20250124", "region_contract": "traditional masked camera full region", "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"]), "selection_expression": expression, "selection_expression_status": ( "Pattern/FLAG provenance is usable; long spatial SLCTEXPR continuation is not treated as complete geometry provenance" ), } ) return rows def build_input_comparison() -> dict[str, dict[str, Any]]: old = {row["label"]: row for row in previous.validate_inputs()} new = {row["label"]: row for row in validate_3b_inputs()} result: dict[str, dict[str, Any]] = {} for label in old: row: dict[str, Any] = { "previous_lineage": ( "4.background" if label != "PN0 (PATTERN=0)" else "3B.background_20250124" ), "traditional_3b_lineage": "3B.background_20250124", "previous_exposure_s": old[label]["exposure_s"], "traditional_3b_exposure_s": new[label]["exposure_s"], "previous_skyarea_arcmin2": old[label]["skyarea_arcmin2"], "traditional_3b_skyarea_arcmin2": new[label]["skyarea_arcmin2"], } for role in ("pha", "bkg", "arf", "rmf"): old_hash = _sha256(old[label][role]) new_hash = _sha256(new[label][role]) row[f"previous_{role}"] = old[label][role] row[f"traditional_3b_{role}"] = new[label][role] row[f"previous_{role}_sha256"] = old_hash row[f"traditional_3b_{role}_sha256"] = new_hash row[f"same_{role}_sha256"] = old_hash == new_hash result[label] = row return result def compute_contract_comparison() -> dict[str, Any]: old_rows = previous.compute_band_statistics(previous.validate_inputs()) new_rows = previous.compute_band_statistics(validate_3b_inputs()) result: dict[str, Any] = {"bands": {}, "differences": {}} for band in sorted({row["band"] for row in old_rows}): old_band = {row["detector"]: row for row in old_rows if row["band"] == band} new_band = {row["detector"]: row for row in new_rows if row["band"] == band} result["bands"][band] = { "previous_mixed": old_band, "traditional_3b": new_band, } result["differences"][band] = {} for label in old_band: old = old_band[label] new = new_band[label] result["differences"][band][label] = { "exposure_ratio_3b_to_previous": new["exposure_s"] / old["exposure_s"], "skyarea_ratio_3b_to_previous": new["skyarea_arcmin2"] / old["skyarea_arcmin2"], "source_counts_ratio_3b_to_previous": new["source_counts"] / old["source_counts"], "scaled_qpb_ratio_3b_to_previous": new["scaled_qpb_counts"] / old["scaled_qpb_counts"], "net_counts_ratio_3b_to_previous": new["net_counts"] / old["net_counts"], "snr_ratio_3b_to_previous": new["snr"] / old["snr"], "net_rate_per_arcmin2_ratio_3b_to_previous": ( new["net_rate_per_arcmin2"] / old["net_rate_per_arcmin2"] ), "approx_sb_flux_ratio_3b_to_previous": ( new["approx_net_sb_flux_density"] / old["approx_net_sb_flux_density"] ), "mean_arf_ratio_3b_to_previous": new["mean_arf_cm2"] / old["mean_arf_cm2"], } result["previous_summary"] = previous._gain_summary(old_rows) result["traditional_3b_summary"] = previous._gain_summary(new_rows) result["fe_l_color_comparison"] = {} for contract in ("previous_mixed", "traditional_3b"): colors: dict[str, dict[str, float]] = {} for label in ("MOS1", "MOS2", "PN0 (PATTERN=0)"): low = result["bands"]["0.7-0.875"][contract][label] high = result["bands"]["0.875-1.05"][contract][label] ratio = ( high["approx_net_sb_flux_density"] / low["approx_net_sb_flux_density"] ) sigma = ratio * np.sqrt( high["variance"] / high["net_counts"] ** 2 + low["variance"] / low["net_counts"] ** 2 ) colors[label] = {"ratio": float(ratio), "sigma": float(sigma)} contract_summary: dict[str, Any] = {"detectors": colors} pn = colors["PN0 (PATTERN=0)"] for label in ("MOS1", "MOS2"): delta = colors[label]["ratio"] - pn["ratio"] delta_sigma = float( np.hypot(colors[label]["sigma"], pn["sigma"]) ) contract_summary[f"{label}_vs_PN0"] = { "delta": float(delta), "delta_sigma": delta_sigma, "z_score": float(delta / delta_sigma), } result["fe_l_color_comparison"][contract] = contract_summary return result def _run_3b_inspector(manifest: list[dict[str, Any]]) -> dict[str, dict[str, int]]: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import sherpa.astro.ui as ui # type: ignore[import-not-found] from xsherpa.inspect import SpectrumInspector # type: ignore[import-not-found] try: plt.style.use(["science", "no-latex"]) except OSError: pass ui.clean() inspector = SpectrumInspector( output_dir=OUT, colors=["#0072B2", "#E69F00", "#D55E00"], ) inspector.load_spectra( [ {"pha": row["pha"], "bkg": row["bkg"], "arf": row["arf"], "rmf": row["rmf"]} for row in manifest ], labels=[row["label"] for row in manifest], energy_range=ENERGY_RANGE, ) variants = {"1col": (3.35, 2.7), "2col": (7.0, 4.5)} ranges = { "r015_3b_fullregion_inspect": ENERGY_RANGE, "r015_3b_fullregion_felzoom_inspect": FEL_ZOOM, } for stem, xlim in ranges.items(): 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, loaded in zip(manifest, inspector.spectra): plot = ui.get_data_plot(loaded["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_stacked_comparisons() -> list[str]: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt try: plt.style.use(["science", "no-latex"]) except OSError: pass panels = { "r015_previous_vs_3b_fullregion_sb_flux_stacked": ( previous.OUT / "r015_epic_inspect_sb_flux_2col_ylog.png", OUT / "r015_3b_fullregion_inspect_sb_flux_2col_ylog.png", "R0-15 EPIC broad comparison | sb_flux | 0.4-3.2 keV | no fit", ), "r015_previous_vs_3b_fullregion_felzoom_sb_flux_stacked": ( previous.OUT / "r015_epic_felzoom_inspect_sb_flux_2col_ylog.png", OUT / "r015_3b_fullregion_felzoom_inspect_sb_flux_2col_ylog.png", "R0-15 EPIC Fe-L comparison | sb_flux | 0.6-1.2 keV | no fit", ), } outputs: list[str] = [] for stem, (old_path, new_path, title) in panels.items(): if not old_path.is_file() or not new_path.is_file(): raise FileNotFoundError(f"Missing stacked input: {old_path} or {new_path}") images = [plt.imread(old_path), plt.imread(new_path)] for width, figsize in {"1col": (3.35, 5.2), "2col": (7.0, 8.2)}.items(): fig, axes = plt.subplots(2, 1, figsize=figsize) for ax, image, panel_title in zip( axes, images, [ "A. Previous mixed contract: 4B MOS + 3B pn0", "B. Traditional 3B masked full region: MOS1 + MOS2 + pn0", ], ): ax.imshow(image) ax.set_axis_off() ax.set_title(panel_title, fontsize=7 if width == "1col" else 9, loc="left") fig.suptitle(title, fontsize=7 if width == "1col" else 10) fig.tight_layout() for ext in ("png", "pdf"): path = OUT / f"{stem}_{width}.{ext}" fig.savefig(path, dpi=300, bbox_inches="tight") outputs.append(str(path)) plt.close(fig) return outputs def _write_csv(comparison: dict[str, Any], path: Path) -> None: rows: list[dict[str, Any]] = [] for band, contracts in comparison["bands"].items(): for contract, detectors in contracts.items(): for detector, values in detectors.items(): rows.append({"contract": contract, "band": band, "detector": detector, **values}) with path.open("w", newline="") as handle: writer = csv.DictWriter(handle, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) def main() -> None: OUT.mkdir(parents=True, exist_ok=True) manifest = validate_3b_inputs() inventory = build_input_comparison() comparison = compute_contract_comparison() grouped = _run_3b_inspector(manifest) stacked = _make_stacked_comparisons() summary = { "purpose": "Previous mixed R0-15 versus traditional 3B masked full-region no-fit comparison", "energy_range_keV": ENERGY_RANGE, "fe_l_zoom_keV": FEL_ZOOM, "grouping": {"method": "group_subtracted_oversampling_SNR", "sampling_rate": SAMPLING_RATE, "snr": SNR}, "traditional_3b_manifest": manifest, "input_comparison": inventory, "contract_comparison": comparison, "traditional_3b_grouped_bins": grouped, "stacked_outputs": stacked, "qualification": ( "No fit. QPB-subtracted first-look only. Traditional 3B and previous 4B MOS processing differ; PN0 is a byte-identical control." ), } (OUT / "r015_previous_vs_3b_input_inventory.json").write_text( json.dumps(inventory, indent=2, allow_nan=False) + "\n" ) _write_csv(comparison, OUT / "r015_previous_vs_3b_band_statistics.csv") (OUT / "r015_previous_vs_3b_summary.json").write_text( json.dumps(summary, indent=2, allow_nan=False) + "\n" ) print(json.dumps({"output_dir": str(OUT), "grouped_bins": grouped}, indent=2)) if __name__ == "__main__": main()