modified: TrackRecon.C implemented and applied beamtilt parameter calcualtion and application

modified:   run_27Al.sh
	new file:   scratch/FitBeamAxis.C code to figure out beamtilt

CLAUDE.md

@@ -1,21 +1,11 @@
 # ANASEN Analysis — working notes for Claude
 
 ## Agent behavior
-- **Architecture & Remotes:** 
-  - `upstream` points to the canonical public repo (read-only reference).
-  - `origin` points to your private sandbox repo (read/write scratchpad).
-  - **NEVER attempt to push to `upstream`.** All upstream synchronization is strictly one-way (`git fetch upstream` -> merge).
-- **Upstream Sync Routine:**
-  - Before applying requested features or edits, sync with `upstream` if requested.
-  - Never pull into a dirty working tree. Use: `git stash` -> `git fetch upstream` -> `git merge upstream/devel_vignesh` -> `git stash pop`.
-  - If a merge conflict occurs in analysis macros (`TrackRecon.C`, `MakeVertex.C`, kinematics scripts, or calibration tables), **STOP immediately**. Do not attempt to guess physics logic; flag the exact conflict lines and await instruction.
-- **Commits & Pushes to `origin`:**
-  - Commit only when a coherent, discrete change is finished or when explicitly asked. Write concise, conventional commit messages.
-  - Push **only** to `origin` (`git push origin <branch>`). Push after completing a requested unit of work so the sandbox stays backed up.
 - ALWAYS output a brief, step-by-step plan before modifying files or running commands.
+- Commit only when explicitly asked. Do not push to remote on your own initiative; pushing is allowed only when the session's instructions name a branch to push to (e.g. a `claude/*` working branch), and never to `devel_vignesh` or any other shared branch.
 - **Targeted edits only.** Never regenerate a whole file to make a local change.
-- **Never read `MakeVertex.C`**
+- **Never read `MakeVertex.C`** (4,003 lines, 240 KB). Same rule as `TrackRecon.C` below if you must touch it: grep, then Read with offset/limit.
-- **Never read `TrackRecon.C` end to end** (~4,250 lines; ~63k tokens each, and it stays in the prefix for the rest of the session). Grep for the symbol, then Read with offset/limit around the hit. If a question genuinely needs the whole file, say so first. Do not explore the codebase open-endedly.
+- **Never read `TrackRecon.C` end to end** (4,268 lines, 218 KB; ~60k tokens, and it stays in the prefix for the rest of the session). Grep for the symbol, then Read with offset/limit around the hit. If a question genuinely needs the whole file, say so first. Do not explore the codebase open-endedly.
 - **Do not spawn subagents** for work that can be done inline — each starts cold and re-derives context already loaded. Worth it only for wide parallel searches where the conclusion is all that's needed.
 - **You cannot verify your own work here.** No test suite, no CI, no data in the repo; validation is visual, via ROOT `TBrowser`, on a machine that has the data. Say explicitly when a change is unverified — never report a physics change as working because it compiled. The dangerous bugs are silent: a flipped sign in a kinematics term, an inverted spline direction, a wrong argument order, an off-by-one in a wire mask. None of them throw. The run completes, the plot renders, and the spectrum is wrong.
 - **Review once per coherent change, not per edit.** Scale depth to blast radius: scratch macros, shell, and docs are shallow; `TrackRecon.C`, `MakeVertex.C`, `Armory/`, or any calibration fit is deep and should tolerate uncertain findings — silent failure is what's being hunted. Run `/security-review` only for file I/O or shell interpolation in the `run_*.sh` drivers; it is not a physics tool.
@@ -34,8 +24,7 @@ All config is `getenv`-driven; the run scripts `export` it. Defaults are in `Tra
 Physics / results-affecting:
 `DATASET` `reactiondata` `CO2percent` `pressure_in_torr` `CATHODE_GAIN` `PC_ENERGY_CALIBRATION`
 `source_vertex` `CUTLIST` `DITHER_SIGMA` `RNG_SEED` `BEAM_AXIS_X` `BEAM_AXIS_Y`
-`timecut_low` `timecut_high` `DISABLE_BAD_ANODE_WIRES` `A1C1_LOWBAND_RFACTOR` `A1C1_CFRAC_SPLIT`
+`timecut_low` `timecut_high` `DISABLE_BAD_ANODE_WIRES` `A1C1_LOWBAND_RFACTOR`
-`A1C1_ANODEE_COFF` `A1C1_ANODEE_REF` `A1C1_MISSING_FMAX`
 `A1C1_Z_SCALE_QQQ` `A1C1_Z_SCALE_SX3` `A1C1_Z_OFF_QQQ` `A1C1_Z_OFF_SX3`
 
 Plumbing: `OUT_DIR` `RUN_NUMBER` `FLUSH_BARRIER` `MAX_RSS_MB` `MEMCHECK_STRIDE`
@@ -44,9 +33,9 @@ Plumbing: `OUT_DIR` `RUN_NUMBER` `FLUSH_BARRIER` `MAX_RSS_MB` `MEMCHECK_STRIDE`
 It multiplies `catima.dedx` when the lookup tables are generated. Changing it invalidates every `.dat` table — re-run `Eloss.py` (the run scripts do this automatically before processing). Tables are keyed `eloss_calculations/<species>_lookup_<E>MeV_<P>torr_<CO2>pc.dat` (e.g. `alpha_lookup_50MeV_250torr_3pc.dat`).
 
 ## Footguns
-- **`HistPlotter` keys `oMap` by histogram NAME only.** The folder argument does not namespace. Two `Fill` calls with the same name in different folders merge silently — no warning, no error, wrong plot. Always make names unique unless required.
+- **`HistPlotter` keys `oMap` by histogram NAME only.** The folder argument does not namespace. Two `Fill` calls with the same name and the *same* folder merge silently — no warning, no error, wrong plot. A same-name/different-folder (or toplevel-vs-folder) clash does print to `stderr` and then merges anyway, but the run scripts send stdout+stderr to `/dev/null`, so in a batch run it is invisible either way. Always make names unique.
-- **Si–PC coincidence sits at negative dt.** The real band is roughly `-450 < t_si - t_pc < -200`, not near zero. Do not "correct" a timing gate toward zero.
+- **Si–PC coincidence sits at negative dt.** The band is centred well below zero, not near it. The one gate of record is `siPcCoincident()` in `TrackRecon.C` (`kSiPcDtMin = -500`, `kSiPcDtMax = 150` ns); anything else is a plotting range, not a cut. Do not "correct" a timing gate toward zero.
-- **RNG**: one shared `anasenRandom` (`TRandom3`, seed 4357, `RNG_SEED` to override) — do not introduce fresh `TRandom3(0)` instances; it breaks reproducibility.
+- **RNG**: in `TrackRecon.C`, one shared `anasenRandom` (`TRandom3`, seed 4357, `RNG_SEED` to override), aliased by reference wherever a local `rnd`/`rand` appears — do not introduce fresh `TRandom3(0)` instances; it breaks reproducibility. `MakeVertex.C` predates this and still has ~12 independent default-seeded `TRandom3` objects; do not assume its output is reproducible against `TrackRecon.C`.
 - **Pass By Reference in Event Loops:** Always iterate via `const auto&` over cluster/hit vectors to prevent expensive struct copies in multi-million event runs.
 
 ## Invariants to check on any physics change

TrackRecon.C

@@ -42,8 +42,8 @@ Int_t colors[40] = {
 bool process_alpha_proton_scattering = false,
      doMiscHistograms = true,
      doRawHistos = false,
-     doPCSX3ClusterAnalysis = false,
+     doPCSX3ClusterAnalysis = true,
-     doPCQQQClusterAnalysis = false,
+     doPCQQQClusterAnalysis = true,
      doOldAnalysis = false,
      BenchMark = false,
      onewire_analysis = true,
@@ -75,6 +75,9 @@ double source_vertex = 53.0,
        a1c1_z_off_sx3 = 2.52614,
        beam_axis_x = 0.0,
        beam_axis_y = 0.0,
+       beam_axis_z0 = 0.0,   // reference z at which (beam_axis_x, beam_axis_y) is specified
+       beam_tilt_x = 0.0,    // dx/dz slope of the beam axis
+       beam_tilt_y = 0.0,    // dy/dz slope of the beam axis
        ta_foil_z_mm = 0.0,
        alpha_source_mev = 5.486;
 
@@ -122,23 +125,111 @@ inline SiPcPid classifyByAnodeDe(double anodeE_MeV)
   return (anodeE_MeV < kProtonAlphaAnodeDeGate_MeV) ? SiPcPid::kProton : SiPcPid::kAlpha;
 }
 
-inline TVector3 beamVertex(const TVector3 &si, const TVector3 &dir)
+inline TVector3 beamAxisOrigin()
 {
-  double d = dir.X() * dir.X() + dir.Y() * dir.Y();
+  return TVector3(beam_axis_x, beam_axis_y, beam_axis_z0);
-  double t = (d > 0.0) ? -((si.X() - beam_axis_x) * dir.X() + (si.Y() - beam_axis_y) * dir.Y()) / d : 0.0;
-  return si + t * dir;
 }
-inline double beamPerp(const TVector3 &p)
+inline TVector3 beamAxisDirUnit()
 {
-  return TMath::Sqrt((p.X() - beam_axis_x) * (p.X() - beam_axis_x) + (p.Y() - beam_axis_y) * (p.Y() - beam_axis_y));
+  return TVector3(beam_tilt_x, beam_tilt_y, 1.0).Unit();
 }
-// A point on the beam axis at height z. Every theta/phi reference point used to be
+
-// spelled beamAxisPoint(z), which silently ignored BEAM_AXIS_X/Y even though
-// Begin() prints them as configured parameters and pcEnergyCalibrationAccumulate
-// already built its source_pos the correct way.
 inline TVector3 beamAxisPoint(double z)
 {
-  return TVector3(beam_axis_x, beam_axis_y, z);
+  double dz = z - beam_axis_z0;
+  return TVector3(beam_axis_x + beam_tilt_x * dz,
+                  beam_axis_y + beam_tilt_y * dz,
+                  z);
+}
+
+inline TVector3 pocaOnTrack(const TVector3 &si, const TVector3 &dir,
+                            const TVector3 &p, const TVector3 &d)
+{
+  TVector3 e = dir;
+  double ee = e.Dot(e);
+  double de = d.Dot(e);
+  double denom = ee - de * de; // = |e|^2 * sin^2(angle between beam and track)
+  double t = 0.0;
+  if (denom > 1e-18 * (ee > 0.0 ? ee : 1.0))
+  {
+    TVector3 w = si - p;
+    double dw = d.Dot(w);
+    double ew = e.Dot(w);
+    t = (de * dw - ew) / denom;
+  }
+  return si + t * e;
+}
+
+inline TVector3 beamVertex(const TVector3 &si, const TVector3 &dir)
+{
+  return pocaOnTrack(si, dir, beamAxisOrigin(), beamAxisDirUnit());
+}
+inline TVector3 beamVertexNominal(const TVector3 &si, const TVector3 &dir)
+{
+  return pocaOnTrack(si, dir, TVector3(0.0, 0.0, 0.0), TVector3(0.0, 0.0, 1.0));
+}
+
+inline double beamPerp(const TVector3 &p)
+{
+  TVector3 p0 = beamAxisOrigin();
+  TVector3 d = beamAxisDirUnit();
+  TVector3 w = p - p0;
+  TVector3 perp = w - d * w.Dot(d);
+  return perp.Mag();
+}
+
+inline void fillBeamProfile(HistPlotter *plotter, const TVector3 &vertex,
+                            const TVector3 &si, const TVector3 &dir, const std::string &tag,
+                            bool axisSafe)
+{
+  const std::string folder = "BeamAxis";
+  const std::string all = "beamAxis_all_";
+  const std::string one = "beamAxis_" + tag + "_";
+  double vx = vertex.X(), vy = vertex.Y(), vz = vertex.Z();
+  double pocaDist = beamPerp(vertex);
+  plotter->Fill2D(all + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder);
+  plotter->Fill2D(all + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder);
+  plotter->Fill1D(all + "pocaDist", 400, 0, 100, pocaDist, folder);
+  // Per-branch copy so anomalies can be traced to a specific reconstruction path.
+  plotter->Fill2D(one + "vertexX_vs_Z", 250, -450, 50, 200, -50, 50, vz, vx, folder);
+  plotter->Fill2D(one + "vertexY_vs_Z", 250, -450, 50, 200, -50, 50, vz, vy, folder);
+  plotter->Fill1D(one + "pocaDist", 400, 0, 100, pocaDist, folder);
+  const double zLo = -440.0, zHi = 40.0;
+  const int nSlice = 16;
+  const double sliceW = (zHi - zLo) / nSlice; // 40 mm
+  if (vz >= zLo && vz < zHi)
+  {
+    int is = static_cast<int>((vz - zLo) / sliceW);
+    if (is >= 0 && is < nSlice)
+    {
+      char buf[8];
+      snprintf(buf, sizeof(buf), "%02d", is);
+      plotter->Fill2D(all + "vertexXY_z" + buf, 200, -50, 50, 200, -50, 50, vx, vy, folder);
+    }
+  }
+
+  double dz = dir.Z();
+  double vzBin = beamVertexNominal(si, dir).Z();
+  if (axisSafe && TMath::Abs(dz) > 1e-6 && vzBin >= zLo && vzBin < zHi)
+  {
+    int k = static_cast<int>((vzBin - zLo) / sliceW);
+    if (k >= 0 && k < nSlice)
+    {
+      double zPlane = zLo + (k + 0.5) * sliceW;
+      double s = (zPlane - si.Z()) / dz;
+      double cx = si.X() + s * dir.X();
+      double cy = si.Y() + s * dir.Y();
+      if (TMath::Abs(cx) < 100.0 && TMath::Abs(cy) < 100.0)
+      {
+        char kbuf[8];
+        snprintf(kbuf, sizeof(kbuf), "%02d", k);
+        plotter->Fill1D(all + "crossX_z" + kbuf, 800, -100, 100, cx, folder);
+        plotter->Fill1D(all + "crossY_z" + kbuf, 800, -100, 100, cy, folder);
+        plotter->Fill1D(one + "crossX_z" + kbuf, 800, -100, 100, cx, folder);
+        plotter->Fill1D(one + "crossY_z" + kbuf, 800, -100, 100, cy, folder);
+      }
+    }
+  }
 }
 
 struct PCPath
@@ -587,7 +678,16 @@ void TrackRecon::Begin(TTree * /*tree*/)
     beam_axis_x = std::atof(getenv("BEAM_AXIS_X"));
   if (getenv("BEAM_AXIS_Y"))
     beam_axis_y = std::atof(getenv("BEAM_AXIS_Y"));
-  std::cout << "Beam-axis origin (x,y) = (" << beam_axis_x << ", " << beam_axis_y << ") mm" << std::endl;
+  if (getenv("BEAM_AXIS_Z0"))
+    beam_axis_z0 = std::atof(getenv("BEAM_AXIS_Z0"));
+  if (getenv("BEAM_TILT_X"))
+    beam_tilt_x = std::atof(getenv("BEAM_TILT_X"));
+  if (getenv("BEAM_TILT_Y"))
+    beam_tilt_y = std::atof(getenv("BEAM_TILT_Y"));
+  std::cout << "Beam axis line: (x,y) = (" << beam_axis_x << ", " << beam_axis_y
+            << ") mm at z0 = " << beam_axis_z0
+            << " mm, tilts (dx/dz, dy/dz) = (" << beam_tilt_x << ", " << beam_tilt_y << ")"
+            << std::endl;
   if (doPCEnergyCalibration)
     std::cout << "PC energy calibration ON: alpha source = " << alpha_source_mev
               << " MeV, source position = (" << beam_axis_x << ", " << beam_axis_y << ", " << source_vertex
@@ -2045,6 +2145,8 @@ void protonAlphaHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_E
         // purely longitudinal direction.
         TVector3 r_rhoMin_fix = beamVertex(x1, v);
         double vertex_z = r_rhoMin_fix.Z();
+        fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "apCoinc",
+                        pcevent.multi1 == 1 && pcevent.multi2 == 2);
         double theta_q = (qqqevent.pos - beamAxisPoint(vertex_z)).Theta();
         // double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta();
         double sinTheta_customV = TMath::Sin(theta_q);
@@ -2485,6 +2587,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
         // purely longitudinal direction.
         TVector3 r_rhoMin_fix = beamVertex(x1, v);
+        fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "sx3a1c2", true); // a1c2 only, gated above
         plotter->Fill1D("VertexRecon_pczfix_sx3", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
         plotter->Fill1D("VertexRecon_pczfix", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
         plotter->Fill1D("pczfix_A1C2_1d_sx3", 600, -200, 200, pcz_fix, "PCZ_Recon");
@@ -3036,6 +3139,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
           // purely longitudinal direction.
           TVector3 r_rhoMin_fix = beamVertex(x1, v);
+          fillBeamProfile(plotter, r_rhoMin_fix, x1, v, "qqqa1c2", true); // a1c2 only, gated above
 
           double sinTheta_customV = TMath::Sin((qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta());
           plotter->Fill2D("dE3_E_CathodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
@@ -3593,6 +3697,8 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
       TVector3 x2f(pcXY.X(), pcXY.Y(), pcz_fix);
       TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos);
       double vertex_z = r_rhoMin_fix.Z();
+      const bool axisSafe = (multi2 == 2) || (multi1 == 2 && multi2 == 0);
+      fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos, "elastic_" + det, axisSafe);
       if (vertex_z < z_entrance)
         return;
       double theta = (sievent.pos - r_rhoMin_fix).Theta();
@@ -3629,7 +3735,7 @@ void protonAlphaElastic_core(HistPlotter *plotter, const std::vector<Event> &Si_
         plotter->Fill1D(rx + "_VertexReconZ" + ejtag + sfx, 800, -400, 400, vertex_z, pmlabel);
         plotter->Fill2D(rx + "_VertexReconXY" + ejtag + sfx, 200, -100, 100, 200, -100, 100, r_rhoMin_fix.X(), r_rhoMin_fix.Y(), pmlabel);
         plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, 10, theta * 180 / M_PI, Efix, pmlabel);
-        plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 360, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel);
+        plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -10, 10, theta * 180 / M_PI, Ex, pmlabel);
         plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 180, -180, 180, 800, -10, 10, sievent.pos.Phi() * 180 / M_PI, Ex, pmlabel);
 
         for (const auto &pcevent : PC_Events)
@@ -3970,6 +4076,9 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
       TVector3 x2f(pcXY.X(), pcXY.Y(), pcz_fix);
       TVector3 r_rhoMin_fix = beamVertex(sievent.pos, x2f - sievent.pos);
       double vertex_z = r_rhoMin_fix.Z();
+      const bool axisSafe = (topo1 == "a1c2fix") || (topo1 == "a2c0") || (topo2 == "a1c1_inband");
+      fillBeamProfile(plotter, r_rhoMin_fix, sievent.pos, x2f - sievent.pos,
+                      "reaction_" + globaltag + "_" + det, axisSafe);
       if (beamPerp(r_rhoMin_fix) > perp_cut || vertex_z < z_entrance)
         return;
 
@@ -4029,7 +4138,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
           plotter->Fill1D(rx + "_Ex_from" + ejtag + t + sfx, 600, -10, 20, Ex, pmlabel);
           plotter->Fill2D(rx + "_VertexReconZ_vs_Ef" + ejtag + t + sfx, 800, -400, 400, 800, 0, ef_max, vertex_z, Efix, pmlabel);
           plotter->Fill2D(rx + "_VertexReconZ_vs_Ex" + ejtag + t + sfx, 800, -400, 400, 600, -10, 20, vertex_z, Ex, pmlabel);
-          plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 360, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel);
+          plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + t + sfx, 720, 0, 180, 600, -10, 20, theta * 180 / M_PI, Ex, pmlabel);
 
           if (ebeam_kin_MeV > 0.0)
             plotter->Fill2D(rx + "_BeamEnergy_ETrack_vs_EKin" + ejtag + t + sfx, 400, 0, beamE0 * 1.5, 400, 0, beamE0 * 1.5,
@@ -4069,6 +4178,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
         }
         plotter->Fill1D(rx + "_pczfix" + sfx, 600, -300, 300, pcz_fix, pmlabel);
         plotter->Fill2D(rx + "_Ef_vs_theta" + ejtag + sfx, 100, 0, 180, 800, 0, ef_max, theta * 180 / M_PI, Efix, pmlabel);
+        plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 720, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel);
         plotter->Fill2D(rx + "_Ex_vs_phi" + ejtag + sfx, 45, -180, 180, 600, -10, 20, phi * 180 / M_PI, Ex, pmlabel);
 
         for (const auto &pcevent : PC_Events)
@@ -4107,7 +4217,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
             plotter->Fill2D(rx + "_dEgasCalib_vs_VertexZ" + ejtag + sfx, 800, -400, 400, 800, 0, 0.6, vertex_z, anodeE_MeV, pmlabel);
             plotter->Fill2D(rx + "_dEgasRaw_vs_VertexZ" + ejtag + sfx, 800, -400, 400, 800, 0, 20000, vertex_z, anodeE, pmlabel);
             plotter->Fill2D(rx + "_dEgasRaw_vs_theta" + ejtag + sfx, 180, 0, 180, 800, 0, 20000, theta * 180 / M_PI, anodeE, pmlabel);
-            plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 360, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel);
+            plotter->Fill2D(rx + "_dEgasCalib_vs_theta" + ejtag + sfx, 720, 0, 180, 800, 0, 0.6, theta * 180 / M_PI, anodeE_MeV, pmlabel);
             plotter->Fill2D(rx + "_dEgasCalib_vs_phi" + ejtag + sfx, 90, -180, 180, 800, 0, 0.6, phi * 180 / M_PI, anodeE_MeV, pmlabel);
             // if (anodeCh >= 0)
             //   plotter->Fill2D(rx + "_dEgasCalib_vs_E" + ejtag + sfx + "_anode" + pad2(anodeCh), 400, 0, ef_max, 800, 0, 0.6, sievent.Energy1, anodeE_MeV, pmlabel);

run_27Al.sh

@@ -11,14 +11,12 @@ export CATHODE_GAIN=3.0
 export source_vertex=-200.0
 export DEDX_SCALE=0.89
 export CUTLIST=cuts_list.txt
-export BEAM_AXIS_X=0
+
-export BEAM_AXIS_Y=0
+export BEAM_AXIS_Z0=-200
-# Tilted-beam parameters: (x, y) at BEAM_AXIS_Z0 (default 0), plus dx/dz and dy/dz slopes.
+export BEAM_AXIS_X=-0.136398
-# Default (zeros) reproduces the original vertical beam bit-for-bit. Fit them
+export BEAM_AXIS_Y=0.034824
-# with scratch/iterate_beam.sh + scratch/FitBeamAxis.C.
+export BEAM_TILT_X=-0.00083798
-export BEAM_AXIS_Z0=${BEAM_AXIS_Z0:-0}
+export BEAM_TILT_Y=0.00191377
-export BEAM_TILT_X=${BEAM_TILT_X:-0}
-export BEAM_TILT_Y=${BEAM_TILT_Y:-0}
 
 echo "Pre-compiling TrackRecon.C safely on a single core..."
 root -q -l -b -e '.L TrackRecon.C++O'
@@ -64,11 +62,11 @@ export -f process_run
         
         echo "Starting parallel processing..."
         # time parallel --bar -j 12 process_run ::: {24..41} 44 45 46 {50..59}
-        time parallel --bar -j 12 process_run ::: {24..41} 
+        # time parallel --bar -j 12 process_run ::: {24..41} 
         # time parallel --bar -j 10 process_run ::: 44 45 46 {50..59}
         # time parallel --bar -j 1 process_run ::: 48 # pc without coincidence
         # mv "${CURRENT_OUT_DIR}/results_run048.root" "Output_27Al_run48/."
-        # time parallel --bar -j 10 process_run ::: 62 63 66 67 68 73 74 {78..89}
+        time parallel --bar -j 10 process_run ::: {24..41} 44 45 46 {50..59} 62 63 66 67 68 73 74 {78..89}
         # time parallel --bar -j 4 process_run ::: 62 63 66 67 68
         # time parallel --bar -j 1 process_run ::: 73
         # time parallel --bar -j 1 process_run ::: 74

scratch/FitBeamAxis.C

@@ -0,0 +1,261 @@
+// FitBeamAxis.C
+// ---------------------------------------------------------------------------
+// Single-pass beam-axis line fit. Replaces the iterative
+// FitBeamAxis.C + iterate_beam.sh scheme, which diverges.
+//
+// WHY THE ITERATIVE SCHEME FAILS
+// ------------------------------
+// FitBeamAxis.C fits beamAxis_all_vertexX_vs_Z and feeds the residual back as a
+// correction to BEAM_AXIS_X / BEAM_TILT_X. But that vertex comes from
+// beamVertex(), which returns the track's point of closest approach to the
+// ASSUMED beam line -- so the reconstructed vertex is pulled onto whatever axis
+// the run was configured with, by construction. Its residual therefore measures
+// acceptance asymmetry and geometry artifacts far more than it measures "how far
+// off is my assumed axis". Feeding that back adds a roughly fixed bias each pass
+// with no restoring force, so the parameters drift monotonically and never
+// converge. Measured behaviour: step L2 grew 1.48 -> 2.65 over 6 passes at full
+// step, and still grew (1.27 -> 1.50 over 10 passes) at relax=0.3 -- damping only
+// scales the drift rate, which is the signature of a wrong-direction update
+// rather than overshoot.
+//
+// WHAT THIS DOES INSTEAD
+// ----------------------
+// TrackRecon.C's fillBeamProfile() also histograms where each track CROSSES a set
+// of fixed z-planes (beamAxis_<tag>_crossX_z**, crossY_z**). A track crossing is a
+// property of the track alone -- it never references the assumed beam axis. Tracks
+// originate on the true beam line and fan out with roughly uniform azimuth, so at
+// each plane the crossing distribution centres on the true beam position at that z.
+// Take a robust centre per plane, line-fit those centres against z, and the axis
+// falls out directly. One pass, no iteration, no feedback.
+//
+// WHICH TRACKS FEED THE PLANES
+// ----------------------------
+// fillBeamProfile() only sends a track to the crossing planes when its PC z is a
+// real measurement: a1c2 (cathode charge division), a2c0 (two wires, geometry
+// fixes z), and a1c1 whose a1c1_cfrac_pcz() solve landed in band. Bare a1c1 and
+// a1c0 are excluded -- their z is a Gaussian dither about the raw wire z, which
+// puts them on the assumed axis by construction and would make this fit return
+// whatever axis the run was configured with. The vertexX_vs_Z and pocaDist
+// diagnostics in the same folder are still filled for every topology, so their
+// entry counts will not match the crossing planes'. Do not fit those.
+//
+// CAVEAT -- azimuthal symmetry. The per-plane centre is only unbiased if
+// acceptance is uniform in phi. Dead channels break that (this analysis already
+// excludes QQQ wedge 48/49, ring 63, and several anode wires). The median is used
+// rather than the mean to blunt the resulting asymmetric tails, but a strongly
+// lopsided acceptance will still pull the centre. If the fitted axis looks
+// implausible, check the crossX/crossY distributions at a few planes for
+// visible asymmetry before trusting the number.
+//
+// USAGE:
+//   root -l -b -q 'scratch/FitBeamAxisCrossings.C("Output_27Al/Output_27Al.root")'
+//   root -l -b -q 'scratch/FitBeamAxisCrossings.C("out.root", -440, 40, "all")'
+//   root -l -b -q 'scratch/FitBeamAxisCrossings.C("out.root", -440, 40, "elastic")'
+//
+// Arguments:
+//   filename  TrackRecon.C output ROOT file
+//   zmin/zmax only planes whose z falls in this window are used in the line fit.
+//             Defaults span every plane; narrow it to exclude planes you have
+//             reason to distrust, not to chase a better chi2.
+//   tag       which branch's crossings to use. "all" pools everything; the
+//             per-branch tags are "reaction_<channel>_<det>" (e.g.
+//             "reaction_m27Alax_qqq", "..._sx3") and "elastic_<DET>" ("elastic_QQQ",
+//             "elastic_SX3"). Fitting the qqq and sx3 tags separately is the check
+//             that a tilt is the beam and not one branch's acceptance: same beam,
+//             different geometry and phi coverage, so they must agree.
+//             Which tags exist at all depends on run configuration -- on a dataset
+//             where doPCSX3ClusterAnalysis / doPCQQQClusterAnalysis /
+//             process_alpha_proton_scattering are false and OUT_DIR is not
+//             Output_p, the sx3a1c2, qqqa1c2, apCoinc and elastic_* tags are never
+//             filled and "all" is exactly the reaction_* tags pooled.
+//
+// NOTE on chi2/ndf: medianAndError reports a STATISTICS-ONLY uncertainty
+// (~8 um at a 4e5-entry plane) while observed plane-to-plane scatter in a
+// well-behaved region is ~0.4 mm, i.e. ~50x larger. chi2/ndf is therefore
+// inflated by a missing systematic term and will not approach 1 even for a
+// correct axis. Judge the fit on whether the medians lie on a straight line,
+// not on chi2/ndf alone.
+//
+// Prints the fitted BEAM_AXIS_X/Y and BEAM_TILT_X/Y as ready-to-paste export
+// lines. These are ABSOLUTE values, not increments -- set them directly in
+// run_27Al.sh; do NOT add them to the current values.
+// ---------------------------------------------------------------------------
+
+#include <TFile.h>
+#include <TH1.h>
+#include <TKey.h>
+#include <TDirectory.h>
+#include <TGraphErrors.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TString.h>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+#include <string>
+#include <vector>
+
+// Must match fillBeamProfile() in TrackRecon.C. The range follows the physical
+// vertex acceptance (z_entrance = -454 mm up to the +100 mm cut), not a
+// symmetric window about 0, and stops well below the QQQ plane at z = 105 --
+// see the comment in fillBeamProfile() for why planes near a detector are
+// unusable.
+static const double kZLo = -440.0;
+static const double kZHi = 40.0;
+static const int kNSlice = 16;
+
+static TH1 *findTH1(TDirectory *dir, const std::string &name)
+{
+  if (TObject *o = dir->Get(name.c_str()))
+    if (o && o->InheritsFrom(TH1::Class()) && !o->InheritsFrom("TH2"))
+      return static_cast<TH1 *>(o);
+  TIter next(dir->GetListOfKeys());
+  while (TKey *key = static_cast<TKey *>(next()))
+  {
+    TObject *obj = key->ReadObj();
+    if (!obj)
+      continue;
+    if (obj->InheritsFrom(TDirectory::Class()))
+    {
+      if (TH1 *h = findTH1(static_cast<TDirectory *>(obj), name))
+        return h;
+    }
+    else if (obj->InheritsFrom(TH1::Class()) && !obj->InheritsFrom("TH2") && name == obj->GetName())
+      return static_cast<TH1 *>(obj);
+  }
+  return nullptr;
+}
+
+// Robust centre of a crossing distribution: the median. Chosen over the mean
+// because the crossing distribution has long tails (tracks originating far from
+// this plane) that are not necessarily symmetric once dead channels bite.
+// Uncertainty on the median ~ 1.253 * sigma / sqrt(N), with sigma taken from the
+// interquartile range (IQR/1.349) so a few wild tails don't inflate it.
+static bool medianAndError(TH1 *h, double &med, double &err, double &nEff, int minEntries = 500)
+{
+  if (!h)
+    return false;
+  double n = h->Integral();
+  if (n < minEntries)
+    return false;
+  double probs[3] = {0.25, 0.5, 0.75};
+  double q[3] = {0, 0, 0};
+  h->GetQuantiles(3, q, probs);
+  med = q[1];
+  double sigma = (q[2] - q[0]) / 1.349; // IQR -> Gaussian-equivalent sigma
+  if (sigma <= 0.0)
+    return false;
+  err = 1.253 * sigma / TMath::Sqrt(n);
+  nEff = n;
+  return true;
+}
+
+static bool fitOneProjection(TFile *f, const std::string &base, const char *coord,
+                             double zmin, double zmax, double z0,
+                             double &axis0, double &tilt, double &axis0_err, double &tilt_err,
+                             double &chi2ndf, int &nUsed)
+{
+  std::vector<double> zs, cs, ces, zes;
+  std::cout << "  plane z (mm)   median (mm)      +/-       entries\n";
+  for (int k = 0; k < kNSlice; k++)
+  {
+    double sliceW = (kZHi - kZLo) / kNSlice;
+    double zPlane = kZLo + (k + 0.5) * sliceW;
+    if (zPlane < zmin || zPlane > zmax)
+      continue;
+    char kbuf[8];
+    snprintf(kbuf, sizeof(kbuf), "%02d", k);
+    TH1 *h = findTH1(f, base + "cross" + coord + "_z" + kbuf);
+    double med = 0, err = 0, nEff = 0;
+    if (!medianAndError(h, med, err, nEff))
+    {
+      std::cout << "  " << zPlane << "\t(skipped -- missing or too few entries)\n";
+      continue;
+    }
+    std::cout << "  " << zPlane << "\t\t" << med << "\t" << err << "\t" << (long)nEff << "\n";
+    zs.push_back(zPlane);
+    zes.push_back(0.0);
+    cs.push_back(med);
+    ces.push_back(err);
+  }
+  nUsed = static_cast<int>(zs.size());
+  if (nUsed < 3)
+  {
+    std::cout << "  [" << coord << "] only " << nUsed << " usable planes -- need >= 3 for a line fit\n";
+    return false;
+  }
+
+  TGraphErrors g(nUsed, zs.data(), cs.data(), zes.data(), ces.data());
+  TF1 fL("fL", Form("[0]+[1]*(x-%.8g)", z0), zmin, zmax);
+  fL.SetParameters(cs[nUsed / 2], 0.0);
+  if (g.Fit(&fL, "QRN") != 0)
+  {
+    std::cout << "  [" << coord << "] line fit failed\n";
+    return false;
+  }
+  axis0 = fL.GetParameter(0);
+  tilt = fL.GetParameter(1);
+  axis0_err = fL.GetParError(0);
+  tilt_err = fL.GetParError(1);
+  chi2ndf = (fL.GetNDF() > 0) ? fL.GetChisquare() / fL.GetNDF() : -1.0;
+  return true;
+}
+
+void FitBeamAxis(const char *filename, double zmin = -440.0, double zmax = 40.0,
+                          const char *tag = "all")
+{
+  TFile *f = TFile::Open(filename, "READ");
+  if (!f || f->IsZombie())
+  {
+    std::cerr << "Cannot open file: " << filename << std::endl;
+    return;
+  }
+
+  double z0 = 0.0;
+  if (const char *s = std::getenv("BEAM_AXIS_Z0"))
+    z0 = std::atof(s);
+
+  std::string base = std::string("beamAxis_") + tag + "_";
+
+  std::cout << "\n============ Beam-axis fit (track crossings) ============\n";
+  std::cout << "file:  " << filename << "\n";
+  std::cout << "tag:   " << tag << "   (histogram prefix: " << base << ")\n";
+  std::cout << "range: " << zmin << " <= z <= " << zmax << " mm\n";
+  std::cout << "z0:    " << z0 << " mm  (BEAM_AXIS_Z0)\n\n";
+
+  double bx = 0, tx = 0, bxe = 0, txe = 0, chi2x = -1;
+  double by = 0, ty = 0, bye = 0, tye = 0, chi2y = -1;
+  int nx = 0, ny = 0;
+
+  std::cout << "--- X crossings ---\n";
+  bool okx = fitOneProjection(f, base, "X", zmin, zmax, z0, bx, tx, bxe, txe, chi2x, nx);
+  std::cout << "\n--- Y crossings ---\n";
+  bool oky = fitOneProjection(f, base, "Y", zmin, zmax, z0, by, ty, bye, tye, chi2y, ny);
+
+  std::cout << "\n--- fitted beam line ---\n";
+  if (okx)
+    std::cout << "  x(z0) = " << bx << " +/- " << bxe << " mm,  dx/dz = " << tx << " +/- " << txe
+              << "  (" << tx * 1000.0 << " mrad, chi2/ndf = " << chi2x << ", " << nx << " planes)\n";
+  if (oky)
+    std::cout << "  y(z0) = " << by << " +/- " << bye << " mm,  dy/dz = " << ty << " +/- " << tye
+              << "  (" << ty * 1000.0 << " mrad, chi2/ndf = " << chi2y << ", " << ny << " planes)\n";
+
+  if (!okx || !oky)
+  {
+    std::cout << "\nOne or both projections failed -- not emitting exports.\n";
+    std::cout << "=========================================================\n\n";
+    f->Close();
+    return;
+  }
+
+  std::cout << "\n--- set these directly (ABSOLUTE values, do not add to current) ---\n";
+  std::cout << "export BEAM_AXIS_X=" << bx << "\n";
+  std::cout << "export BEAM_AXIS_Y=" << by << "\n";
+  std::cout << "export BEAM_AXIS_Z0=" << z0 << "\n";
+  std::cout << "export BEAM_TILT_X=" << tx << "\n";
+  std::cout << "export BEAM_TILT_Y=" << ty << "\n";
+  std::cout << "=========================================================\n\n";
+
+  f->Close();
+}