File indexing completed on 2026-10-08 08:36:03
0001
0002 """Compare the FullSet G1 reproduction and assemble calibration PDF plot books.
0003
0004 Run from examples/yall/fullset-g1-repro after:
0005 source env.tcsh
0006 python3 compare_and_report.py
0007
0008 Uses only Python's standard library plus pdfunite.
0009 """
0010
0011 from __future__ import annotations
0012
0013 import argparse
0014 import csv
0015 import json
0016 import math
0017 import os
0018 import re
0019 import shutil
0020 import statistics
0021 import subprocess
0022 import sys
0023 from decimal import Decimal
0024 from pathlib import Path
0025
0026 COLUMNS = (
0027 "cell_id", "layer", "row", "column", "module",
0028 "ped_mean_h", "ped_sigma_h", "ped_mean_l", "ped_sigma_l",
0029 "mip_scale_h", "mip_width_h", "mip_scale_l", "mip_width_l",
0030 "lghg_corr", "lghg_corr_offset", "hglg_corr",
0031 "hglg_corr_offset_toa", "bc",
0032 )
0033 INT_FIELDS = {"cell_id", "layer", "row", "column", "module", "bc"}
0034 STAGES = ("pedestal", "mip", "refine1", "refine2", "refine3", "refine4", "refine5", "final")
0035
0036
0037 def args():
0038 here = Path(__file__).resolve().parent
0039 repo = here.parents[2]
0040 p = argparse.ArgumentParser(description=__doc__)
0041 p.add_argument("--work", type=Path, default=None,
0042 help="FullSet G1 workflow output directory")
0043 p.add_argument("--reference", type=Path,
0044 default=repo / "calibrations/TB2026/calib_SPS-H2_FullSetG_1.txt")
0045 p.add_argument("--set-name", default="FullSetG_1")
0046 p.add_argument("--out", type=Path, default=None)
0047 p.add_argument("--no-pdf", action="store_true")
0048 return p.parse_args()
0049
0050
0051 def find_work(explicit):
0052 """Find a completed G1 output tree without trusting stale example state."""
0053 if explicit is not None:
0054 return explicit
0055
0056 expected = "fullset-g1-repro"
0057 final_rel = Path("final/calib_Final_Muon_FullSetG_1_calib.txt")
0058 candidates = []
0059
0060
0061 if os.environ.get("LFHCAL_EXAMPLE") == expected:
0062 value = os.environ.get("LFHCAL_EXAMPLE_WORK")
0063 if value:
0064 candidates.append(Path(value))
0065
0066
0067 value = os.environ.get("LFHCAL_WORK")
0068 if value:
0069 candidates.append(Path(value) / expected)
0070
0071
0072
0073 user = os.environ.get("USER") or os.environ.get("LOGNAME")
0074 if user:
0075 candidates.append(Path(f"/gpfs01/star/scratch/{user}/lfhcal") / expected)
0076
0077 seen = set()
0078 for candidate in candidates:
0079 candidate = candidate.expanduser()
0080 key = str(candidate)
0081 if key in seen:
0082 continue
0083 seen.add(key)
0084 if (candidate / final_rel).is_file():
0085 return candidate
0086
0087 tried = "\n ".join(str(x) for x in candidates) or "(no candidates)"
0088 raise ValueError(
0089 "could not locate a completed FullSet G1 output tree.\n"
0090 "Pass it explicitly, for example:\n"
0091 " python3 compare_and_report.py "
0092 "--work /gpfs01/star/scratch/$USER/lfhcal/fullset-g1-repro\n"
0093 f"Tried:\n {tried}"
0094 )
0095
0096
0097 def parse_calib(path):
0098 rows = {}
0099 payload = path.read_bytes()
0100 for line_no, line in enumerate(payload.decode("utf-8").splitlines(), 1):
0101 s = line.strip()
0102 if not s or s.startswith(("#", "%")) or s.startswith("RunNr:"):
0103 continue
0104 fields = s.split()
0105 if not fields or not fields[0].isdigit():
0106 continue
0107 if len(fields) != len(COLUMNS):
0108 raise ValueError(f"{path}:{line_no}: expected 18 columns, got {len(fields)}")
0109 row = {}
0110 for name, token in zip(COLUMNS, fields):
0111 value = Decimal(token)
0112 if name in INT_FIELDS:
0113 if value != int(value):
0114 raise ValueError(f"{path}:{line_no}: non-integer {name}")
0115 value = int(value)
0116 row[name] = value
0117 cell = row["cell_id"]
0118 if cell in rows:
0119 raise ValueError(f"{path}: duplicate cell {cell}")
0120 rows[cell] = row
0121 if not rows:
0122 raise ValueError(f"{path}: no calibration records")
0123 return rows, payload
0124
0125
0126 def one_calib(directory):
0127 files = sorted(directory.glob("*_calib.txt"))
0128 if len(files) != 1:
0129 raise ValueError(f"expected one *_calib.txt in {directory}, found {len(files)}")
0130 return files[0]
0131
0132
0133 def pct(ref, got):
0134 return float(100 * (got - ref) / ref) if ref > 0 and got > 0 else None
0135
0136
0137 def percentile(values, q):
0138 if not values:
0139 return None
0140 x = sorted(values)
0141 pos = (len(x) - 1) * q
0142 lo, hi = math.floor(pos), math.ceil(pos)
0143 return x[lo] + (x[hi] - x[lo]) * (pos - lo)
0144
0145
0146 def stats(values):
0147 if not values:
0148 return {"n": 0}
0149 a = [abs(x) for x in values]
0150 return {
0151 "n": len(values),
0152 "median_abs_pct": statistics.median(a),
0153 "p95_abs_pct": percentile(a, .95),
0154 "rms_pct": math.sqrt(statistics.mean(x*x for x in values)),
0155 "max_abs_pct": max(a),
0156 }
0157
0158
0159 def write_csv(path, rows):
0160 if not rows:
0161 return
0162 fields = list(dict.fromkeys(k for row in rows for k in row))
0163 with path.open("w", newline="", encoding="utf-8") as f:
0164 w = csv.DictWriter(f, fieldnames=fields)
0165 w.writeheader()
0166 w.writerows(rows)
0167
0168
0169 def compare(work, reference_path, out, set_name):
0170 ref, ref_bytes = parse_calib(reference_path)
0171 tables = {"reference": ref}
0172 payloads = {"reference": ref_bytes}
0173 paths = {"reference": reference_path}
0174
0175 for stage in STAGES:
0176 path = one_calib(work / stage)
0177 tables[stage], payloads[stage] = parse_calib(path)
0178 paths[stage] = path
0179
0180 ids = sorted(ref)
0181 for stage, table in tables.items():
0182 if set(table) != set(ref):
0183 raise ValueError(f"{stage}: cell-ID set differs from reference")
0184 bad = [i for i in ids if any(
0185 table[i][x] != ref[i][x] for x in ("layer", "row", "column", "module")
0186 )]
0187 if bad:
0188 raise ValueError(f"{stage}: geometry mismatch, first cells {bad[:10]}")
0189
0190 final = tables["final"]
0191 final_rows = []
0192 for i in ids:
0193 a, b = ref[i], final[i]
0194 row = {k: a[k] for k in ("cell_id", "module", "layer", "row", "column")}
0195 for field in ("mip_scale_h", "mip_width_h",
0196 "ped_mean_h", "ped_sigma_h", "ped_mean_l", "ped_sigma_l"):
0197 row[f"{field}_reference"] = a[field]
0198 row[f"{field}_reproduced"] = b[field]
0199 row[f"{field}_delta"] = b[field] - a[field]
0200 if field.startswith("mip_"):
0201 row[f"{field}_pct_delta"] = pct(a[field], b[field])
0202 final_rows.append(row)
0203
0204 def metric(table, field):
0205 d = {i: pct(ref[i][field], table[i][field]) for i in ids}
0206 d = {i: v for i, v in d.items() if v is not None}
0207 s = stats(list(d.values()))
0208 if d:
0209 s["max_pct_cell"] = max(d, key=lambda i: abs(d[i]))
0210 s["exact_equal"] = sum(ref[i][field] == table[i][field] for i in d)
0211 s["within_0.001_pct"] = sum(abs(v) <= .001 for v in d.values())
0212 s["within_0.01_pct"] = sum(abs(v) <= .01 for v in d.values())
0213 return s
0214
0215 scale = metric(final, "mip_scale_h")
0216 width = metric(final, "mip_width_h")
0217
0218 stage_rows = []
0219 for stage in ("mip", "refine1", "refine2", "refine3", "refine4", "refine5", "final"):
0220 s = metric(tables[stage], "mip_scale_h")
0221 stage_rows.append({"stage": stage, **s})
0222
0223 ranked_scale = []
0224 for i in ids:
0225 d = pct(ref[i]["mip_scale_h"], final[i]["mip_scale_h"])
0226 if d is not None:
0227 ranked_scale.append((abs(d), d, i))
0228 ranked_scale.sort(reverse=True)
0229
0230 unavailable = [i for i in ids if ref[i]["mip_scale_h"] <= 0 or final[i]["mip_scale_h"] <= 0]
0231 restored = [i for i in ids
0232 if tables["mip"][i]["mip_scale_h"] <= 0
0233 and tables["refine1"][i]["mip_scale_h"] > 0]
0234 bc_mismatch = [i for i in ids if ref[i]["bc"] != final[i]["bc"]]
0235
0236 comparison = out / "comparison"
0237 comparison.mkdir(parents=True, exist_ok=True)
0238 write_csv(comparison / "final_comparison.csv", final_rows)
0239 write_csv(comparison / "stage_summary.csv", stage_rows)
0240
0241 summary = {
0242 "set": set_name,
0243 "reference": str(reference_path),
0244 "reproduced": str(paths["final"]),
0245 "cells": len(ids),
0246 "mip_scale_h": scale,
0247 "mip_width_h": width,
0248 "bc_mismatch_ids": bc_mismatch,
0249 "unavailable_mip_ids": unavailable,
0250 "restored_at_refine1_ids": restored,
0251 "final_identical_to_refine5_bytes": payloads["final"] == payloads["refine5"],
0252 }
0253 (comparison / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
0254
0255 lines = [
0256 f"{set_name} calibration comparison",
0257 f"Reference: {reference_path}",
0258 f"Reproduced: {paths['final']}",
0259 f"Cells: {len(ids)}",
0260 "",
0261 "HG MIP scale",
0262 f" positive pairs: {scale['n']}",
0263 f" exact at exported precision: {scale.get('exact_equal', 0)}",
0264 f" median |relative difference|: {scale['median_abs_pct']:.9f}%",
0265 f" 95th percentile: {scale['p95_abs_pct']:.9f}%",
0266 f" RMS: {scale['rms_pct']:.9f}%",
0267 f" maximum: {scale['max_abs_pct']:.9f}% at cell {scale.get('max_pct_cell')}",
0268 f" within 0.001%: {scale.get('within_0.001_pct', 0)}/{scale['n']}",
0269 f" within 0.01%: {scale.get('within_0.01_pct', 0)}/{scale['n']}",
0270 "",
0271 "Stored HG MIP width",
0272 f" median |relative difference|: {width['median_abs_pct']:.9f}%",
0273 f" maximum: {width['max_abs_pct']:.9f}% at cell {width.get('max_pct_cell')}",
0274 "",
0275 f"BC mismatches: {bc_mismatch}",
0276 f"Unavailable/nonpositive HG MIP cells: {unavailable}",
0277 f"Recovered at refine1: {restored}",
0278 f"Final byte-identical to refine5: {summary['final_identical_to_refine5_bytes']}",
0279 "",
0280 "Convergence vs reference (HG MIP scale)",
0281 ]
0282 for row in stage_rows:
0283 lines.append(
0284 f" {row['stage']:7s} n={row['n']:3d} "
0285 f"median |d|={row['median_abs_pct']:.9f}% "
0286 f"max |d|={row['max_abs_pct']:.9f}%"
0287 )
0288
0289 lines += ["", "Largest final HG MIP-scale residuals"]
0290 for _, d, i in ranked_scale[:15]:
0291 a, b = ref[i], final[i]
0292 lines.append(
0293 f" cell {i:4d} mod={a['module']} layer={a['layer']} "
0294 f"row={a['row']} col={a['column']} "
0295 f"ref={float(a['mip_scale_h']):.6f} "
0296 f"repro={float(b['mip_scale_h']):.6f} "
0297 f"d={d:+.9f}%"
0298 )
0299
0300 if len(ranked_scale) > 1:
0301 ss = sum(d*d for _, d, _ in ranked_scale)
0302 rms_without_worst = math.sqrt((ss - ranked_scale[0][1]**2) / (len(ranked_scale) - 1))
0303 lines += [
0304 "",
0305 f"HG MIP-scale RMS excluding the single worst cell: {rms_without_worst:.9f}%",
0306 ]
0307
0308 text = "\n".join(lines) + "\n"
0309 (comparison / "comparison.txt").write_text(text)
0310 print(text, end="")
0311
0312
0313 def natural(paths):
0314 def key(path):
0315 return [int(x) if x.isdigit() else x.lower()
0316 for x in re.split(r"(\d+)", path.name)]
0317 return sorted(paths, key=key)
0318
0319
0320 def exact(directory, names):
0321 files, missing = [], []
0322 for name in names:
0323 path = directory / name
0324 (files if path.is_file() else missing).append(path if path.is_file() else name)
0325 return files, missing
0326
0327
0328 def mip_inputs(directory, initial):
0329 suffix = "_2nd" if initial else ""
0330 head = [
0331 f"HG_FWHMMip{suffix}.pdf",
0332 f"HG_GaussSigMip{suffix}.pdf",
0333 f"HG_LandMPVMip{suffix}.pdf",
0334 f"HG_LandSigMip{suffix}.pdf",
0335 f"HG_MaxMip{suffix}.pdf",
0336 "HGscaleChi2VsLayer.pdf",
0337 ]
0338 tail = [
0339 "MipTriggXY.pdf",
0340 "MuonTriggers.pdf",
0341 "SNRTriggVsLayer.pdf",
0342 "SuppressionNoise.pdf",
0343 "SuppressionSignal.pdf",
0344 ]
0345 files, missing = exact(directory, head)
0346 mip_layers = natural(directory.glob("MIP_HG_Layer*.pdf"))
0347 if mip_layers:
0348 files += mip_layers
0349 else:
0350 missing.append("MIP_HG_Layer*.pdf")
0351 got, miss = exact(directory, tail)
0352 files += got
0353 missing += miss
0354 trig_layers = natural(directory.glob("TriggPrimitive_Layer*.pdf"))
0355 if trig_layers:
0356 files += trig_layers
0357 else:
0358 missing.append("TriggPrimitive_Layer*.pdf")
0359 return files, missing
0360
0361
0362 def unite(output, inputs):
0363 if not inputs:
0364 return
0365 output.parent.mkdir(parents=True, exist_ok=True)
0366 if output.exists():
0367 output.unlink()
0368 if len(inputs) == 1:
0369 shutil.copyfile(inputs[0], output)
0370 return
0371 exe = shutil.which("pdfunite")
0372 if exe:
0373 subprocess.run([exe, *map(str, inputs), str(output)], check=True)
0374 return
0375
0376 exe = shutil.which("qpdf")
0377 if exe:
0378 subprocess.run(
0379 [exe, "--empty", "--pages", *map(str, inputs), "--", str(output)],
0380 check=True,
0381 )
0382 return
0383
0384 exe = shutil.which("gs")
0385 if exe:
0386 subprocess.run(
0387 [
0388 exe, "-q", "-dBATCH", "-dNOPAUSE", "-sDEVICE=pdfwrite",
0389 f"-sOutputFile={output}", *map(str, inputs),
0390 ],
0391 check=True,
0392 )
0393 return
0394
0395 raise RuntimeError("no PDF merger found; need pdfunite, qpdf, or gs in PATH")
0396
0397
0398 def make_pdfs(work, out, set_name):
0399 report_dir = out / "pdf"
0400 report_dir.mkdir(parents=True, exist_ok=True)
0401 plots = work / "plots"
0402 muon = f"Muon_{set_name}"
0403
0404 stages = [
0405 ("mip", "Initial", True),
0406 ("refine1", "ImpR", False),
0407 ("refine2", "Imp2R", False),
0408 ("refine3", "Imp3R", False),
0409 ("refine4", "Imp4R", False),
0410 ("refine5", "Imp5R", False),
0411 ]
0412 made = {}
0413 for stage, label, initial in stages:
0414 directory = plots / stage / muon
0415 if not directory.is_dir():
0416 print(f"PDF: skip {stage}; no {directory}", file=sys.stderr)
0417 continue
0418 files, missing = mip_inputs(directory, initial)
0419 if missing:
0420 print(f"PDF: {stage} missing: {', '.join(map(str, missing))}", file=sys.stderr)
0421 if files:
0422 output = report_dir / f"SummaryMipCalibration_{label}_{set_name}.pdf"
0423 unite(output, files)
0424 made[stage] = output
0425 print(f"PDF: wrote {output} from {len(files)} plots")
0426
0427 if "refine5" in made:
0428 final = report_dir / f"SummaryMipCalibration_Final_{set_name}.pdf"
0429 shutil.copyfile(made["refine5"], final)
0430 made["final"] = final
0431 print(f"PDF: wrote {final}")
0432
0433 for stage, label in (("pedestal", "Pedestal"), ("transfer", "Transfer")):
0434 root = plots / stage
0435 files = natural(root.rglob("*.pdf")) if root.is_dir() else []
0436 if files:
0437 output = report_dir / f"Summary{label}_{set_name}.pdf"
0438 unite(output, files)
0439 made[stage] = output
0440 print(f"PDF: wrote {output} from {len(files)} plots")
0441
0442 parts = [made[x] for x in ("pedestal", "transfer", "final") if x in made]
0443 if parts:
0444 output = report_dir / f"CalibrationPlotBook_{set_name}.pdf"
0445 unite(output, parts)
0446 print(f"PDF: wrote {output}")
0447
0448
0449 def main():
0450 a = args()
0451 try:
0452 work = find_work(a.work).resolve()
0453 out = (a.out or work / "report").resolve()
0454 reference = a.reference.resolve()
0455 out.mkdir(parents=True, exist_ok=True)
0456 compare(work, reference, out, a.set_name)
0457 if not a.no_pdf:
0458 make_pdfs(work, out, a.set_name)
0459 except (OSError, ValueError, RuntimeError, subprocess.CalledProcessError) as e:
0460 print(f"ERROR: {e}", file=sys.stderr)
0461 return 1
0462 print(f"Report directory: {out}")
0463 return 0
0464
0465
0466 if __name__ == "__main__":
0467 raise SystemExit(main())
