diff --git a/CHANGELOG.md b/CHANGELOG.md index 18872a2bed..79b7643f2d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,9 @@ it in future. * Add FileSummary to run_fixedTarget.py to save all the options for reference (#1140) * Add new 2026_04_01_SHiP_MainSpectrometerField_V13.root fieldmap +* Make particle gun polar angle and multiplicity configurable via `--thetaMin`, `--thetaMax`, `--nTracks` +* Add `run_tracking_scan.py` to sweep tracking benchmark over angle and multiplicity grids +* Add charge-ID efficiency metric and `--mixCharges` flag for particle/antiparticle generation ### Changed * Make artificial retina the baseline option for pattern recognition diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 6cc917a70e..6b9c92bfa6 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -70,6 +70,10 @@ pg_parser.add_argument( "--Dy", dest="Dy", type=float, help="size of the full uniform spread of PG ypos: (Vy - Dy/2, Vy + Dy/2)" ) +pg_parser.add_argument("--thetaMin", type=float, default=0, help="Minimum polar angle [deg] (default: 0)") +pg_parser.add_argument("--thetaMax", type=float, default=0, help="Maximum polar angle [deg] (default: 0)") +pg_parser.add_argument("--nTracks", type=int, default=1, help="Number of particles per event (default: 1)") +pg_parser.add_argument("--mixCharges", action="store_true", help="Generate equal mix of particle and antiparticle") # === End of PG commands === # === Genie subcommand === genie_parser = subparsers.add_parser("Genie", help="Genie for reading and processing neutrino interactions") @@ -517,23 +521,33 @@ # -----Particle Gun----------------------- if options.command == "PG": - myPgun = ROOT.FairBoxGenerator(options.pID, 1) - myPgun.SetPRange(options.Estart, options.Eend) - myPgun.SetPhiRange(0, 360) # // Azimuth angle range [degree] - myPgun.SetThetaRange(0, 0) # // Polar angle in lab system range [degree] - if options.multiplePG: - # multiple PG sources in the x-y plane; z is always the same! - myPgun.SetBoxXYZ( - (options.Vx - options.Dx / 2) * u.cm, - (options.Vy - options.Dy / 2) * u.cm, - (options.Vx + options.Dx / 2) * u.cm, - (options.Vy + options.Dy / 2) * u.cm, - options.Vz * u.cm, - ) + if options.mixCharges: + pdg_particle = ROOT.TDatabasePDG.Instance().GetParticle(abs(options.pID)) + if pdg_particle is None or pdg_particle.Charge() == 0: + print(f"WARNING: pID {options.pID} is neutral or unknown, disabling mixCharges") + options.mixCharges = False + if options.mixCharges and options.nTracks > 1: + pids = [(abs(options.pID), options.nTracks // 2), (-abs(options.pID), (options.nTracks + 1) // 2)] else: - # point source - myPgun.SetXYZ(options.Vx * u.cm, options.Vy * u.cm, options.Vz * u.cm) - primGen.AddGenerator(myPgun) + pids = [(options.pID, options.nTracks)] + for pid, mult in pids: + myPgun = ROOT.FairBoxGenerator(pid, mult) + myPgun.SetPRange(options.Estart, options.Eend) + myPgun.SetPhiRange(0, 360) # // Azimuth angle range [degree] + myPgun.SetThetaRange(options.thetaMin, options.thetaMax) # Polar angle in lab system [degree] + if options.multiplePG: + # multiple PG sources in the x-y plane; z is always the same! + myPgun.SetBoxXYZ( + (options.Vx - options.Dx / 2) * u.cm, + (options.Vy - options.Dy / 2) * u.cm, + (options.Vx + options.Dx / 2) * u.cm, + (options.Vy + options.Dy / 2) * u.cm, + options.Vz * u.cm, + ) + else: + # point source + myPgun.SetXYZ(options.Vx * u.cm, options.Vy * u.cm, options.Vz * u.cm) + primGen.AddGenerator(myPgun) # -----muon DIS Background------------------------ if options.mudis: ut.checkFileExists(inputFile) diff --git a/macro/run_tracking_benchmark.py b/macro/run_tracking_benchmark.py index 8f55ee3206..53fa1d3e52 100644 --- a/macro/run_tracking_benchmark.py +++ b/macro/run_tracking_benchmark.py @@ -32,6 +32,9 @@ parser.add_argument("--Dx", type=float, default=200.0, help="Position spread in x [cm] (default: 200)") parser.add_argument("--Dy", type=float, default=300.0, help="Position spread in y [cm] (default: 300)") parser.add_argument("--nTracks", type=int, default=1, help="Tracks per event (default: 1)") +parser.add_argument("--thetaMin", type=float, default=0, help="Minimum polar angle [deg] (default: 0)") +parser.add_argument("--thetaMax", type=float, default=0, help="Maximum polar angle [deg] (default: 0)") +parser.add_argument("--mixCharges", action="store_true", help="Generate equal mix of particle and antiparticle") parser.add_argument("--tag", default="benchmark", help="Output file tag (default: benchmark)") parser.add_argument("--output-json", default=None, help="JSON metrics output path") parser.add_argument("--seed", type=int, default=42, help="Random seed (default: 42)") @@ -48,23 +51,28 @@ if not os.path.exists(options.outputDir): os.makedirs(options.outputDir) +options.outputDir = os.path.abspath(options.outputDir) tag = options.tag -sim_file = f"{options.outputDir}/sim_{tag}.root" -geo_file = f"{options.outputDir}/geo_{tag}.root" -reco_file = f"{options.outputDir}/sim_{tag}_rec.root" -json_file = options.output_json or f"{options.outputDir}/tracking_metrics.json" -histo_file = f"{options.outputDir}/tracking_benchmark_histos.root" +sim_file = os.path.join(options.outputDir, f"sim_{tag}.root") +geo_file = os.path.join(options.outputDir, f"geo_{tag}.root") +reco_file = os.path.join(options.outputDir, f"sim_{tag}_rec.root") +json_file = ( + os.path.abspath(options.output_json) + if options.output_json + else os.path.join(options.outputDir, "tracking_metrics.json") +) +histo_file = os.path.join(options.outputDir, "tracking_benchmark_histos.root") fairship = os.environ.get("FAIRSHIP", "") -def run_phase(description: str, cmd: list[str]) -> None: +def run_phase(description: str, cmd: list[str], **kwargs: str) -> None: """Run a subprocess phase, raising on failure.""" print("=" * 60) print(f"{description}") print("=" * 60) - result = subprocess.run(cmd, check=False) + result = subprocess.run(cmd, check=False, **kwargs) if result.returncode != 0: print(f"FAILED: {description} (exit code {result.returncode})") sys.exit(result.returncode) @@ -107,7 +115,15 @@ def run_phase(description: str, cmd: list[str]) -> None: str(options.Dx), "--Dy", str(options.Dy), + "--nTracks", + str(options.nTracks), + "--thetaMin", + str(options.thetaMin), + "--thetaMax", + str(options.thetaMax), ] +if options.mixCharges: + sim_cmd.append("--mixCharges") run_phase("Phase 1: Simulation", sim_cmd) @@ -134,7 +150,7 @@ def run_phase(description: str, cmd: list[str]) -> None: if options.debug > 0: reco_cmd.append("--Debug") -run_phase("Phase 2: Reconstruction", reco_cmd) +run_phase("Phase 2: Reconstruction", reco_cmd, cwd=options.outputDir) if not os.path.exists(reco_file): print(f"ERROR: Reconstruction output {reco_file} not found") diff --git a/macro/run_tracking_scan.py b/macro/run_tracking_scan.py new file mode 100644 index 0000000000..4674f7c054 --- /dev/null +++ b/macro/run_tracking_scan.py @@ -0,0 +1,387 @@ +#!/usr/bin/env python3 +# SPDX-License-Identifier: LGPL-3.0-or-later +# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration + +"""Scan tracking benchmark over a grid of polar angles and multiplicities. + +Calls run_tracking_benchmark.py for each (theta, nTracks) point, collects +metrics into a combined JSON, and produces summary plots with matplotlib. + +Example usage: + # 1D theta scan at single-track multiplicity + python macro/run_tracking_scan.py -n 200 --scan-mode theta + + # 1D multiplicity scan at theta=0 + python macro/run_tracking_scan.py -n 200 --scan-mode nTracks + + # Full 2D grid + python macro/run_tracking_scan.py -n 500 --scan-mode grid + + # Custom theta values, keep intermediate files + python macro/run_tracking_scan.py -n 200 --thetas 0,5,15,30 --keep-files +""" + +from __future__ import annotations + +import glob +import json +import os +import subprocess +import sys +from argparse import ArgumentParser +from collections import defaultdict + +parser = ArgumentParser(description="Scan tracking benchmark over angle and multiplicity") +parser.add_argument("-n", "--nEvents", type=int, default=500, help="Events per scan point (default: 500)") +parser.add_argument( + "--thetas", + default="0,1,2,3,5,7,10,15,20,30", + help="Comma-separated fixed theta values in degrees (default: 0,1,2,3,5,7,10,15,20,30)", +) +parser.add_argument( + "--nTracksList", + default="1,2,5,10,20", + help="Comma-separated multiplicities (default: 1,2,5,10,20)", +) +parser.add_argument( + "--scan-mode", + choices=["theta", "nTracks", "grid"], + default="grid", + help="Scan mode: theta (fixed theta per point), nTracks (vary multiplicity at thetaRange), grid (full 2D)", +) +parser.add_argument( + "--thetaRange", + default=None, + help="Theta range 'min,max' in degrees for nTracks scan mode (default: 0,10). " + "Tracks are generated uniformly in this range. Ignored in theta scan mode.", +) +parser.add_argument("--pID", type=int, default=13, help="Particle PDG ID (default: 13, mu-)") +parser.add_argument("--Estart", type=float, default=1.0, help="Start of energy range in GeV (default: 1)") +parser.add_argument("--Eend", type=float, default=100.0, help="End of energy range in GeV (default: 100)") +parser.add_argument("--seed", type=int, default=42, help="Base random seed (default: 42)") +parser.add_argument("-o", "--outputDir", default="./scan_results", help="Base output directory") +parser.add_argument("--mixCharges", action="store_true", help="Generate equal mix of particle and antiparticle") +parser.add_argument("--keep-files", action="store_true", help="Keep intermediate sim/reco ROOT files") +parser.add_argument( + "--debug", + type=int, + default=0, + choices=range(0, 4), + help="FairLogger verbosity: 0=info, 1=debug, 2=debug1, 3=debug2", +) + +options = parser.parse_args() + +thetas = [float(x) for x in options.thetas.split(",")] +n_tracks_list = [int(x) for x in options.nTracksList.split(",")] + +# Parse theta range for nTracks mode +if options.thetaRange: + tr_min, tr_max = (float(x) for x in options.thetaRange.split(",")) +else: + tr_min, tr_max = 0.0, 10.0 + +# Build scan grid: each point is (thetaMin, thetaMax, nTracks, label_theta) +# In theta mode: fixed theta per point (thetaMin == thetaMax) +# In nTracks mode: uniform theta range, vary multiplicity +# In grid mode: fixed theta per point, vary both +if options.scan_mode == "theta": + scan_points = [(t, t, 1, t) for t in thetas] +elif options.scan_mode == "nTracks": + scan_points = [(tr_min, tr_max, n, f"{tr_min}-{tr_max}") for n in n_tracks_list] +else: # grid + scan_points = [(t, t, n, t) for t in thetas for n in n_tracks_list] + +os.makedirs(options.outputDir, exist_ok=True) + +# Dx=200, Dy=300 cm (benchmark defaults) -> area in m^2 +GUN_AREA_M2 = 200 * 300 * 1e-4 + +fairship = os.environ.get("FAIRSHIP", "") +benchmark_script = ( + os.path.join(fairship, "macro", "run_tracking_benchmark.py") if fairship else "macro/run_tracking_benchmark.py" +) + +all_results: list[dict] = [] +n_total = len(scan_points) + +for i, (theta_min, theta_max, n_tracks, label_theta) in enumerate(scan_points): + tag = f"scan_theta{label_theta}_nTr{n_tracks}" + point_dir = os.path.join(options.outputDir, tag) + json_path = os.path.join(point_dir, "tracking_metrics.json") + histo_path = os.path.join(point_dir, "tracking_benchmark_histos.root") + + theta_desc = f"{theta_min}°" if theta_min == theta_max else f"[{theta_min},{theta_max}]°" + print(f"\n{'=' * 60}") + print(f"Scan point {i + 1}/{n_total}: theta={theta_desc}, nTracks={n_tracks}") + print(f"{'=' * 60}") + + cmd = [ + sys.executable, + benchmark_script, + "-n", + str(options.nEvents), + "--thetaMin", + str(theta_min), + "--thetaMax", + str(theta_max), + "--nTracks", + str(n_tracks), + "--tag", + tag, + "--seed", + str(options.seed + i), + "-o", + point_dir, + "--output-json", + json_path, + "--pID", + str(options.pID), + "--Estart", + str(options.Estart), + "--Eend", + str(options.Eend), + "--debug", + str(options.debug), + ] + if options.mixCharges: + cmd.append("--mixCharges") + + result = subprocess.run(cmd, check=False) + if result.returncode != 0: + print(f"WARNING: scan point theta={theta_desc}, nTracks={n_tracks} failed (exit code {result.returncode})") + continue + + try: + with open(json_path) as f: + metrics = json.load(f) + benchmark_metrics = metrics["tracking_benchmark"] + except (FileNotFoundError, json.JSONDecodeError, KeyError) as e: + print(f"WARNING: failed to read results from {json_path}: {e}") + continue + + all_results.append( + { + "theta": label_theta, + "thetaMin": theta_min, + "thetaMax": theta_max, + "nTracks": n_tracks, + "density_per_m2": round(n_tracks / GUN_AREA_M2, 4), + "metrics": benchmark_metrics, + } + ) + + # Clean up intermediate ROOT files + if not options.keep_files: + for pattern in [f"sim_{tag}.root", f"geo_{tag}.root", f"params_{tag}.root", f"sim_{tag}_rec*.root"]: + for fpath in glob.glob(os.path.join(point_dir, pattern)): + os.remove(fpath) + print(f" Cleaned up {fpath}") + +# Save combined results +combined = { + "scan_config": { + "scan_mode": options.scan_mode, + "thetas": thetas, + "nTracksList": n_tracks_list, + "nEvents": options.nEvents, + "pID": options.pID, + "Estart": options.Estart, + "Eend": options.Eend, + "seed": options.seed, + }, + "scan_results": all_results, +} + +combined_json = os.path.join(options.outputDir, "scan_results.json") +with open(combined_json, "w") as f: + json.dump(combined, f, indent=2) +print(f"\nCombined results saved to {combined_json}") + +# Save TGraphErrors to ROOT file +import ROOT + +ROOT.gROOT.SetBatch(True) + +metric_names = [ + "efficiency", + "ghost_rate", + "clone_rate", + "dp_over_p_sigma", + "n_reconstructible", + "charge_id_efficiency", +] +metric_titles = { + "efficiency": "Tracking efficiency", + "ghost_rate": "Ghost rate", + "clone_rate": "Clone rate", + "dp_over_p_sigma": "#Deltap/p resolution (#sigma)", + "n_reconstructible": "Reconstructible tracks", + "charge_id_efficiency": "Charge-ID efficiency", +} + +# Group results — only include numeric theta values for "vs theta" graphs +by_nTracks: dict[int, list[tuple[float, dict]]] = defaultdict(list) +by_theta: dict[str, list[tuple[int, dict]]] = defaultdict(list) + +for r in all_results: + theta_label = r["theta"] + if isinstance(theta_label, (int, float)): + by_nTracks[r["nTracks"]].append((float(theta_label), r["metrics"])) + by_theta[str(theta_label)].append((r["nTracks"], r["metrics"])) + +summary_root = os.path.join(options.outputDir, "scan_summary.root") +f_out = ROOT.TFile.Open(summary_root, "recreate") + +for metric in metric_names: + for n_tr, points in sorted(by_nTracks.items()): + points.sort() + valid = [(x, m) for x, m in points if m[metric]["value"] is not None] + if not valid: + continue + gr = ROOT.TGraphErrors(len(valid)) + gr.SetName(f"g_{metric}_vs_theta_nTr{n_tr}") + gr.SetTitle(f"{metric_titles[metric]}, nTracks={n_tr};#theta [deg];{metric_titles[metric]}") + for j, (theta_val, m) in enumerate(valid): + gr.SetPoint(j, theta_val, m[metric]["value"]) + gr.SetPointError(j, 0, m[metric].get("uncertainty") or 0) + gr.Write() + + for theta_label, points in sorted(by_theta.items()): + points.sort() + valid = [(x, m) for x, m in points if m[metric]["value"] is not None] + if not valid: + continue + gr = ROOT.TGraphErrors(len(valid)) + gr.SetName(f"g_{metric}_vs_nTracks_theta{theta_label}") + gr.SetTitle(f"{metric_titles[metric]}, #theta={theta_label}#circ;nTracks;{metric_titles[metric]}") + for j, (n_tr_val, m) in enumerate(valid): + gr.SetPoint(j, n_tr_val, m[metric]["value"]) + gr.SetPointError(j, 0, m[metric].get("uncertainty") or 0) + gr.Write() + +f_out.Close() +print(f"ROOT summary saved to {summary_root}") + +# Generate matplotlib plots +try: + import matplotlib.pyplot as plt + import numpy as np + + # Axis limits: rates are bounded to [0, max], efficiency to [min, 1.05] + metric_ylimits: dict[str, tuple[float | None, float | None]] = { + "efficiency": (None, 1.05), + "ghost_rate": (0, None), + "clone_rate": (0, None), + "dp_over_p_sigma": (0, None), + "n_reconstructible": (0, None), + "charge_id_efficiency": (None, 1.05), + } + + metric_ylabels = { + "efficiency": "Efficiency", + "ghost_rate": "Ghost rate", + "clone_rate": "Clone rate", + "dp_over_p_sigma": r"$\Delta p / p$ resolution ($\sigma$)", + "n_reconstructible": f"Reconstructible tracks (out of {options.nEvents})", + "charge_id_efficiency": "Charge-ID efficiency", + } + + def _save_fig(fig: plt.Figure, base_path: str) -> None: + """Save figure as both PDF and PNG.""" + fig.savefig(base_path + ".pdf") + fig.savefig(base_path + ".png", dpi=150) + + def _get_values(points: list, metric: str) -> tuple[list, list, list]: + """Extract x, y, yerr from sorted scan points, skipping None values.""" + x, y, yerr = [], [], [] + for p in points: + val = p[1][metric]["value"] + if val is None: + continue + x.append(p[0]) + y.append(val) + yerr.append(p[1][metric].get("uncertainty") or 0) + return x, y, yerr + + n_curves_theta = max(len(by_nTracks), 1) + n_curves_mult = max(len(by_theta), 1) + colours_theta = plt.cm.viridis(np.linspace(0, 0.9, n_curves_theta)) + colours_mult = plt.cm.plasma(np.linspace(0, 0.9, n_curves_mult)) + + for metric in metric_names: + ymin, ymax = metric_ylimits.get(metric, (None, None)) + ylabel = metric_ylabels.get(metric, metric) + + # Plots vs theta + if len(by_nTracks) > 0 and any(len(pts) > 1 for pts in by_nTracks.values()): + fig, ax = plt.subplots(figsize=(8, 5)) + for idx, (n_tr, points) in enumerate(sorted(by_nTracks.items())): + points.sort() + x, y, yerr = _get_values(points, metric) + if not x: + continue + ax.errorbar( + x, + y, + yerr=yerr, + marker="o", + linestyle="none", + capsize=3, + label=f"nTracks={n_tr}", + color=colours_theta[idx], + ) + ax.set_xlabel(r"$\theta$ [deg]") + ax.set_ylabel(ylabel) + ax.set_title(f"{ylabel} vs polar angle") + if ymin is not None: + ax.set_ylim(bottom=ymin) + if ymax is not None: + ax.set_ylim(top=ymax) + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + _save_fig(fig, os.path.join(options.outputDir, f"scan_{metric}_vs_theta")) + plt.close(fig) + + # Plots vs nTracks + if len(by_theta) > 0 and any(len(pts) > 1 for pts in by_theta.values()): + fig, ax = plt.subplots(figsize=(8, 5)) + for idx, (theta_label, points) in enumerate(sorted(by_theta.items())): + points.sort() + x, y, yerr = _get_values(points, metric) + if not x: + continue + try: + theta_legend = rf"$\theta$={float(theta_label):.0f}°" + except ValueError: + theta_legend = rf"$\theta \in [{theta_label}]$°" + ax.errorbar( + x, + y, + yerr=yerr, + marker="s", + linestyle="none", + capsize=3, + label=theta_legend, + color=colours_mult[idx], + ) + ax.set_xlabel("Tracks per event") + ax.set_ylabel(ylabel) + ax.set_title(f"{ylabel} vs multiplicity (gun area = {GUN_AREA_M2:.1f} m²)") + if ymin is not None: + ax.set_ylim(bottom=ymin) + if ymax is not None: + ax.set_ylim(top=ymax) + ax.legend() + ax.grid(True, alpha=0.3) + fig.tight_layout() + _save_fig(fig, os.path.join(options.outputDir, f"scan_{metric}_vs_nTracks")) + plt.close(fig) + + print(f"Plots saved to {options.outputDir}/scan_*.{{pdf,png}}") + +except ImportError: + print("matplotlib not available — skipping plots (ROOT summary file still produced)") + +print(f"\nScan complete: {len(all_results)}/{n_total} points succeeded.") diff --git a/python/tracking_benchmark.py b/python/tracking_benchmark.py index 76bdbe7bf0..bded6d8725 100644 --- a/python/tracking_benchmark.py +++ b/python/tracking_benchmark.py @@ -12,6 +12,7 @@ from __future__ import annotations import json +import logging import math from typing import Any @@ -19,8 +20,10 @@ ROOT.gROOT.SetBatch(True) +logger = logging.getLogger(__name__) -def wilson_interval(k: int, n: int) -> float: + +def wilson_interval(k: int, n: int) -> float | None: """Wilson score interval half-width for a binomial proportion. Parameters @@ -32,11 +35,12 @@ def wilson_interval(k: int, n: int) -> float: Returns ------- - float - Half-width of the 68% Wilson score interval (~1 sigma). + float or None + Half-width of the 68% Wilson score interval (~1 sigma), + or None if n == 0 (undefined proportion). """ if n == 0: - return 0.0 + return None z = 1.0 # 1-sigma p = k / n denom = 1 + z**2 / n @@ -44,6 +48,11 @@ def wilson_interval(k: int, n: int) -> float: return spread +def _round_or_none(value: float | None, ndigits: int = 6) -> float | None: + """Round a value, passing through None unchanged.""" + return round(value, ndigits) if value is not None else None + + class TrackingBenchmark: """Compute tracking benchmark metrics from simulation and reconstruction files. @@ -199,6 +208,8 @@ def compute_metrics(self) -> dict[str, Any]: h_chi2ndf = ROOT.TH1D("h_chi2ndf", "#chi^{2}/ndf;#chi^{2}/ndf;Entries", 100, 0, 20) h_p_truth = ROOT.TH1D("h_p_truth", "p_{truth} (reconstructible);p [GeV/c];Entries", 50, 0, 120) h_p_matched = ROOT.TH1D("h_p_matched", "p_{truth} (matched);p [GeV/c];Entries", 50, 0, 120) + h_p_correct_charge = ROOT.TH1D("h_p_correct_charge", "p_{truth} (correct charge);p [GeV/c];Entries", 50, 0, 120) + h_p_wrong_charge = ROOT.TH1D("h_p_wrong_charge", "p_{truth} (wrong charge);p [GeV/c];Entries", 50, 0, 120) # Counters n_reconstructible = 0 @@ -206,6 +217,8 @@ def compute_metrics(self) -> dict[str, Any]: n_total_reco = 0 n_matched_reco = 0 # reco tracks passing purity cut n_clone_reco = 0 # extra matches beyond the first for same MC particle + n_correct_charge = 0 + n_wrong_charge = 0 for i_event in range(n_events): self.sim_tree.GetEvent(i_event) @@ -287,31 +300,57 @@ def compute_metrics(self) -> dict[str, Any]: h_dty.Fill(ty_reco - ty_t) h_p_matched.Fill(p_truth) + + # Charge identification + q_reco = fitted_state.getCharge() + mc_pdg = self.sim_tree.MCTrack[mc_id].GetPdgCode() + mc_particle = self.PDG.GetParticle(mc_pdg) + if mc_particle: + q_truth = mc_particle.Charge() / 3.0 + if q_reco * q_truth > 0: + n_correct_charge += 1 + h_p_correct_charge.Fill(p_truth) + else: + n_wrong_charge += 1 + h_p_wrong_charge.Fill(p_truth) except Exception: - pass + logger.debug( + "Failed to extract fitted state for mc_id=%d, p_truth=%.3f", + mc_id, + p_truth, + exc_info=True, + ) n_matched_mc += len(matched_mc_this_event) - # Compute metrics + # Compute metrics (None when denominator is zero) n_ghost_reco = n_total_reco - n_matched_reco - efficiency = n_matched_mc / n_reconstructible if n_reconstructible > 0 else 0.0 + efficiency = n_matched_mc / n_reconstructible if n_reconstructible > 0 else None efficiency_unc = wilson_interval(n_matched_mc, n_reconstructible) - clone_rate = n_clone_reco / n_matched_reco if n_matched_reco > 0 else 0.0 + clone_rate = n_clone_reco / n_matched_reco if n_matched_reco > 0 else None clone_rate_unc = wilson_interval(n_clone_reco, n_matched_reco) - ghost_rate = n_ghost_reco / n_total_reco if n_total_reco > 0 else 0.0 + ghost_rate = n_ghost_reco / n_total_reco if n_total_reco > 0 else None ghost_rate_unc = wilson_interval(n_ghost_reco, n_total_reco) # Fit dp/p with Gaussian - dp_p_sigma = h_dp_over_p.GetRMS() - dp_p_sigma_unc = h_dp_over_p.GetRMSError() - if h_dp_over_p.GetEntries() > 20: - fit_result = h_dp_over_p.Fit("gaus", "QS") - if fit_result and int(fit_result) == 0: - dp_p_sigma = fit_result.Parameter(2) - dp_p_sigma_unc = fit_result.ParError(2) + dp_p_sigma: float | None = None + dp_p_sigma_unc: float | None = None + if h_dp_over_p.GetEntries() > 0: + dp_p_sigma = h_dp_over_p.GetRMS() + dp_p_sigma_unc = h_dp_over_p.GetRMSError() + if h_dp_over_p.GetEntries() > 20: + fit_result = h_dp_over_p.Fit("gaus", "QS") + if fit_result and int(fit_result) == 0: + dp_p_sigma = fit_result.Parameter(2) + dp_p_sigma_unc = fit_result.ParError(2) + + # Charge-ID efficiency + n_charge_total = n_correct_charge + n_wrong_charge + charge_id_eff = n_correct_charge / n_charge_total if n_charge_total > 0 else None + charge_id_eff_unc = wilson_interval(n_correct_charge, n_charge_total) self.metrics = { "tracking_benchmark": { @@ -319,23 +358,23 @@ def compute_metrics(self) -> dict[str, Any]: "n_reconstructible": {"value": int(n_reconstructible), "compare": "exact"}, "n_total_reco": {"value": int(n_total_reco), "compare": "exact"}, "efficiency": { - "value": round(efficiency, 6), - "uncertainty": round(efficiency_unc, 6), + "value": _round_or_none(efficiency), + "uncertainty": _round_or_none(efficiency_unc), "compare": "statistical", }, "clone_rate": { - "value": round(clone_rate, 6), - "uncertainty": round(clone_rate_unc, 6), + "value": _round_or_none(clone_rate), + "uncertainty": _round_or_none(clone_rate_unc), "compare": "statistical", }, "ghost_rate": { - "value": round(ghost_rate, 6), - "uncertainty": round(ghost_rate_unc, 6), + "value": _round_or_none(ghost_rate), + "uncertainty": _round_or_none(ghost_rate_unc), "compare": "statistical", }, "dp_over_p_sigma": { - "value": round(dp_p_sigma, 6), - "uncertainty": round(dp_p_sigma_unc, 6), + "value": _round_or_none(dp_p_sigma), + "uncertainty": _round_or_none(dp_p_sigma_unc), "compare": "statistical", }, "dx_rms": { @@ -358,6 +397,13 @@ def compute_metrics(self) -> dict[str, Any]: "uncertainty": round(h_dty.GetRMSError(), 6), "compare": "statistical", }, + "charge_id_efficiency": { + "value": _round_or_none(charge_id_eff), + "uncertainty": _round_or_none(charge_id_eff_unc), + "compare": "statistical", + }, + "n_correct_charge": {"value": int(n_correct_charge), "compare": "exact"}, + "n_wrong_charge": {"value": int(n_wrong_charge), "compare": "exact"}, } } @@ -371,6 +417,8 @@ def compute_metrics(self) -> dict[str, Any]: "h_chi2ndf": h_chi2ndf, "h_p_truth": h_p_truth, "h_p_matched": h_p_matched, + "h_p_correct_charge": h_p_correct_charge, + "h_p_wrong_charge": h_p_wrong_charge, } return self.metrics @@ -389,24 +437,38 @@ def save_histograms(self, output_path: str) -> None: f_out.Close() print(f"Histograms saved to {output_path}") + @staticmethod + def _fmt(value: float | None, uncertainty: float | None, fmt: str = ".4f") -> str: + """Format value +/- uncertainty, handling None.""" + if value is None: + return "undefined (no data)" + if uncertainty is None: + return f"{value:{fmt}}" + return f"{value:{fmt}} +/- {uncertainty:{fmt}}" + def print_summary(self) -> None: """Print a human-readable summary of the benchmark results.""" if not self.metrics: print("No metrics computed yet. Call compute_metrics() first.") return m = self.metrics["tracking_benchmark"] + fmt = self._fmt print("\n=== Tracking Benchmark Summary ===") print(f" Events: {m['n_events']['value']}") print(f" Reconstructible: {m['n_reconstructible']['value']}") print(f" Total reco: {m['n_total_reco']['value']}") - print(f" Efficiency: {m['efficiency']['value']:.4f} +/- {m['efficiency']['uncertainty']:.4f}") - print(f" Clone rate: {m['clone_rate']['value']:.4f} +/- {m['clone_rate']['uncertainty']:.4f}") - print(f" Ghost rate: {m['ghost_rate']['value']:.4f} +/- {m['ghost_rate']['uncertainty']:.4f}") - print(f" dp/p sigma: {m['dp_over_p_sigma']['value']:.6f} +/- {m['dp_over_p_sigma']['uncertainty']:.6f}") - print(f" dx RMS: {m['dx_rms']['value']:.4f} +/- {m['dx_rms']['uncertainty']:.4f} cm") - print(f" dy RMS: {m['dy_rms']['value']:.4f} +/- {m['dy_rms']['uncertainty']:.4f} cm") - print(f" dtx RMS: {m['dtx_rms']['value']:.6f} +/- {m['dtx_rms']['uncertainty']:.6f}") - print(f" dty RMS: {m['dty_rms']['value']:.6f} +/- {m['dty_rms']['uncertainty']:.6f}") + print(f" Efficiency: {fmt(m['efficiency']['value'], m['efficiency']['uncertainty'])}") + print(f" Clone rate: {fmt(m['clone_rate']['value'], m['clone_rate']['uncertainty'])}") + print(f" Ghost rate: {fmt(m['ghost_rate']['value'], m['ghost_rate']['uncertainty'])}") + print(f" dp/p sigma: {fmt(m['dp_over_p_sigma']['value'], m['dp_over_p_sigma']['uncertainty'], '.6f')}") + print(f" dx RMS: {fmt(m['dx_rms']['value'], m['dx_rms']['uncertainty'])} cm") + print(f" dy RMS: {fmt(m['dy_rms']['value'], m['dy_rms']['uncertainty'])} cm") + print(f" dtx RMS: {fmt(m['dtx_rms']['value'], m['dtx_rms']['uncertainty'], '.6f')}") + print(f" dty RMS: {fmt(m['dty_rms']['value'], m['dty_rms']['uncertainty'], '.6f')}") + print( + f" Charge-ID eff: {fmt(m['charge_id_efficiency']['value'], m['charge_id_efficiency']['uncertainty'])}" + ) + print(f" Correct/wrong: {m['n_correct_charge']['value']}/{m['n_wrong_charge']['value']}") print("==================================\n") def __del__(self) -> None: