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