From c0c49be0c9cef37e4c332549f0bfe2cef504012b Mon Sep 17 00:00:00 2001
From: vsitaraman <vs19g@fsu.edu>
Date: Tue, 18 Aug 2026 15:14:18 -0400
Subject: [PATCH] 	modified:   TrackRecon.C consolidated the ditehr
 sigma, the timing cuts and reset the beam xy to (0,0), beamenergy gets set to
 0.001 instead of returning for values that are beyond calculated catima
 range.     Moved some diagnostic plots because they were out of range for the
 boleans to actually catch them, lifted the si time ditehring so that it gets
 applied everywhere for smootehr time plots. 	modified:  
 eloss_calculations/Eloss.py introduced a external scaling factr to allow
 eloss calculations to be paired with the run scripts to figure out the right
 dEdx scaling 	modified:   run_27Al.sh

---
 TrackRecon.C                | 481 +++++++++++++++++++++---------------
 eloss_calculations/Eloss.py |   8 +-
 run_27Al.sh                 |  20 +-
 3 files changed, 306 insertions(+), 203 deletions(-)

diff --git a/TrackRecon.C b/TrackRecon.C
index 8f7224d..fea2e60 100644
--- a/TrackRecon.C
+++ b/TrackRecon.C
@@ -61,9 +61,12 @@ bool process_alpha_proton_scattering = false,
 // 1175mm instead of 1105, plus some part of the window actually lies outside the chamber
 double source_vertex = 53.0,
        z_entrance = -174.3 - 9.7 - 270.0,
-       dither_sigma = 8.0,
-       dither_sigma_c0 = 16.0,
-       cathode_gain = 1.0,
+       dither_sigma = 8.0, // single pcz dither width; A1C0/A1C1/A1C2 all use this.
+                           // dither_sigma_c0 (=16.0, always consumed as /2.0) was
+                           // removed -- numerically identical at the default but it
+                           // desynced from dither_sigma the moment DITHER_SIGMA was
+                           // set, giving the QQQ and SX3 twins a silent 2x difference.
+    cathode_gain = 1.0,
        a1c1_cfrac_split = 0.0,
        a1c1_missing_fmax = 2.0,
        a1c1_lowband_rfactor = 0.0,
@@ -78,7 +81,6 @@ double source_vertex = 53.0,
 
 // --- Immutable Constants ---
 const double qqq_z = 105.0,
-             anode_gain = 1.5146e-5,
              sx3_phi_pitch = 6.5 * (M_PI / 180.0),
              qqq_wedge_pitch = (87.0 / 16.0) * (M_PI / 180.0),
              qqq_ring_pitch = 48.0 / 16.0;
@@ -87,6 +89,25 @@ int co2pc = 3;      // default to 3% CO2; also selects the Eloss table pc suffix
 int pressure = 250; // gas pressure (torr) for the Eloss-table filenames;
                     // overridable via the pressure_in_torr env var.
 
+// One analysis-wide RNG. Previously every dithering/smearing site declared its
+// own `static TRandom3 x(0)`, and ROOT reads seed 0 as "seed from a TUUID" -- so
+// each of the 11 generators picked a fresh stream on every run and the same input
+// file produced different dithered histograms each time, making it impossible to
+// separate a real change from dither noise. Fixed default seed, overridable via
+// RNG_SEED when an independent stream is genuinely wanted.
+TRandom3 anasenRandom(4357);
+
+// Si <-> PC time coincidence. Kept in one place because this gate was previously
+// spelled five different ways (`< 0`, `< 150`, `< -200`, `> 150`-reject,
+// `!(< 150)`-reject) across 15 sites, which is how a sign inversion went unnoticed.
+// One-sided by design: the real coincidence band sits well below zero (see the
+// DelT_Vs_*ECal diagnostics), so only the late side needs rejecting.
+constexpr double kSiPcDtMax = 150.0;
+inline bool siPcCoincident(double t_si, double t_pc)
+{
+  return (t_si - t_pc) < kSiPcDtMax;
+}
+
 inline TVector3 beamVertex(const TVector3 &si, const TVector3 &dir)
 {
   double d = dir.X() * dir.X() + dir.Y() * dir.Y();
@@ -97,6 +118,14 @@ inline double beamPerp(const TVector3 &p)
 {
   return TMath::Sqrt((p.X() - beam_axis_x) * (p.X() - beam_axis_x) + (p.Y() - beam_axis_y) * (p.Y() - beam_axis_y));
 }
+// 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);
+}
 
 struct PCPath
 {
@@ -144,7 +173,7 @@ static const double a1c1_k_17F[7] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
 static const double a1c1_cfmin_27Al[7] = {0.42, 0.42, 0.42, 0.40, 0.42, 0.43, 0.43};
 static const double a1c1_k_27Al[7] = {0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06};
 
-//low band for 17F data
+// low band for 17F data
 
 static const double a1c1_cfmin2_17F[7] = {0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10};
 static const double a1c1_k2_17F[7] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
@@ -474,10 +503,15 @@ void TrackRecon::Begin(TTree * /*tree*/)
   if (getenv("DITHER_SIGMA"))
   {
     dither_sigma = std::atof(getenv("DITHER_SIGMA"));
-    dither_sigma_c0 = dither_sigma;
     std::cout << "Dither Sigma set to " << dither_sigma << " mm" << std::endl;
   }
 
+  if (getenv("RNG_SEED"))
+    anasenRandom.SetSeed(std::atoi(getenv("RNG_SEED")));
+  std::cout << "RNG seed = " << anasenRandom.GetSeed()
+            << " (fixed by default so dithered/smeared histograms are reproducible;"
+            << " set RNG_SEED to vary it, RNG_SEED=0 for a per-run random stream)" << std::endl;
+
   if (getenv("CATHODE_GAIN"))
     cathode_gain = std::atof(getenv("CATHODE_GAIN"));
 
@@ -491,10 +525,9 @@ void TrackRecon::Begin(TTree * /*tree*/)
     std::cout << (excludeBadAnodeWires ? "(active)" : "(list defined but not active)") << 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
-              << ") mm -- appends raw calibration points to pc_calib_raw/ in Terminate()" << std::endl;
+  // (the PC-energy-calibration banner is printed once, further down, after
+  // beam_axis_x/y have been read from the environment -- printing it here too
+  // duplicated the line and reported the pre-override beam axis.)
 
   const double *cfmin_src = a1c1_cfmin_17F;
   const double *k_src = a1c1_k_17F;
@@ -842,7 +875,7 @@ inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events,
     {
       if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) > phi_win)
         return;
-      if (sievent.Time1 - pcevent.Time1 > 150) // time coincidence
+      if (!siPcCoincident(sievent.Time1, pcevent.Time1))
         return;
       double theta = (sievent.pos - source_pos).Theta();
       if (theta <= 0.0 || !std::isfinite(theta))
@@ -917,7 +950,10 @@ inline void pcEnergyCalibrationAccumulateProton(const std::vector<Event> &PC_Eve
     double beam_energy_at_vertex = evalElossForward(MeV_to_cm_p_spl, cm_to_MeVp_spl, initial_energy, beam_path_length);
     beam_energy_at_vertex = applyTaFoilEloss(beam_energy_at_vertex, vertex.Z());
     if (beam_energy_at_vertex <= 0.0)
-      return;
+      beam_energy_at_vertex = 0.001; // clamp rather than drop, matching protonMiscHistograms
+                                     // and reaction_ax_core: a ranged-out beam should show up
+                                     // at the bottom of the spectrum, not vanish and look like
+                                     // the end of the data. Gate it away downstream.
     apkin_a.setValues(mass_1H, mass_4He, mass_4He, mass_1H, beam_energy_at_vertex / mass_1H);
 
     double predicted_alpha_E = predictElasticEnergy(apkin_a, theta * 180.0 / M_PI);
@@ -1056,8 +1092,20 @@ Bool_t TrackRecon::Process(Long64_t entry)
   qqq.CalIndex();
   pc.CalIndex();
 
-  static TRandom3 rnd_qqq(0);
-  static TRandom3 rnd_sx3(0);
+  // Cathode wires 16-23 sit on a different module and read ~300 ns late. This used
+  // to be applied further down, inside the gain-matching loop and only for hits with
+  // pc.e > 50, which meant every timing plot filled before that point (the QQQ and
+  // SX3 dt diagnostics, PCCQQQTimeCut) compared against uncorrected cathode times
+  // while anodeT/cathodeT and the cluster analyses saw corrected ones. Applied once
+  // here, before any consumer, so a cathode dt means the same thing everywhere.
+  for (int i = 0; i < pc.multi; i++)
+  {
+    if (pc.index[i] >= 24 && pc.index[i] - 24 > 15)
+      pc.t[i] -= 300;
+  }
+
+  TRandom3 &rnd_qqq = anasenRandom; // shared stream (see anasenRandom)
+  TRandom3 &rnd_sx3 = anasenRandom;
 
   std::vector<Event> SX3_Events;
   if (sx3.multi > 1)
@@ -1162,7 +1210,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
           plotter->Fill1D("sx3backs_calib", 800, 0, 8192, backE, "hCalSX3");
 
           // plotter->Fill2D("SX3CartesianPlot", 200, -100, 100, 200, -100, 100, 88.0*TMath::Cos(phi_n),88.0*TMath::Sin(phi_n), "hCalSX3");
-          plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, 88.0 * TMath::Cos(phi_n), 88.0 * TMath::Sin(phi_n), "hCalSX3");
+          plotter->Fill2D("SX3CartesianPlot" + std::to_string(id), 200, -100, 100, 200, -100, 100, rho_at_strip * TMath::Cos(phi_n), rho_at_strip * TMath::Sin(phi_n), "hCalSX3");
         }
         if (diagnostic_tplots)
         {
@@ -1328,8 +1376,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
             }
             if (diagnostic_eplots)
             {
-              // if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al
-              if (tRing - static_cast<double>(pc.t[k]) < 150) // proton tests, 27Al
+              if (siPcCoincident(tRing, static_cast<double>(pc.t[k])))
               {
                 PCAQQQTimeCut = true;
                 plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
@@ -1340,7 +1387,10 @@ Bool_t TrackRecon::Process(Long64_t entry)
 
           if (pc.index[k] >= 24 && pc.e[k] > 10)
           {
-            if (tRing - static_cast<double>(pc.t[k]) < -200)
+            // unified onto the shared gate (was < -200); the -300 ns correction for
+            // cathode wires 16-23 is now applied before this point, so the anode and
+            // cathode differences are on the same footing.
+            if (siPcCoincident(tRing, static_cast<double>(pc.t[k])))
               PCCQQQTimeCut = true;
             if (diagnostic_tplots)
             {
@@ -1369,7 +1419,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
 
   PCQQQTimeCut = PCAQQQTimeCut && PCCQQQTimeCut;
 #ifdef RAW_HISTOS
-  plotter->Fill1D("QQQ_Multiplicity", 10, 0, 10, qqqCount, "hRawQQQ");
+  plotter->Fill1D("QQQ_Multiplicity", 11, -0.5, 10.5, qqqCount, "hRawQQQ");
 #endif
   aWireEvents.clear();
   aWireEvents.reserve(24);
@@ -1408,8 +1458,8 @@ Bool_t TrackRecon::Process(Long64_t entry)
       }
       else
       {
-        if (pc.index[i] - 24 > 15)
-          pc.t[i] -= 300;
+        // (the -300 ns correction for cathode wires 16-23 is applied once, up near
+        // pc.CalIndex(), so every consumer sees the same corrected time)
         cathodeT = static_cast<double>(pc.t[i]);
         cathodeIndex = pc.index[i] - 24;
         // cWireEvents[pc.index[i] - 24] = std::tuple(pc.index[i] - 24, pc.e[i], static_cast<double>(pc.t[i]));
@@ -1417,49 +1467,14 @@ Bool_t TrackRecon::Process(Long64_t entry)
       }
     }
 
-    if (anodeT != -99999 && cathodeT != 99999)
-    {
-      for (int j = 0; j < qqq.multi; j++)
-      {
-        if (diagnostic_tplots)
-        {
-          plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
-          plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, qqq.ch[j], "hTiming");
-          plotter->Fill2D("PC_Time_vs_AIndex_qqq", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
-          plotter->Fill2D("PC_Time_vs_CIndex_qqq", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
-          // plotter->Fill1D("PC_Time_A" + std::to_string(anodeIndex) + "_C" + std::to_string(cathodeIndex), 200, -1000, 1000, anodeT - cathodeT, "TimingPC");
-        }
-      }
-
-      for (int j = 0; j < sx3.multi; j++)
-      {
-        if (diagnostic_tplots)
-        {
-          plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
-          // plotter->Fill2D("PC_Time_Vs_SX3_ch", 200, -2000, 2000, 16 * 8, 0, 16 * 8, anodeT - cathodeT, sx3.ch[j], "hTiming");
-          plotter->Fill2D("PC_Time_vs_AIndex_sx3", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, anodeIndex, "hTiming");
-          plotter->Fill2D("PC_Time_vs_CIndex_sx3", 200, -2000, 2000, 24, 0, 24, anodeT - cathodeT, cathodeIndex, "hTiming");
-        }
-      }
-      for (const auto &sx3event : SX3_Events)
-      {
-        bool TCC = sx3event.Time1 - cathodeT < 0;
-        bool TCA = sx3event.Time1 - anodeT < 0;
-        // plotter->Fill2D("sx3_z_phi_awire"+std::to_string(anodeIndex)+"_TC"+std::to_string(TCA), 400,-100,100, 200, -200,200,sx3event.pos.Z(), sx3event.pos.Phi()*180/M_PI );
-        // plotter->Fill2D("sx3_z_phi_cwire"+std::to_string(cathodeIndex)+"_TC"+std::to_string(TCC), 400,-100,100, 200, -200,200,sx3event.pos.Z(), sx3event.pos.Phi()*180/M_PI );
-      }
-
-      if (diagnostic_tplots)
-      {
-        plotter->Fill1D("PC_Time", 200, -2000, 2000, anodeT - cathodeT, "hTiming");
-      }
-    }
-
     for (int j = i + 1; j < pc.multi; j++)
     {
 #ifdef RAW_HISTOS
       plotter->Fill2D("PC_Coincidence_Matrix", 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
-      plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
+      // only tag with anodeT-cathodeT once both are real: with the sentinels still in
+      // place the difference is -199998, which piled the early hits into the "_1" plot
+      if (anodeT != -99999 && cathodeT != 99999)
+        plotter->Fill2D("PC_Coincidence_Matrix_anodeMinusCathode_lt_-200_" + std::to_string(anodeT - cathodeT < -200), 48, 0, 48, 48, 0, 48, pc.index[i], pc.index[j], "hRawPC");
 #endif
 
       if (diagnostic_eplots)
@@ -1468,6 +1483,33 @@ Bool_t TrackRecon::Process(Long64_t entry)
       }
     }
   }
+
+  // anodeT - cathodeT is one number per event, so it gets filled once per event.
+  // This block used to sit inside the pc.multi loop above *and* wrap an inner loop
+  // over qqq.multi / sx3.multi, so the identical value was filled O(pc.multi x
+  // qqq.multi) times -- inflating the statistics and weighting every event by its
+  // own multiplicity. PC_Time_Vs_QQQ_ch genuinely needs the channel loop, so it
+  // keeps one; the rest do not.
+  if (diagnostic_tplots && anodeT != -99999 && cathodeT != 99999)
+  {
+    double pcDT = anodeT - cathodeT;
+    plotter->Fill1D("PC_Time", 200, -2000, 2000, pcDT, "hTiming");
+    if (qqq.multi > 0)
+    {
+      plotter->Fill1D("PC_Time_qqq", 200, -2000, 2000, pcDT, "hTiming");
+      plotter->Fill2D("PC_Time_vs_AIndex_qqq", 200, -2000, 2000, 24, -0.5, 23.5, pcDT, anodeIndex, "hTiming");
+      plotter->Fill2D("PC_Time_vs_CIndex_qqq", 200, -2000, 2000, 24, -0.5, 23.5, pcDT, cathodeIndex, "hTiming");
+      for (int j = 0; j < qqq.multi; j++)
+        plotter->Fill2D("PC_Time_Vs_QQQ_ch", 200, -2000, 2000, 16 * 8, -0.5, 16 * 8 - 0.5, pcDT, qqq.ch[j], "hTiming");
+    }
+    if (sx3.multi > 0)
+    {
+      plotter->Fill1D("PC_Time_sx3", 200, -2000, 2000, pcDT, "hTiming");
+      plotter->Fill2D("PC_Time_vs_AIndex_sx3", 200, -2000, 2000, 24, -0.5, 23.5, pcDT, anodeIndex, "hTiming");
+      plotter->Fill2D("PC_Time_vs_CIndex_sx3", 200, -2000, 2000, 24, -0.5, 23.5, pcDT, cathodeIndex, "hTiming");
+    }
+  }
+
   anodeHits.clear();
   cathodeHits.clear();
   corrcatMax.clear();
@@ -1555,9 +1597,25 @@ Bool_t TrackRecon::Process(Long64_t entry)
               anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w);
           }
           PCEventCalibrated.Energy1 = anodeCalibSum;
-          // Cathode uses the single max wire (cpMaxE) -- indexed by z, so it's
-          // already phi-consistent; leave it as-is.
-          PCEventCalibrated.Energy2 = pcEnergySlope[24 + PCEvent.Cathodech] * cpMaxE;
+          // Cathode uses the single max-energy wire (cpMaxE). That wire is NOT
+          // necessarily cCluster[0], which is all PCEvent.Cathodech records, so
+          // pcEnergySlope[24 + Cathodech] was applying the wrong wire's constant to
+          // cpMaxE for every multi-wire cathode cluster -- i.e. for A1C2, the primary
+          // topology. GetPseudoWire tracks the max energy but not its index, so find
+          // it here rather than change that signature for its five call sites.
+          int cMaxWire = PCEvent.Cathodech;
+          double cMaxE = -1.0;
+          for (const auto &w : cCluster)
+          {
+            if (std::get<1>(w) > cMaxE)
+            {
+              cMaxE = std::get<1>(w);
+              cMaxWire = std::get<0>(w);
+            }
+          }
+          PCEventCalibrated.Energy2 = (cMaxWire >= 0 && cMaxWire < 24)
+                                          ? pcEnergySlope[24 + cMaxWire] * cpMaxE
+                                          : cpMaxE;
           PC_Events_calibrated.push_back(PCEventCalibrated);
         }
       }
@@ -1592,11 +1650,12 @@ Bool_t TrackRecon::Process(Long64_t entry)
       {
         for (const auto &si : sis)
         {
-          if (!(si.Time1 - apTSMaxE < 150)) // loose time coincidence (benchmark sign)
+          if (!siPcCoincident(si.Time1, apTSMaxE))
             continue;
           TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire, si.pos.Phi());
           double dphi = TMath::Abs(si.pos.DeltaPhi(pc));
-          if (dphi <= TMath::Pi() / 4.0 && dphi < bestDphi)
+          double phi_win = isQQQ ? TMath::Pi() / 4.0 : TMath::Pi() / 3.0; // per-detector, as elsewhere
+          if (dphi <= phi_win && dphi < bestDphi)
           {
             bestDphi = dphi;
             bestSi = &si;
@@ -1636,7 +1695,6 @@ Bool_t TrackRecon::Process(Long64_t entry)
         ev.Cathodech = -1;
         PC_Events_calibrated.push_back(ev);
       }
-
     }
   }
 
@@ -1648,23 +1706,34 @@ Bool_t TrackRecon::Process(Long64_t entry)
 
   //////Timing stuff for F data
 
-  static TRandom3 rnd(0); // seeded once (random seed via TUUID), not per event
+  TRandom3 &rnd = anasenRandom;
   if (dataset == "17F" && reactiondata)
   {
+    // misc.ch is a property of the event, not of any Si hit -- filled inside the
+    // per-QQQ and per-SX3 loops it was multiplied by the Si multiplicity.
+    for (int j = 0; j < misc.multi; j++)
+    {
+      if (QQQ_Events.size())
+        plotter->Fill1D("channels_misc_qqq", 20, -0.5, 19.5, misc.ch[j], "misc");
+      if (SX3_Events.size())
+        plotter->Fill1D("channels_misc_sx3", 20, -0.5, 19.5, misc.ch[j], "misc");
+    }
     int ctr = 0;
     for (const auto &qqqevent : QQQ_Events)
     {
       double ts_rf = -987654321;
       double ts_needle = -987654321;
       double ts_mcp = -987654321;
-      double ts_qqq = static_cast<double>(qqqevent.Time1) + (rnd.Uniform(16.0) - 8.0);
+      // Time1 already carries the +/-8 clock dither applied when the QQQ Event was
+      // built (tRing), so re-dithering here widened this folder's timing by sqrt(2)
+      // relative to every other timing plot in the analysis.
+      double ts_qqq = static_cast<double>(qqqevent.Time1);
       bool found_rf = false;
       bool found_mcp = false;
       bool found_needle = false;
       bool qqq_inner_ring = (qqqevent.ch1 % 16) < 8;
       for (int j = 0; j < misc.multi; j++)
       {
-        plotter->Fill1D("channels_misc_qqq", 20, 0, 20, misc.ch[j], "misc");
         if (misc.ch[j] == 2)
         { // Needle
           plotter->Fill2D("needle_vs_qqqE", 800, 0, 16384, 800, 0, 10, misc.e[j], qqqevent.Energy1, "misc");
@@ -1704,13 +1773,13 @@ Bool_t TrackRecon::Process(Long64_t entry)
       double ts_rf = -987654321;
       double ts_needle = -987654321;
       double ts_mcp = -987654321;
-      double ts_sx3 = static_cast<double>(sx3event.Time1) + (rnd.Uniform(16.0) - 8.0);
+      // as with ts_qqq: det.ts was already dithered when the SX3 Event was built
+      double ts_sx3 = static_cast<double>(sx3event.Time1);
       bool found_rf = false;
       bool found_mcp = false;
       bool found_needle = false;
       for (int j = 0; j < misc.multi; j++)
       {
-        plotter->Fill1D("channels_misc_sx3", 20, 0, 20, misc.ch[j], "misc");
         if (misc.ch[j] == 2)
         { // Needle
           plotter->Fill2D("needle_vs_sx3E", 800, 0, 16384, 800, 0, 10, misc.e[j], sx3event.Energy1, "misc");
@@ -1756,15 +1825,24 @@ Bool_t TrackRecon::Process(Long64_t entry)
     pcCalibratedHistograms(plotter, QQQ_Events, SX3_Events, PC_Events_calibrated);
 
   a1c1CalibDiagnostic(plotter, PC_Events);                            // <-- new, unconditional
-  pcVertexByWireGeometry(plotter, QQQ_Events, SX3_Events, PC_Events_calibrated); // <-- new, unconditional
+  pcVertexByWireGeometry(plotter, QQQ_Events, SX3_Events, PC_Events); // <-- new, unconditional
 
-  auto hasPCCoincidence = [&](const TVector3 &pos)
+  // phi_win matches every other Si-PC match in this file: SX3 sits at a longer lever
+  // arm (rho ~88mm vs the PC anode at 37mm) than QQQ, so its true phi spread is wider
+  // -- pi/4 for both detectors under-counted real SX3-PC coincidences. Also gated on
+  // siPcCoincident(): this used to be the one Si-PC match in the file with a phi
+  // window but no time gate, so "withPC" included phi-aligned but time-accidental
+  // pairs.
+  auto hasPCCoincidence = [&](const Event &sievent, bool isQQQ)
   {
+    double phi_win = isQQQ ? TMath::Pi() / 4.0 : TMath::Pi() / 3.0;
     for (const auto &pcevent : PC_Events)
     {
       if (pcevent.multi1 < 1)
         continue;
-      if (TMath::Abs(pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0)
+      if (!siPcCoincident(sievent.Time1, pcevent.Time1))
+        continue;
+      if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) <= phi_win)
         return true;
     }
     return false;
@@ -1772,13 +1850,13 @@ Bool_t TrackRecon::Process(Long64_t entry)
   for (const auto &qqqevent : QQQ_Events)
   {
     plotter->Fill1D("siE_qqq_calibrated_all", 800, 0, 15, qqqevent.Energy1, "siE");
-    bool coinc = hasPCCoincidence(qqqevent.pos);
+    bool coinc = hasPCCoincidence(qqqevent, true);
     plotter->Fill1D(coinc ? "siE_qqq_calibrated_withPC" : "siE_qqq_calibrated_noPC", 800, 0, 15, qqqevent.Energy1, "siE");
   }
   for (const auto &sx3event : SX3_Events)
   {
     plotter->Fill1D("siE_sx3_calibrated_all", 800, 0, 15, sx3event.Energy1, "siE");
-    bool coinc = hasPCCoincidence(sx3event.pos);
+    bool coinc = hasPCCoincidence(sx3event, false);
     plotter->Fill1D(coinc ? "siE_sx3_calibrated_withPC" : "siE_sx3_calibrated_noPC", 800, 0, 15, sx3event.Energy1, "siE");
   }
 
@@ -1806,11 +1884,11 @@ Bool_t TrackRecon::Process(Long64_t entry)
   plotter->Fill2D("ac_vs_cc", 20, 0, 20, 20, 0, 20, aClusters.size(), cClusters.size(), "wiremult");
   for (const auto &cluster : aClusters)
   {
-    plotter->Fill1D("aClusters" + std::to_string(aClusters.size()), 20, -5, 15, cluster.size(), "wiremult");
+    plotter->Fill1D("aClusters" + std::to_string(aClusters.size()), 20, -0.5, 19.5, cluster.size(), "wiremult");
   }
   for (const auto &cluster : cClusters)
   {
-    plotter->Fill1D("cClusters" + std::to_string(cClusters.size()), 20, -5, 15, cluster.size(), "wiremult");
+    plotter->Fill1D("cClusters" + std::to_string(cClusters.size()), 20, -0.5, 19.5, cluster.size(), "wiremult");
   }
 
   if (cClusters.size() && aClusters.size())
@@ -1894,13 +1972,15 @@ void protonAlphaHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_E
         TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
         TVector3 x1(qqqevent.pos);
         TVector3 v = x2f - x1;
-        double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
-        TVector3 r_rhoMin_fix = x1 + t_minimum * v;
+        // beamVertex() instead of an inline projection: the hand-rolled version pinned the
+        // 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);
         double vertex_z = r_rhoMin_fix.Z();
-        double theta_q = (qqqevent.pos - TVector3(0, 0, vertex_z)).Theta();
+        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);
-        double theta_s = (sx3event.pos - TVector3(0, 0, vertex_z)).Theta();
+        double theta_s = (sx3event.pos - beamAxisPoint(vertex_z)).Theta();
         // double theta_s = (sx3event.pos - r_rhoMin_fix).Theta();
         double sinTheta_s = TMath::Sin(theta_s);
         // if(vertex_z<0 || vertex_z>100) continue;
@@ -2003,6 +2083,7 @@ void a1c1CalibDiagnostic(HistPlotter *plotter, const std::vector<Event> &PC_Even
       double fracPos = TMath::Abs(fracPos_signed);
 
       plotter->Fill1D(Form("A1C1Calib_cfrac_cell%d", cell), 220, -0.05, 1.05, cfrac, "A1C1Calib");
+      plotter->Fill2D(Form("A1C1Calib_cfrac_vs_AnodeE_cell%d", cell), 220, -0.05, 1.05, 800, 0, 40000, cfrac, pcevent.Energy1, "A1C1Calib");
       plotter->Fill2D("A1C1Calib_cfrac_vs_cellFrac", 120, 0, 1.2, 220, -0.05, 1.05, fracPos, cfrac, "A1C1Calib");
       plotter->Fill2D(Form("A1C1Calib_cfrac_vs_cellFrac_signed_cell%d", cell), 240, -1.2, 1.2, 220, -0.05, 1.05,
                       fracPos_signed, cfrac, "A1C1Calib");
@@ -2013,7 +2094,7 @@ void a1c1CalibDiagnostic(HistPlotter *plotter, const std::vector<Event> &PC_Even
 
 void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events_calibrated)
 {
-  static TRandom3 rand(0); // seeded once, not per call -- dithers A1C0's Z below
+  TRandom3 &rand = anasenRandom; // dithers A1C0's Z below
 
   auto fillFor = [&](const std::vector<Event> &sis, bool isQQQ)
   {
@@ -2047,7 +2128,7 @@ void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_
         else if (pcevent.multi1 == 1 && pcevent.multi2 == 1) // A1C1
           pcz = a1c1_cfrac_pcz(pcevent, si.pos, a1c1_inband);
         else if (pcevent.multi1 == 1 && pcevent.multi2 == 0) // A1C0
-          pcz = rand.Gaus(pcevent.pos.Z(), dither_sigma_c0 / 2.0);
+          pcz = rand.Gaus(pcevent.pos.Z(), dither_sigma);
         else // A2C0 (multi1==2, multi2==0) -- undithered by design
           pcz = pcevent.pos.Z();
 
@@ -2055,7 +2136,7 @@ void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_
         TVector3 vtx = beamVertex(si.pos, x2 - si.pos);
         if (beamPerp(vtx) > perp_max)
           continue;
-        if (vtx.Z() < -173.6 || vtx.Z() > 100)
+        if (vtx.Z() < z_entrance || vtx.Z() > 100)
           continue;
 
         std::string topo = "_a" + std::to_string(pcevent.multi1) + "c" + std::to_string(pcevent.multi2);
@@ -2078,7 +2159,7 @@ void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_
 
 void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events_calibrated)
 {
-  static TRandom3 rand(0); // seeded once, not per call -- for Si-side pixel/strip dithering below
+  TRandom3 &rand = anasenRandom; // for Si-side pixel/strip dithering below
 
   for (const auto &pcevent : PC_Events_calibrated)
   {
@@ -2154,7 +2235,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
         for (const auto &qqqevent : QQQ_Events)
         {
           bool phicut = TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
-          bool timecut = (qqqevent.Time1 - pcevent.Time1) < 150;
+          bool timecut = siPcCoincident(qqqevent.Time1, pcevent.Time1);
           if (!(phicut && timecut))
             continue;
 
@@ -2177,8 +2258,8 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
         }
         for (const auto &sx3event : SX3_Events)
         {
-          bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
-          bool timecut = (sx3event.Time1 - pcevent.Time1) < 150;
+          bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 3.0; // wider lever arm than QQQ, see phi_win elsewhere
+          bool timecut = siPcCoincident(sx3event.Time1, pcevent.Time1);
           if (!(phicut && timecut))
             continue;
 
@@ -2209,7 +2290,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
                           const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters)
 {
 
-  static TRandom3 rand(0);
+  TRandom3 &rand = anasenRandom;
 
   // --- GENUINE A1C0 events:
   if (BenchMark && aClusters.size() == 1 && cClusters.size() == 0)
@@ -2220,9 +2301,9 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
     double anodeTS = std::get<3>(aPw);
     for (const auto &sx3event : SX3_Events)
     {
-      bool PCSX3TimeCut = (sx3event.Time1 - anodeTS < 150);
+      bool PCSX3TimeCut = siPcCoincident(sx3event.Time1, anodeTS);
       TVector3 pc = a1c0_wirePos(apwire_bm, sx3event.pos.Phi(), false);
-      bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pc)) <= TMath::Pi() / 4.0;
+      bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pc)) <= TMath::Pi() / 3.0; // per-detector phi_win, see pcCalibratedHistograms
       if (!(phicut && PCSX3TimeCut))
         continue;
       double smeared_phi = sx3event.pos.Phi() + rand.Uniform(-sx3_phi_pitch / 2.0, sx3_phi_pitch / 2.0);
@@ -2232,9 +2313,9 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       TVector3 vtx0 = beamVertex(sx3event.pos, pc - sx3event.pos);
       TVector3 vtx1 = beamVertex(smeared_sx3, pc_hybrid - smeared_sx3);
 
-      if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
+      if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance))
         continue;
-      double sx3theta = TMath::ATan2(88.0, sx3event.pos.Z() - source_vertex);
+      double sx3theta = TMath::ATan2(sx3event.pos.Perp(), sx3event.pos.Z() - source_vertex); // true per-strip rho, not nominal 88
       double pczguess = 37.0 / TMath::Tan(sx3theta) + source_vertex;
       plotter->Fill1D("Benchmark_SX3_VertexZ_trueA1C0", 800, -400, 400, vtx0.Z(), "A1C0True_SX3");
       plotter->Fill1D("Benchmark_SX3_VertexZ_trueA1C0_Hybrid", 800, -400, 400, vtx1.Z(), "A1C0True_SX3");
@@ -2248,20 +2329,16 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
   for (const auto &pcevent : PC_Events)
   {
-    bool PCSX3TimeCut = false;
-    bool PCASX3TimeCut = false;
-    bool PCCSX3TimeCut = false;
     for (const auto &sx3event : SX3_Events)
     {
       plotter->Fill1D("dt_pcA_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "Timing");
       plotter->Fill1D("dt_pcC_sx3B" + std::to_string(sx3event.ch2), 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "Timing");
-      if (sx3event.Time1 - pcevent.Time1 < 0) //-150 for alphas
-        PCASX3TimeCut = 1;
-      if (sx3event.Time1 - pcevent.Time2 < 0) //-200 for alphas
-        PCCSX3TimeCut = 1;
-      PCSX3TimeCut = PCASX3TimeCut && PCCSX3TimeCut;
 
-      bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
+      bool PCASX3TimeCut = siPcCoincident(sx3event.Time1, pcevent.Time1);
+      bool PCCSX3TimeCut = siPcCoincident(sx3event.Time1, pcevent.Time2);
+      bool PCSX3TimeCut = PCASX3TimeCut && PCCSX3TimeCut;
+
+      bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 3.0; // per-detector phi_win, see pcCalibratedHistograms
 
       plotter->Fill1D("dt_pcA_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time1, "Timing");
       plotter->Fill1D("dt_pcC_sx3B", 640, -2000, 2000, sx3event.Time1 - pcevent.Time2, "Timing");
@@ -2289,8 +2366,8 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       if (pcevent.multi1 == 2 && pcevent.multi2 == 0)
         plotter->Fill2D("dE_E_Cathodesx3B_a2c0", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
 
-      plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < 150), 100, -200, 200, 100, -200, 200, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
-      plotter->Fill1D("sx3phi_minus_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < 150), 100, -180, 180, (sx3event.pos.DeltaPhi(pcevent.pos)) * 180 / M_PI, "Kinematics_Angles");
+      plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(siPcCoincident(sx3event.Time1, pcevent.Time1)), 100, -200, 200, 100, -200, 200, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
+      plotter->Fill1D("sx3phi_minus_pcphi" + std::to_string(siPcCoincident(sx3event.Time1, pcevent.Time1)), 100, -180, 180, (sx3event.pos.DeltaPhi(pcevent.pos)) * 180 / M_PI, "Kinematics_Angles");
 
       if (PCSX3TimeCut)
       {
@@ -2301,7 +2378,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         plotter->Fill2D("pcz_vs_pcphi_TimeCut", 600, -200, 200, 120, -200, 200, pcevent.pos.Z(), pcevent.pos.Phi() * 180 / M_PI, "PCZ_Recon");
       }
 
-      double sx3rho = 88.0;
+      double sx3rho = sx3event.pos.Perp(); // 88/cos(alpha_n) as built, not nominal 88
       double sx3z = sx3event.pos.Z();
       double pcz = pcevent.pos.Z();
       double calcsx3theta = TMath::ATan2(sx3rho - z_to_crossover_rho(pcz), sx3z - pcz);
@@ -2315,8 +2392,10 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
       TVector3 x2(pcevent.pos), x1(sx3event.pos);
       TVector3 v = x2 - x1;
-      double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
-      TVector3 vector_closest_to_z_sx3 = x1 + t_minimum * v;
+      // beamVertex() instead of an inline projection: the hand-rolled version pinned the
+      // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
+      // purely longitudinal direction.
+      TVector3 vector_closest_to_z_sx3 = beamVertex(x1, v);
       plotter->Fill1D("VertexReconZ_SX3" + std::to_string(PCSX3TimeCut), 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "Vertex_Reconstruction");
       plotter->Fill1D("VertexReconZ_SX3", 600, -1300, 1300, vector_closest_to_z_sx3.Z(), "Vertex_Reconstruction");
       plotter->Fill2D("VertexReconXY_SX3" + std::to_string(PCSX3TimeCut), 100, -100, 100, 100, -100, 100, vector_closest_to_z_sx3.X(), vector_closest_to_z_sx3.Y(), "Vertex_Reconstruction");
@@ -2334,8 +2413,10 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
         TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
         TVector3 v = x2f - x1;
-        double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
-        TVector3 r_rhoMin_fix = x1 + t_minimum * v;
+        // beamVertex() instead of an inline projection: the hand-rolled version pinned the
+        // 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);
         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");
@@ -2349,7 +2430,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         plotter->Fill1D("pczfix-sx3pczguess_A1C2", 200, -100, 100, pcz_fix - pczguess, "Residuals");
         plotter->Fill2D("pczfix_vs_sx3pczguess_A1C2_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "PCZ_Recon");
 
-        double sinTheta_customV = TMath::Sin((sx3event.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta());
+        double sinTheta_customV = TMath::Sin((sx3event.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta());
         plotter->Fill2D("dE3_E_CathodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
         plotter->Fill2D("dE3_E_AnodeSX3_A1C2_TC" + std::to_string(PCSX3TimeCut) + "_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
 
@@ -2359,15 +2440,17 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         //   plotter->Fill2D("dE3_E_CathodeSX3B_A1C2_(vertex_fix_z/100)=" + std::to_string(floor(r_rhoMin_fix.Z() / 100.0)), 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
         // }
 
-        if (pcevent.multi1 == 1 && pcevent.multi2 == 3)
-        {
-          plotter->Fill2D("pcz_vs_sx3pczguess_A1C3", 600, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "PCZ_Recon");
-        }
+        // (a `multi2 == 3` branch used to be nested here inside `multi2 == 2`, so
+        // pcz_vs_sx3pczguess_A1C3 could never be filled; removed rather than moved,
+        // since A1C2 is the only cathode topology with a z model.)
 
         plotter->Fill2D("pcz_vs_sx3pczguess_int", 600, -200, 200, 600, -200, 200, pcz_guess_int, pcevent.pos.Z(), "PCZ_Recon");
         // plotter->Fill2D("pcz_vs_sx3pczguess_strip" + std::to_string(sx3event.ch2), 300, -200, 200, 600, -200, 200, pczguess, pcevent.pos.Z(), "PCZ_Recon");
 
-        bool sx3PhiCut = (TMath::Abs(sx3event.pos.Phi() - pcevent.pos.Phi()) < 45.0 * M_PI / 180.);
+        // was a raw Phi() difference, which ignores the +/-pi wrap (a +170/-170 pair
+        // is 20 deg apart but scored as 340). 45 deg is pi/4, so this is exactly the
+        // `phicut` computed above -- reuse it instead of keeping two names for it.
+        bool sx3PhiCut = phicut;
 
         plotter->Fill1D("pcz_sx3Coinc_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, sx3z, "PCZ_Recon");
         plotter->Fill2D("pcz_vs_sx3z_phiCut" + std::to_string(sx3PhiCut) + "_TC" + std::to_string(PCSX3TimeCut), 300, 0, 200, 600, -400, 400, sx3z, pcevent.pos.Z(), "PCZ_Recon");
@@ -2376,11 +2459,18 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
         // plotter->Fill2D("pcdEA_vs_sx3z", 300, 0, 200, 800, 0, 20000, sx3z, pcevent.Energy1, "Kinematics_Angles");
         // plotter->Fill2D("pcdEA_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "Kinematics_Angles");
-        plotter->Fill2D("pcdEA_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z");
+        // gated on the (now per-pair) PCSX3TimeCut to match the QQQ twin: these four names
+        // are shared with PCQQQClusterAnalysis, whose fills sit inside `if (timecut)`, so
+        // leaving the SX3 side ungated made each merged histogram half-gated.
+        if (PCSX3TimeCut)
+          plotter->Fill2D("pcdEA_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z");
         // plotter->Fill2D("pcdEC_vs_sx3z", 300, 0, 200, 800, 0, 20000, sx3z, pcevent.Energy2, "Kinematics_Angles");
         // plotter->Fill2D("pcdEC_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "Kinematics_Angles");
-        plotter->Fill2D("pcdEC_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy2, "PCdE_vs_Z");
-        plotter->Fill2D("pcdEACSum_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z");
+        // gated on the (now per-pair) PCSX3TimeCut to match the QQQ twin: these four names
+        // are shared with PCQQQClusterAnalysis, whose fills sit inside `if (timecut)`, so
+        // leaving the SX3 side ungated made each merged histogram half-gated.
+        if (PCSX3TimeCut)
+          plotter->Fill2D("pcdEC_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy2, "PCdE_vs_Z");
 
         // plotter->Fill2D("pcdEA_vs_sx3z" + std::to_string(sx3event.ch2), 300, 0, 200, 800, 0, 20000, sx3z, pcevent.Energy1, "Kinematics_Angles");
         // plotter->Fill2D("pcdEA_vs_sx3pczguess" + std::to_string(sx3event.ch2), 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "Kinematics_Angles");
@@ -2395,7 +2485,11 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
       // plotter->Fill2D("pcdEA_vs_sx3pczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z");
       plotter->Fill2D("pcdEA_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z");
-      plotter->Fill2D("pcdEA_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z");
+      // gated on the (now per-pair) PCSX3TimeCut to match the QQQ twin: these four names
+      // are shared with PCQQQClusterAnalysis, whose fills sit inside `if (timecut)`, so
+      // leaving the SX3 side ungated made each merged histogram half-gated.
+      if (PCSX3TimeCut)
+        plotter->Fill2D("pcdEA_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy1, "PCdE_vs_Z");
       // plotter->Fill2D("pcdEA_vs_pczfix" + std::to_string(pcevent.multi1) + "A" + std::to_string(pcevent.multi1) + "C", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z");
       // plotter->Fill2D("pcdEC_vs_sx3pczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z");
       if (pcevent.multi1 == 1)
@@ -2411,8 +2505,11 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         plotter->Fill2D("pcdEC_vs_sx3pczguess_C2", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z");
       }
       plotter->Fill2D("pcdEC_vs_sx3pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z");
-      plotter->Fill2D("pcdEC_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z");
-      plotter->Fill2D("pcdEACSum_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, (pcevent.Energy2 + pcevent.Energy1) / 2, "PCdE_vs_Z");
+      // gated on the (now per-pair) PCSX3TimeCut to match the QQQ twin: these four names
+      // are shared with PCQQQClusterAnalysis, whose fills sit inside `if (timecut)`, so
+      // leaving the SX3 side ungated made each merged histogram half-gated.
+      if (PCSX3TimeCut)
+        plotter->Fill2D("pcdEC_vs_pczguess", 600, -200, 200, 800, 0, 20000, pczguess, pcevent.Energy2, "PCdE_vs_Z");
 
       // plotter->Fill2D("pcdEC_vs_pczfix" + std::to_string(pcevent.multi1) + "A" + std::to_string(pcevent.multi1) + "C", 800, 0, 20000, 600, -200, 200, pcevent.Energy2, pcz_fix, "PCdE_vs_Z");
 
@@ -2495,7 +2592,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         {
           TVector3 pc = a1c0_wirePos(apwire_bm, phi_use, false);
           TVector3 vtx0 = vertexFrom(si_point, pc);
-          if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
+          if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance))
             return;
           double pcz = dither ? rand.Gaus(pc.Z(), dither_sigma) : pc.Z();
           TVector3 vtx = vertexFrom(si_point, TVector3(pc.X(), pc.Y(), pcz));
@@ -2543,7 +2640,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
             {
               double pcz_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi()).Z();
-              double theta_ref = (sx3event.pos - TVector3(0, 0, vtx_ref.Z())).Theta() * 180. / M_PI;
+              double theta_ref = (sx3event.pos - beamAxisPoint(vtx_ref.Z())).Theta() * 180. / M_PI;
               plotter->Fill2D("Benchmark_SX3_PCZ_A1C0_minus_ref_vs_theta", 180, 0, 180, 400, -200, 200, theta_ref, pcz_a1c0 - pcz_ref, "Benchmark_SX3_ref");
               plotter->Fill2D("Benchmark_PCZ_A1C0_minus_ref_vs_theta", 180, 0, 180, 400, -200, 200, theta_ref, pcz_a1c0 - pcz_ref, "Benchmark_AnodeOnly");
 
@@ -2618,7 +2715,6 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
                   plotter->Fill2D("AnodeEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 40000, 1 - TMath::Abs(pcz_ref - zc) / half, pcevent.Energy1);
                   plotter->Fill2D("CathodeEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 40000, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2);
                   plotter->Fill2D("FracEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 10, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2 / pcevent.Energy1);
-                  plotter->Fill2D("SumEnergy_vs_CellSX3", 120, 0, 1.2, 800, 0, 10, TMath::Abs(pcz_ref - zc) / half, (pcevent.Energy1 + pcevent.Energy2) / 2);
 
                   if (half > 0.0)
                   {
@@ -2695,7 +2791,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
             {
               TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, sx3event.pos.Phi());
               TVector3 vtx0 = vertexFrom(sx3event.pos, pc);
-              if (vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6)
+              if (vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance)
               {
                 fillSuite("A1C1asA1C0", pc.Z(), vtx0, "A1C1True_SX3");
                 plotter->Fill2D("Benchmark_SX3_PCZ_A1C1asA1C0_vs_sx3pczguess", 400, -200, 200, 400, -200, 200, pczguess, pc.Z(), "A1C1True_SX3");
@@ -2711,7 +2807,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events,
                           const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters)
 {
-  static TRandom3 rand(0);
+  TRandom3 &rand = anasenRandom;
 
   // --- GENUINE A1C0 events (QQQ twin):
   if (BenchMark && aClusters.size() == 1 && cClusters.size() == 0)
@@ -2722,7 +2818,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
     double anodeTS = std::get<3>(aPw);
     for (const auto &qqqevent : QQQ_Events)
     {
-      bool timecut = (qqqevent.Time1 - anodeTS < 150);
+      bool timecut = siPcCoincident(qqqevent.Time1, anodeTS);
       double smeared_phi = qqqevent.pos.Phi() + rand.Uniform(-qqq_wedge_pitch / 2.0, qqq_wedge_pitch / 2.0);
       TVector3 pc = a1c0_wirePos(apwire_bm, smeared_phi, true);
 
@@ -2733,13 +2829,13 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       double smeared_rho = qqqevent.pos.Perp() + rand.Uniform(-qqq_ring_pitch / 2.0, qqq_ring_pitch / 2.0);
       TVector3 smeared_qqq(smeared_rho * TMath::Cos(smeared_phi), smeared_rho * TMath::Sin(smeared_phi), qqqevent.pos.Z());
       // A1C0 hybrid z (shared with the SX3 twin block + miscHistograms_oneWire).
-      TVector3 pc_hybrid = a1c0_hybrid_pcz(apwire_bm, smeared_phi, true, dither_sigma_c0 / 2.0, rand);
+      TVector3 pc_hybrid = a1c0_hybrid_pcz(apwire_bm, smeared_phi, true, dither_sigma, rand);
       TVector3 vtx0 = beamVertex(qqqevent.pos, pc - qqqevent.pos);
       TVector3 vtx1 = beamVertex(smeared_qqq, pc_hybrid - smeared_qqq);
 
-      if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
+      if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance))
         continue;
-      double qqqTheta = (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta();
+      double qqqTheta = (qqqevent.pos - beamAxisPoint(source_vertex)).Theta();
       double pcz_guess_37 = 37. / TMath::Tan(qqqTheta) + source_vertex;
       plotter->Fill1D("Benchmark_QQQ_VertexZ_trueA1C0", 800, -400, 400, vtx0.Z(), "A1C0True_QQQ");
       plotter->Fill1D("Benchmark_QQQ_VertexZ_trueA1C0_Hybrid", 800, -400, 400, vtx1.Z(), "A1C0True_QQQ");
@@ -2754,9 +2850,6 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
   for (const auto &pcevent : PC_Events)
   {
 
-    plotter->Fill2D("pcdEACAvg_vs_anodechannel", 24, 0, 23, 800, 0, 20000, pcevent.Anodech, (pcevent.Energy1 + pcevent.Energy2) / 2);
-    plotter->Fill2D("pcdEACAvg_vs_cathodechannel", 24, 0, 23, 800, 0, 20000, pcevent.Cathodech, (pcevent.Energy1 + pcevent.Energy2) / 2);
-
     for (const auto &qqqevent : QQQ_Events)
     {
       plotter->Fill1D("dt_pcA_qqqR", 640, -2000, 2000, qqqevent.Time1 - pcevent.Time1, "Timing");
@@ -2764,16 +2857,18 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       plotter->Fill1D("dt_pcC_qqqW", 640, -2000, 2000, qqqevent.Time2 - pcevent.Time2, "Timing");
       plotter->Fill2D("phiPC_vs_phiQQQ", 100, -200, 200, 100, -200, 200, qqqevent.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
 
-      double qqqTheta = (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta();
+      double qqqTheta = (qqqevent.pos - beamAxisPoint(source_vertex)).Theta();
       double sinTheta = TMath::Sin(qqqTheta);
 
       TVector3 x2(pcevent.pos);
       TVector3 x1(qqqevent.pos);
       TVector3 v = x2 - x1;
-      double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
-      TVector3 r_rhoMin = x1 + t_minimum * v;
+      // beamVertex() instead of an inline projection: the hand-rolled version pinned the
+      // beam axis at (0,0), silently ignoring BEAM_AXIS_X/Y, and had no guard for a
+      // purely longitudinal direction.
+      TVector3 r_rhoMin = beamVertex(x1, v);
 
-      bool timecut = (qqqevent.Time1 - pcevent.Time1 < 150);
+      bool timecut = siPcCoincident(qqqevent.Time1, pcevent.Time1);
       bool lowercut_cath = pcevent.Energy2 * sinTheta < 1 && (qqqevent.Energy2 < 5.0 || qqqevent.Energy1 < 5.0);
 
       bool phicut = TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
@@ -2787,7 +2882,9 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       if (timecut)
       {
         plotter->Fill2D("dE_E_AnodeQQQR", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
-        plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy2, pcevent.Energy2, "PID_dE_E");
+        // ring energy (Energy1), matching every _a1c*/_a2c* subset below -- this was the
+        // only cathode-PID plot on the wedge energy, so its subsets never summed to it
+        plotter->Fill2D("dE_E_CathodeQQQR", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
         if (pcevent.multi1 == 1 && pcevent.multi2 == 2)
         {
           plotter->Fill2D("dE_E_AnodeQQQR_a1c2", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
@@ -2798,20 +2895,23 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           plotter->Fill2D("dE_E_AnodeQQQR_a2c1", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
           plotter->Fill2D("dE_E_CathodeQQQR_a2c1", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
         }
+        // These three topologies were copied from the SX3 twin without renaming, so they
+        // filled the "...sx3B..." histograms with QQQ data. HistPlotter keys its object
+        // map on the name alone, so both detectors landed in one histogram each.
         if (pcevent.multi1 == 1 && pcevent.multi2 == 1)
         {
-          plotter->Fill2D("dE_E_Anodesx3B_a1c1", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
-          plotter->Fill2D("dE_E_Cathodesx3B_a1c1", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
+          plotter->Fill2D("dE_E_AnodeQQQR_a1c1", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
+          plotter->Fill2D("dE_E_CathodeQQQR_a1c1", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
         }
         if (pcevent.multi1 == 1 && pcevent.multi2 == 0)
         {
-          plotter->Fill2D("dE_E_Anodesx3B_a1c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
-          plotter->Fill2D("dE_E_Cathodesx3B_a1c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
+          plotter->Fill2D("dE_E_AnodeQQQR_a1c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
+          plotter->Fill2D("dE_E_CathodeQQQR_a1c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
         }
         if (pcevent.multi1 == 2 && pcevent.multi2 == 0)
         {
-          plotter->Fill2D("dE_E_Anodesx3B_a2c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
-          plotter->Fill2D("dE_E_Cathodesx3B_a2c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
+          plotter->Fill2D("dE_E_AnodeQQQR_a2c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
+          plotter->Fill2D("dE_E_CathodeQQQR_a2c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
         }
         if (phicut)
         {
@@ -2838,14 +2938,16 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
         plotter->Fill1D("phiQQQ_minus_phiPC_TimeCut_QQQ" + std::to_string(qqqevent.ch1 / 16), 180, -180, 180, qqqevent.pos.DeltaPhi(pcevent.pos) * 180 / M_PI, "Kinematics_Angles");
 
-        plotter->Fill2D("Etot2_theta_AnodeQQQR", 75, 0, 90, 300, 0, 15, qqqTheta * 180 / M_PI, qqqevent.Energy1 + pcevent.Energy1 * anode_gain * sinTheta, "Kinematics_Angles");
+        // (Etot2_theta_AnodeQQQR removed with the anode_gain constant: it was the only
+        // consumer of a second, hardcoded anode ADC->MeV factor that competed with the
+        // per-wire pcEnergySlope[] table. Use the PC_Events_calibrated path instead.)
 
         plotter->Fill2D("dE_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, qqqTheta * 180 / M_PI, pcevent.Energy2, "Kinematics_Angles");
         plotter->Fill2D("dE2_theta_CathodeQQQR", 75, 0, 90, 800, 0, 10000, qqqTheta * 180 / M_PI, pcevent.Energy2 * sinTheta, "Kinematics_Angles");
         plotter->Fill2D("dE2_theta_CathodeQQQR_zoomin", 60, 0, 30, 800, 0, 3000, qqqTheta * 180 / M_PI, pcevent.Energy2 * sinTheta, "Kinematics_Angles");
 
-        plotter->Fill2D("dE_phi_AnodeQQQR", 90, -180, 180, 800, 0, 40000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy1, "Kinematics_Angles");
-        plotter->Fill2D("dE_phi_CathodeQQQR", 90, -180, 180, 800, 0, 40000, (qqqevent.pos - TVector3(0, 0, source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2, "Kinematics_Angles");
+        plotter->Fill2D("dE_phi_AnodeQQQR", 90, -180, 180, 800, 0, 40000, (qqqevent.pos - beamAxisPoint(source_vertex)).Phi() * 180 / M_PI, pcevent.Energy1, "Kinematics_Angles");
+        plotter->Fill2D("dE_phi_CathodeQQQR", 90, -180, 180, 800, 0, 40000, (qqqevent.pos - beamAxisPoint(source_vertex)).Phi() * 180 / M_PI, pcevent.Energy2, "Kinematics_Angles");
 
         plotter->Fill1D("PCZ", 800, -200, 200, pcevent.pos.Z(), "PCZ_Recon");
 
@@ -2862,11 +2964,13 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           double pcz_fix = a1c2_zfix(pcevent.pos.Z());
           TVector3 x2f(pcevent.pos.X(), pcevent.pos.Y(), pcz_fix);
           TVector3 v = x2f - x1;
-          double t_minimum = -1.0 * (x1.X() * v.X() + x1.Y() * v.Y()) / (v.X() * v.X() + v.Y() * v.Y());
-          TVector3 r_rhoMin_fix = x1 + t_minimum * v;
+          // beamVertex() instead of an inline projection: the hand-rolled version pinned the
+          // 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);
 
-          double sinTheta_customV = TMath::Sin((qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta());
-          plotter->Fill2D("dE3_E_CathodeQQQW_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
+          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");
           plotter->Fill2D("dE3_E_AnodeQQQR_A1C2_TC1_PC" + std::to_string(phicut), 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
 
           plotter->Fill1D("VertexRecon_pczfix_qqq", 800, -300, 300, r_rhoMin_fix.Z(), "Vertex_Reconstruction");
@@ -2890,24 +2994,21 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           plotter->Fill2D("pcdEA_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy1, "PCdE_vs_Z");
           plotter->Fill2D("pcdEC_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, pcevent.Energy2, "PCdE_vs_Z");
 
-          plotter->Fill2D("pcdEACSum_vs_pczfix", 600, -200, 200, 800, 0, 20000, pcz_fix, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z");
-
           double path_length = pathLengthCm(qqqevent.pos, r_rhoMin_fix);
           double qqqEfix = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, qqqevent.Energy1, path_length);
           double qqqEfix_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, qqqevent.Energy1, path_length);
 
-          plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqevent.Energy1, "Kinematics_Angles");
+          plotter->Fill2D("E_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut), 180, 0, 180, 600, 0, 15, (qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqevent.Energy1, "Kinematics_Angles");
           if (lowercut_cath)
-            plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix_p, "Kinematics_Angles");
+            plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 180, 0, 180, 600, 0, 15, (qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix_p, "Kinematics_Angles");
           else
           {
             std::string zcut = "_" + std::to_string((TMath::Abs(r_rhoMin_fix.Z()) < 180));
-            plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath) + zcut, 180, 0, 180, 600, 0, 15, (qqqevent.pos - TVector3(0, 0, r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix, "Kinematics_Angles");
+            plotter->Fill2D("Ef_thetaf_AnodeQQQR_TC1_PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath) + zcut, 180, 0, 180, 600, 0, 15, (qqqevent.pos - beamAxisPoint(r_rhoMin_fix.Z())).Theta() * 180 / M_PI, qqqEfix, "Kinematics_Angles");
           }
 
           plotter->Fill2D("dE3_Ef_AnodeQQQR_TC1" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 40000, qqqEfix, pcevent.Energy1 * sinTheta_customV, "PID_dE_E");
           plotter->Fill2D("dE3_Ef_CathodeQQQR_TC1PC" + std::to_string(phicut) + "_pidlow" + std::to_string(lowercut_cath), 600, 0, 15, 800, 0, 10000, qqqEfix, pcevent.Energy2 * sinTheta_customV, "PID_dE_E");
-          plotter->Fill2D("pcdEACAvg_A1C2_vs_qqqE", 400, 0, 10, 1600, 0, 40000, qqqevent.Energy1, (pcevent.Energy1 + pcevent.Energy2) / 2);
         }
 
         // plotter->Fill2D("pcdEA_vs_qqqpczguess_A" + std::to_string(pcevent.multi1) + "C" + std::to_string(pcevent.multi2), 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy1, "PCdE_vs_Z");
@@ -2918,11 +3019,6 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         plotter->Fill2D("pcdEC_vs_qqqpczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy2, "PCdE_vs_Z");
         plotter->Fill2D("pcdEC_vs_pczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, pcevent.Energy2, "PCdE_vs_Z");
 
-        plotter->Fill2D("pcdEACAvg_vs_qqqpczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z");
-        plotter->Fill2D("pcdEACAvg_vs_pczguess", 600, -200, 200, 800, 0, 20000, pcz_guess_37, (pcevent.Energy1 + pcevent.Energy2) / 2, "PCdE_vs_Z");
-
-        plotter->Fill2D("pcdEACAvg_vs_qqqE", 400, 0, 10, 1600, 0, 40000, qqqevent.Energy1, (pcevent.Energy1 + pcevent.Energy2) / 2);
-
         // plotter->Fill2D("pcdEC_vs_pczfix" + std::to_string(pcevent.multi1) + "A" + std::to_string(pcevent.multi1) + "C", 800, 0, 20000, 600, -200, 200, pcevent.Energy2, pcz_fix, "PCdE_vs_Z");
 
         //-----------------------Benchmarking Method for Source Runs (QQQ)------------------------//
@@ -2999,9 +3095,9 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           {
             TVector3 pc = a1c0_wirePos(apwire_bm, phi_use, true);
             TVector3 vtx0 = vertexFrom(si_point, pc);
-            if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6))
+            if (!(vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance))
               return;
-            double pcz = dither ? rand.Gaus(pc.Z(), dither_sigma_c0 / 2.0) : pc.Z();
+            double pcz = dither ? rand.Gaus(pc.Z(), dither_sigma) : pc.Z();
             TVector3 vtx = vertexFrom(si_point, TVector3(pc.X(), pc.Y(), pcz));
             fillSuite(tag, pcz, vtx, benchBranch);
             fillVsRef(tag, pcz, vtx, pcz_ref, vtx_ref);
@@ -3048,7 +3144,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 
               {
                 double pcz_a1c0 = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi()).Z();
-                double theta_ref = (qqqevent.pos - TVector3(0, 0, vtx_ref.Z())).Theta() * 180. / M_PI;
+                double theta_ref = (qqqevent.pos - beamAxisPoint(vtx_ref.Z())).Theta() * 180. / M_PI;
                 plotter->Fill2D("Benchmark_QQQ_PCZ_A1C0_minus_ref_vs_theta", 180, 0, 180, 400, -200, 200, theta_ref, pcz_a1c0 - pcz_ref, "Benchmark_QQQ_ref");
                 plotter->Fill2D("Benchmark_PCZ_A1C0_minus_ref_vs_theta", 180, 0, 180, 400, -200, 200, theta_ref, pcz_a1c0 - pcz_ref, "Benchmark_AnodeOnly");
 
@@ -3122,7 +3218,6 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
                     plotter->Fill2D("AnodeEnergy_vs_CellQQQ", 120, 0, 1.2, 800, 0, 40000, 1 - TMath::Abs(pcz_ref - zc) / half, pcevent.Energy1);
                     plotter->Fill2D("CathodeEnergy_vs_CellQQQ", 120, 0, 1.2, 800, 0, 40000, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2);
                     plotter->Fill2D("FracEnergy_vs_CellQQQ", 120, 0, 1.2, 1200, 0, 20, TMath::Abs(pcz_ref - zc) / half, pcevent.Energy2 / pcevent.Energy1);
-                    plotter->Fill2D("SumEnergy_vs_CellQQQ", 120, 0, 1.2, 800, 0, 40000, TMath::Abs(pcz_ref - zc) / half, (pcevent.Energy2 + pcevent.Energy1) / 2);
 
                     if (half > 0.0)
                     {
@@ -3197,7 +3292,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
             {
               TVector3 pc = pwinstance.getClosestWirePosAtWirePhi(apwire_bm, qqqevent.pos.Phi());
               TVector3 vtx0 = vertexFrom(qqqevent.pos, pc);
-              if (vtx0.Perp() <= 6.0 && vtx0.Z() >= -173.6)
+              if (vtx0.Perp() <= 6.0 && vtx0.Z() >= z_entrance)
               {
                 fillSuite("A1C1asA1C0", pc.Z(), vtx0, "A1C1True_QQQ");
                 plotter->Fill2D("Benchmark_QQQ_PCZ_A1C1asA1C0_vs_qqqpczguess", 400, -200, 200, 400, -200, 200, pcz_guess_int, pc.Z(), "A1C1True_QQQ");
@@ -3206,7 +3301,7 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
           }
         }
         double qqqrho = qqqevent.pos.Perp();
-        double qqqz = (qqqevent.pos - TVector3(0, 0, source_vertex)).Z();
+        double qqqz = (qqqevent.pos - beamAxisPoint(source_vertex)).Z();
         double tan_theta = qqqrho / qqqz;
         double pcz_guess_int2 = z_to_crossover_rho(pcevent.pos.Z()) / tan_theta + source_vertex;
         plotter->Fill2D("pczguess_vs_pc_int2", 180, 0, 200, 150, 0, 200, pcz_guess_int2, pcevent.pos.Z(), "PCZ_Recon");
@@ -3217,8 +3312,8 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
         plotter->Fill2D("pczguess_vs_pc_int3", 180, 0, 200, 150, 0, 200, pcz_guess_int3, pcevent.pos.Z(), "PCZ_Recon");
 
         double pcz_guess = pcz_guess_int;
-        plotter->Fill2D("pctheta_vs_qqqtheta_sv", 180, -200, 200, 180, -200, 200, qqqTheta * 180 / M_PI, (pcevent.pos - TVector3(0, 0, source_vertex)).Theta() * 180 / M_PI, "Kinematics_Angles");
-        plotter->Fill2D("pctheta_vs_qqqtheta_rmz", 180, -200, 200, 180, -200, 200, (qqqevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, (pcevent.pos - TVector3(0, 0, r_rhoMin.Z())).Theta() * 180 / M_PI, "Kinematics_Angles");
+        plotter->Fill2D("pctheta_vs_qqqtheta_sv", 180, -200, 200, 180, -200, 200, qqqTheta * 180 / M_PI, (pcevent.pos - beamAxisPoint(source_vertex)).Theta() * 180 / M_PI, "Kinematics_Angles");
+        plotter->Fill2D("pctheta_vs_qqqtheta_rmz", 180, -200, 200, 180, -200, 200, (qqqevent.pos - beamAxisPoint(r_rhoMin.Z())).Theta() * 180 / M_PI, (pcevent.pos - beamAxisPoint(r_rhoMin.Z())).Theta() * 180 / M_PI, "Kinematics_Angles");
         plotter->Fill2D("pctheta_vs_qqqtheta_rm", 180, -200, 200, 180, -200, 200, (qqqevent.pos - r_rhoMin).Theta() * 180 / M_PI, (pcevent.pos - r_rhoMin).Theta() * 180 / M_PI, "Kinematics_Angles");
         plotter->Fill2D("pczguess_vs_pc_phi=" + std::to_string(qqqevent.pos.Phi() * 180. / M_PI), 300, 0, 200, 150, 0, 200, pcz_guess, pcevent.pos.Z(), "Z_Reconstruction");
       }
@@ -3291,7 +3386,7 @@ void TrackRecon::OldAnalysis()
   // position estimate) are untouched.
 
   if (anodeHits.size() > 0 && cathodeHits.size() > 0)
-    plotter->Fill2D("AHits_vs_CHits", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
+    plotter->Fill2D("AHits_vs_CHits", 13, -0.5, 12.5, 7, -0.5, 6.5, anodeHits.size(), cathodeHits.size(), "hRawPC");
 
   // make another plot with nearest neighbour constraint
   bool hasNeighbourAnodes = false;
@@ -3330,12 +3425,12 @@ void TrackRecon::OldAnalysis()
   if (anodeHits.size() > 0 && cathodeHits.size() > 0)
   {
 #ifdef RAW_HISTOS
-    plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
-    plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
+    plotter->Fill2D("AHits_vs_CHits_NA" + std::to_string(hasNeighbourAnodes), 13, -0.5, 12.5, 7, -0.5, 6.5, anodeHits.size(), cathodeHits.size(), "hRawPC");
+    plotter->Fill2D("AHits_vs_CHits_NC" + std::to_string(hasNeighbourCathodes), 13, -0.5, 12.5, 7, -0.5, 6.5, anodeHits.size(), cathodeHits.size(), "hRawPC");
 
     if (hasNeighbourAnodes && hasNeighbourCathodes)
     {
-      plotter->Fill2D("AHits_vs_CHits_NN", 12, 0, 11, 6, 0, 5, anodeHits.size(), cathodeHits.size(), "hRawPC");
+      plotter->Fill2D("AHits_vs_CHits_NN", 13, -0.5, 12.5, 7, -0.5, 6.5, anodeHits.size(), cathodeHits.size(), "hRawPC");
     }
 #endif
   }
@@ -3359,14 +3454,14 @@ void TrackRecon::OldAnalysis()
       plotter->Fill2D("CMax_over_Anode_vs_CathodeID", 24, 0, 24, 200, 0, 2.0, cIDMax, ratio, folder);
     }
   }
-  plotter->Fill1D("Correlated_Cathode_MaxAnode", 6, 0, 5, corrcatMax.size(), "hGMPC");
-  plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 6, 0, 5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
-  plotter->Fill1D("AnodeHits", 12, 0, 11, anodeHits.size(), "hGMPC");
-  plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 12, 0, 11, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
+  plotter->Fill1D("Correlated_Cathode_MaxAnode", 7, -0.5, 6.5, corrcatMax.size(), "hGMPC");
+  plotter->Fill2D("Correlated_Cathode_VS_MaxAnodeEnergy", 7, -0.5, 6.5, 2000, 0, 30000, corrcatMax.size(), aEMax, "hGMPC");
+  plotter->Fill1D("AnodeHits", 13, -0.5, 12.5, anodeHits.size(), "hGMPC");
+  plotter->Fill2D("AnodeMaxE_vs_AnodeHits", 13, -0.5, 12.5, 2000, 0, 30000, anodeHits.size(), aEMax, "hGMPC");
 
   if (anodeHits.size() < 1)
   {
-    plotter->Fill1D("NoAnodeHits_CathodeHits", 6, 0, 5, cathodeHits.size(), "hGMPC");
+    plotter->Fill1D("NoAnodeHits_CathodeHits", 7, -0.5, 6.5, cathodeHits.size(), "hGMPC");
   }
 
   for (const auto &cwevent : cWireEvents)
@@ -3382,7 +3477,7 @@ void TrackRecon::OldAnalysis()
 void miscHistograms_oneWire(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters)
 {
   // consider the 'proton-like' QQQ branch seen in a,p data
-  static TRandom3 rand(0); // seeded once (random seed via TUUID), not per call
+  TRandom3 &rand = anasenRandom;
   double initial_energy = 6.89;
 
   Kinematics apkin_a(mass_1H, mass_4He, mass_4He, mass_1H, initial_energy / mass_1H);
@@ -3402,10 +3497,10 @@ void miscHistograms_oneWire(HistPlotter *plotter, const std::vector<Event> &QQQ_
       int a_number = acluster.size();
       TVector3 pc_closest = pwinstance.getClosestWirePosAtWirePhi(apwire, qqqevent.pos.Phi());
       plotter->Fill1D("dt_anode_interp_qqq", 800, -2000, 2000, qqqevent.Time1 - apTSMaxE, "ainterp_noc");
-      if (qqqevent.Time1 - apTSMaxE < 150)
+      if (siPcCoincident(qqqevent.Time1, apTSMaxE))
       {
         bool phicut = TMath::Abs(qqqevent.pos.DeltaPhi(pc_closest)) <= TMath::Pi() / 4.0;
-        TVector3 pc_hybrid = a1c0_hybrid_pcz(apwire, qqqevent.pos.Phi(), true, dither_sigma_c0 / 2.0, rand);
+        TVector3 pc_hybrid = a1c0_hybrid_pcz(apwire, qqqevent.pos.Phi(), true, dither_sigma, rand);
         TVector3 r_rhoMin_fix = beamVertex(qqqevent.pos, pc_hybrid - qqqevent.pos);
 
         double theta_q = (qqqevent.pos - r_rhoMin_fix).Theta();
@@ -3413,7 +3508,7 @@ void miscHistograms_oneWire(HistPlotter *plotter, const std::vector<Event> &QQQ_
 
         if (beamPerp(r_rhoMin_fix) > 6.0)
           continue;
-        if (r_rhoMin_fix.Z() < -173.6 || r_rhoMin_fix.Z() > 100)
+        if (r_rhoMin_fix.Z() < z_entrance || r_rhoMin_fix.Z() > 100)
           continue;
         if (!phicut)
           continue;
@@ -3476,7 +3571,7 @@ void miscHistograms_oneWire(HistPlotter *plotter, const std::vector<Event> &QQQ_
 void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events)
 {
   // consider the 'proton-like' QQQ branch seen in a,p data
-  static TRandom3 rand(0);
+  TRandom3 &rand = anasenRandom;
   double initial_energy = 6.89;
 
   for (const auto &qqqevent : QQQ_Events)
@@ -3576,7 +3671,7 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       double sinTheta_customV = TMath::Sin(theta_q);
       // if(beamPerp(r_rhoMin_fix)>6) continue;
       bool cathode_alpha_select = (pcevent.Energy2 > 1400);
-      if (vertex_z < -173.6 || vertex_z > 100)
+      if (vertex_z < z_entrance || vertex_z > 100)
         continue;
 
       double beam_path_length_q = TMath::Abs(vertex_z - z_entrance) * 0.1;
@@ -3594,7 +3689,7 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
       double pa_dist_mm = (qqqevent.pos - r_rhoMin_fix).Mag();
       double pa_pathfraction = (pa_dist_mm > 0.0) ? pa_dl_cm * 10.0 / pa_dist_mm : 0.0;
       double pcz_guess_int = z_to_crossover_rho(pcevent.pos.Z()) /
-                                 TMath::Tan((qqqevent.pos - TVector3(0, 0, source_vertex)).Theta()) +
+                                 TMath::Tan((qqqevent.pos - beamAxisPoint(source_vertex)).Theta()) +
                              source_vertex;
 
       // Calibrated anode energy, same lookup reaction_ax_core uses for its dEgasCalib plots.
@@ -3722,7 +3817,7 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
 void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events)
 {
   // consider the 'proton-like' QQQ branch seen in a,p data
-  static TRandom3 rand(0); // seeded once (random seed via TUUID), not per call
+  TRandom3 &rand = anasenRandom;
   double initial_energy = 6.89;
 
   for (const auto &sx3event : SX3_Events)
@@ -3751,7 +3846,7 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector<Event> &QQ
 
       if (beamPerp(r_rhoMin_fix) > 10.0)
         continue;
-      if (vertex_z < -173.6 || vertex_z > 100)
+      if (vertex_z < z_entrance || vertex_z > 100)
         continue; // same beam-region acceptance as the QQQ branch
       double theta_s = (sx3event.pos - r_rhoMin_fix).Theta();
       double sinTheta_customV = TMath::Sin(theta_s);
@@ -3888,7 +3983,7 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector<Event> &QQ
 
       if (beamPerp(r_rhoMin) > 10.0)
         continue;
-      if (vertex_z < -173.6 || vertex_z > 100)
+      if (vertex_z < z_entrance || vertex_z > 100)
         continue; // same beam-region acceptance as the A1C2/A1C1 loops above
 
       double theta_s = (sx3event.pos - r_rhoMin).Theta();
@@ -3989,7 +4084,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
                              const AAEjectileMasses &ej_m, const std::string &globaltag)
 {
   const std::string sfx = "_" + det + globaltag;
-  static TRandom3 rand(0);
+  TRandom3 &rand = anasenRandom;
   for (const auto &sievent : Si_Events)
   {
     if (sievent.Energy1 < si_ecut)
@@ -4105,7 +4200,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, 100, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel);
+        plotter->Fill2D(rx + "_Ex_vs_theta" + ejtag + sfx, 180, 0, 180, 800, -20, 20, theta * 180 / M_PI, Ex, pmlabel);
         plotter->Fill1D(rx + "_VertexReconZ" + sfx, 800, -400, 400, vertex_z, pmlabel);
 
         forEachTier(topo1, topo2, methodGroup, plot_with_tag);
@@ -4232,7 +4327,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
       if (anodeCh_a0 < 0 || anodeCh_a0 >= 24)
         anodeCh_a0 = -1;
 
-      double pcz_a0 = isA2C0 ? pc.Z() : rand.Gaus(pc.Z(), dither_sigma_c0 / 2.0);
+      double pcz_a0 = isA2C0 ? pc.Z() : rand.Gaus(pc.Z(), dither_sigma);
       reconstructAndFill(pcz_a0, pc, apSumE, -1.0, anodeE_MeV_a0, -1.0, a0tag, "", anodeCh_a0);
     }
   }
diff --git a/eloss_calculations/Eloss.py b/eloss_calculations/Eloss.py
index f24f307..f4c3bda 100644
--- a/eloss_calculations/Eloss.py
+++ b/eloss_calculations/Eloss.py
@@ -1,5 +1,9 @@
 import pycatima as catima
 import numpy as np
+import os
+
+DEDX_SCALE = float(os.environ.get("DEDX_SCALE", "1.0"))
+print(f"Using dEdX scale factor: {DEDX_SCALE}")
 
 # --- 1. Constants ---
 P_TORR = 250
@@ -32,7 +36,7 @@ fine_step_cm = 0.03     # row spacing near the Bragg peak
 fine_zone_frac = 0.085   # fraction of the *total* range treated as "near the peak"
 
 def generate_lookup(z, mass_u, e_start_mev, label):
-    filename = f"{label}_lookup_{e_start_mev}MeV_{P_TORR}torr_{P_CO2}pc.dat"
+    filename = f"/home/vsitaraman/ANASEN_analysis/eloss_calculations/{label}_lookup_{e_start_mev}MeV_{P_TORR}torr_{P_CO2}pc.dat"
     header = f"Energy(MeV) \tmg/cm2 \tcm\nStarting Energy: {e_start_mev} MeV"
 
     # Pass 1: integrate at full precision just to find the total range (needed
@@ -78,7 +82,7 @@ def generate_lookup(z, mass_u, e_start_mev, label):
         projectile.T(e_u)
         # dedx returns MeV / (g/cm2)
         if(mass_u >=10.0):
-            loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2 * 0.89
+            loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2 * DEDX_SCALE
         else:
             loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2
 
diff --git a/run_27Al.sh b/run_27Al.sh
index 6bc8d4a..a7aac53 100644
--- a/run_27Al.sh
+++ b/run_27Al.sh
@@ -6,14 +6,15 @@ export OUT_DIR="Output_27Al"
 export reactiondata=1
 export CO2percent=3
 export pressure_in_torr=250
-export CATHODE_GAIN=4.0
+export CATHODE_GAIN=3.0
 export source_vertex=-200.0
-export BEAM_AXIS_X=-15
-export BEAM_AXIS_Y=-5
+export DEDX_SCALE=0.90
+# export BEAM_AXIS_X=-15
+# export BEAM_AXIS_Y=-5
 export CUTLIST=cuts_list.txt
 
 # Clean up previous runs
-rm -f Output_27Al/*.root
+rm -f Output_27Al_$DEDX/*.root
 
 echo "Pre-compiling TrackRecon.C safely on a single core..."
 root -q -l -b -e '.L TrackRecon.C++O'
@@ -22,9 +23,9 @@ process_run() {
     local wrun=$(printf "%03d" "$1")
     local prefix="${PREFIX:-Run_}"
     local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root"
-    local out="Output_27Al/results_run${wrun}.root"
+    local out="Output_27Al_$DEDX_SCALE/results_run${wrun}.root"
 
-    mkdir -p Output_27Al
+    mkdir -p Output_27Al_$DEDX_SCALE
 
     root -q -l -b -x "$infile" \
          -e "tree->Process(\"TrackRecon.C+\", \"${out}\")" > /dev/null 2>&1
@@ -38,9 +39,11 @@ process_run() {
 
 export -f process_run
 
+echo "Running Eloss.py with a scaling parameter of $DEDX_SCALE"
+python3 eloss_calculations/Eloss.py
 echo "Starting parallel processing..."
 
-# time parallel --bar -j 8 process_run ::: {24..41} 
+time parallel --bar -j 8 process_run ::: {24..41} 
 time parallel --bar -j 3 process_run ::: 44 45 46
 time parallel --bar -j 8 process_run ::: {50..59}
 # time parallel --bar -j 4 process_run ::: 62 63 66 67 68
@@ -49,7 +52,7 @@ time parallel --bar -j 8 process_run ::: {50..59}
 # time parallel --bar -j 4 process_run ::: {78..89}
 
 echo "Merging files..."
-hadd -k -j 4 Output_27Al/output_27Al.root Output_27Al/results_run*.root
+hadd -k -j 4 Output_27Al_$DEDX_SCALE/output_27Al.root Output_27Al_$DEDX_SCALE/results_run*.root
 
 # rootbrowse Output_27Al/output_27Al.root
 
@@ -68,4 +71,5 @@ unset A1C1_Z_OFF_SX3
 unset BEAM_AXIS_X
 unset BEAM_AXIS_Y
 unset CUTLIST
+unset DEDX_SCALE
 echo "Script execution finished."
\ No newline at end of file