Make particle gun parameters configurable and add tracking scan

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

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)

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")