modified: .gitignore

modified:   TrackRecon.C
	new file:   pc_energy_calibration.dat
	modified:   pccal/fit_pc_energy_calibration.C
	modified:   run_27Al.sh
	modified:   run_tr.sh

.gitignore

@@ -36,3 +36,4 @@ Armory/EventBuilder
 EventBuilder
 anasen_fem/anode_to_cathode_1d_17.43.csv
 CLAUDE.md
+pc_calib_raw/

TrackRecon.C

@@ -44,15 +44,16 @@ bool process_alpha_proton_scattering = false,
      doMiscHistograms = true,
      doPCSX3ClusterAnalysis = true,
      doPCQQQClusterAnalysis = true,
-     doOldAnalysis = false,
+     doOldAnalysis = true,
      do27AlapAnalysis = false,
      BenchMark = true,
      onwire_analysis = true,
      diagnostic_eplots = true,
-     diagnostic_tplots = false,
+     diagnostic_tplots = true,
      reactiondata = false,
      doPCEnergyCalibration = false,
-     ta_foil_run = false;
+     ta_foil_run = false,
+     source_run = false;
 
 // --- Geometry, Calibration, & Model Variables ---
 double source_vertex = 53.0,
@@ -153,7 +154,17 @@ static std::vector<int> a1c1_dead_cathode_17F = {};    // 0,13,15 can be recover
 static std::vector<int> a1c1_dead_anode_27Al = {0, 12, 19};
 static std::vector<int> a1c1_dead_cathode_27Al = {13};
 
-std::vector<std::pair<double, double>> pcCalibData[48];
+// Calibration points are streamed straight to disk as they're generated
+// (opened in Begin(), written by pcCalibWritePoint(), closed in Terminate())
+// instead of buffered in memory for the whole run -- a full run can produce
+// many millions of points, and holding them all in a std::vector until
+// Terminate() grows unbounded and can exhaust memory on long runs.
+std::ofstream pcCalibOutFile;
+inline void pcCalibWritePoint(int wire, double adc, double dE_gas)
+{
+  if (pcCalibOutFile.is_open())
+    pcCalibOutFile << wire << " " << adc << " " << dE_gas << "\n";
+}
 
 std::vector<int> *a1c1_dead_anode = &a1c1_dead_anode_17F; // active set, chosen in Begin()
 std::vector<int> *a1c1_dead_cathode = &a1c1_dead_cathode_17F;
@@ -507,20 +518,27 @@ void TrackRecon::Begin(TTree * /*tree*/)
     std::cout << "Analyzing dataset as reactiondata" << std::endl;
   }
 
-  if (getenv("RUN_NUMBER"))
+  // Run classification is by OUT_DIR, not RUN_NUMBER: run numbers collide
+  // across datasets/blocks (e.g. 17F's alpha+gas block and 27Al's proton
+  // block both use runs 18-21), so a RUN_NUMBER-keyed lookup misclassifies
+  // whichever dataset didn't originally populate kTaFoilRuns.
+  std::string outdir = getenv("OUT_DIR") ? getenv("OUT_DIR") : "";
+  ta_foil_run = (outdir == "Output_p");
+  source_run = (outdir == "Output_a");
+  if (ta_foil_run && getenv("RUN_NUMBER"))
   {
     int run_number = std::atoi(getenv("RUN_NUMBER"));
     for (const auto &r : kTaFoilRuns)
     {
       if (r.run == run_number)
       {
-        ta_foil_run = true;
         ta_foil_z_mm = r.z_mm;
         break;
       }
     }
   }
-  std::cout << "Ta foil: " << (ta_foil_run ? ("present, run in proton-scattering campaign, z=" + std::to_string(ta_foil_z_mm) + " mm") : std::string("not applicable (RUN_NUMBER unset or not a proton-scattering run)")) << std::endl;
+  std::cout << "OUT_DIR=" << outdir << " -> ta_foil_run=" << ta_foil_run
+             << " (z=" << ta_foil_z_mm << " mm), source_run=" << source_run << std::endl;
 
   // if (getenv("PC_ENERGY_CALIBRATION"))
   //   doPCEnergyCalibration = std::atoi(getenv("PC_ENERGY_CALIBRATION")) != 0;
@@ -530,6 +548,25 @@ void TrackRecon::Begin(TTree * /*tree*/)
   if (getenv("source_vertex"))
     source_vertex = (double)std::atof(std::string(getenv("source_vertex")).c_str());
 
+  if (doPCEnergyCalibration)
+  {
+    gSystem->mkdir("pc_calib_raw", kTRUE); // kTRUE = create parents, no-op if it exists
+    // Tag source-run vs proton-run files distinctly (src_/ap_ prefix) so the
+    // aggregator can treat them differently: source runs sit at a handful of
+    // known, fixed positions -- each run is a single clean, additive
+    // calibration point per wire, not thousands of samples of the same thing --
+    // while proton runs sample a continuously varying vertex/angle within a
+    // single run, where the per-event spread IS the useful signal and must
+    // stay pooled at the individual-point level.
+    std::string runTypeTag = source_run ? "src_" : (ta_foil_run ? "ap_" : "other_");
+    std::string tag = runTypeTag + (getenv("RUN_NUMBER") ? std::string("run") + getenv("RUN_NUMBER")
+                                                          : dataset + "_pid" + std::to_string(getpid()));
+    std::string outname = "pc_calib_raw/points_" + tag + ".dat";
+    pcCalibOutFile.open(outname);
+    pcCalibOutFile << std::scientific << std::setprecision(6);
+    std::cout << "PC energy calibration: streaming raw points to " << outname << std::endl;
+  }
+
   if (getenv("CO2percent"))
     co2pc = std::atoi(getenv("CO2percent"));
   std::cout << "CO2 percent set to  " << co2pc << std::endl;
@@ -787,7 +824,7 @@ inline double evalElossForward(TSpline3 *fwd, TSpline3 *inv, double E, double pa
   return e;
 }
 
-inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events)
+inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &QQQ_Events)
 {
   const TVector3 source_pos(beam_axis_x, beam_axis_y, source_vertex);
   for (const auto &pcevent : PC_Events)
@@ -795,26 +832,66 @@ inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events)
     if (!(pcevent.multi1 >= 1 && pcevent.multi2 >= 1))
       continue;
 
-    TVector3 interaction = pcevent.pos;
+    // Use an Si-hit-derived z (pczguess) instead of the PC's own cfrac
-    if (pcevent.multi1 == 1 && pcevent.multi2 == 1)
+    // sub-cell reconstruction: after depositing dE_gas in the PC, the alpha
+    // continues on and strikes SX3 or QQQ, giving a genuinely INDEPENDENT
+    // trajectory-angle measurement from the known, fixed source position --
+    // same pczguess pattern used elsewhere in this file for A1C0 benchmarks
+    // (each detector keeps its own established formula/geometry constants
+    // for consistency with the rest of the file). This sidesteps the
+    // PC-only reconstruction entirely (no more circular dependence on the
+    // PC's own charge-division z, and no more dependence on anode/cathode
+    // multiplicity for z precision), so the push gates below only need to
+    // protect ADC purity, not z precision. Try SX3 first, then QQQ.
+    bool foundSi = false;
+    double pcz = 0.0;
+    for (const auto &sx3event : SX3_Events)
     {
-      bool inband = false;
+      if (!(std::isfinite(sx3event.Time1) && std::isfinite(pcevent.Time1)) || TMath::Abs(sx3event.Time1 - pcevent.Time1) > 150.0)
-      double pcz = a1c1_cfrac_pcz(pcevent, source_pos, inband);
-      if (!inband)
         continue;
-      interaction.SetZ(pcz);
+      if (TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) > TMath::Pi() / 3.0)
+        continue;
+      double sx3theta = TMath::ATan2(88.0, sx3event.pos.Z() - source_vertex);
+      pcz = 37.0 / TMath::Tan(sx3theta) + source_vertex;
+      foundSi = true;
+      break;
     }
+    if (!foundSi)
+    {
+      for (const auto &qqqevent : QQQ_Events)
+      {
+        if (!(std::isfinite(qqqevent.Time1) && std::isfinite(pcevent.Time1)) || TMath::Abs(qqqevent.Time1 - pcevent.Time1) > 150.0)
+          continue;
+        if (TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) > TMath::Pi() / 4.0)
+          continue;
+        double qqqTheta = (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta();
+        pcz = 37.0 / TMath::Tan(qqqTheta) + source_vertex;
+        foundSi = true;
+        break;
+      }
+    }
+    if (!foundSi)
+      continue;
 
+    TVector3 interaction(pcevent.pos.X(), pcevent.pos.Y(), pcz);
     double path_length = pathLengthCm(source_pos, interaction);
     double e_remaining = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, pc_calib_alpha_source_mev, path_length);
     double dE_gas = pc_calib_alpha_source_mev - e_remaining;
 
     if (!std::isfinite(dE_gas) || dE_gas <= 0.0)
       continue;
-    if (pcevent.Anodech >= 0 && pcevent.Anodech < 24)
+    // Energy1/Energy2 are per-wire ADC. When multi1/multi2 > 1 the charge is
-      pcCalibData[pcevent.Anodech].push_back({pcevent.Energy1, dE_gas});
+    // shared across several wires, so a single wire's Energy here is only
-    if (pcevent.Cathodech >= 0 && pcevent.Cathodech < 24)
+    // its (event-dependent) share of the total -- not comparable to the
-      pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, dE_gas});
+    // full predicted dE_gas. Restrict each side to clean, single-wire
+    // clusters so ADC and dE_gas both refer to the whole deposited charge.
+    // z now comes from SX3 independently of multi1/multi2, so anode and
+    // cathode purity are gated independently (mirrors the proton
+    // accumulator's tryEvent).
+    if (pcevent.multi1 == 1 && pcevent.Anodech >= 0 && pcevent.Anodech < 24)
+      pcCalibWritePoint(pcevent.Anodech, pcevent.Energy1, dE_gas);
+    if (pcevent.multi2 >= 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24)
+      pcCalibWritePoint(24 + pcevent.Cathodech, pcevent.Energy2, dE_gas);
   }
 }
 
@@ -848,7 +925,7 @@ inline double predictElasticEnergy(Kinematics &kin, double angle3_deg, double t3
 // reading the PC's own (uncalibrated) energy, predicts the alpha's energy
 // from elastic kinematics at the reconstructed angle/vertex and walks that
 // forward through the gas to a predicted dE_gas -- same target quantity as
-// the source calibration, merged into the same pcCalibData[wire] clouds.
+// the source calibration, streamed into the same per-wire output file.
 inline void pcEnergyCalibrationAccumulateProton(const std::vector<Event> &PC_Events, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events)
 {
   if (!ta_foil_run)
@@ -902,9 +979,9 @@ inline void pcEnergyCalibrationAccumulateProton(const std::vector<Event> &PC_Eve
     // Single-wire only, so each wire's cloud maps its own charge (see the
     // source-run accumulator for the same reasoning).
     if (pcevent.multi1 == 1 && pcevent.Anodech >= 0 && pcevent.Anodech < 24)
-      pcCalibData[pcevent.Anodech].push_back({pcevent.Energy1, dE_gas});
+      pcCalibWritePoint(pcevent.Anodech, pcevent.Energy1, dE_gas);
-    if (pcevent.Cathodech >= 0 && pcevent.Cathodech < 24)
+    if (pcevent.multi2 == 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24)
-      pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, dE_gas});
+      pcCalibWritePoint(24 + pcevent.Cathodech, pcevent.Energy2, dE_gas);
   };
 
   for (const auto &pcevent : PC_Events)
@@ -1396,6 +1473,38 @@ Bool_t TrackRecon::Process(Long64_t entry)
 
   for (const auto &aCluster : aClusters)
   {
+    // A2 anode charge-sharing ratio: for exactly-two-fired-wire (multi-anode)
+    // clusters, the ratio of the smaller to the larger wire's energy. This
+    // characterises how the deposited charge splits between the two straddled
+    // anode wires (~0 = one wire dominates, ~1 = even split) and is filled once
+    // per anode cluster (before the cathode loop) so it isn't inflated by the
+    // cathode multiplicity.
+    if (aCluster.size() == 2)
+    {
+      double ae0 = std::get<1>(aCluster[0]);
+      double ae1 = std::get<1>(aCluster[1]);
+      double alo = std::min(ae0, ae1);
+      double ahi = std::max(ae0, ae1);
+      if (ahi > 0.0)
+      {
+        double aratio = alo / ahi;
+        plotter->Fill1D("A2_anode_ratio", 120, 0, 1.2, aratio, "hGMPC");
+        plotter->Fill2D("A2_anode_ratio_vs_sum", 800, 0, 40000, 120, 0, 1.2, ae0 + ae1, aratio, "hGMPC");
+        plotter->Fill2D("A2_anode_ratio_vs_lowerIndex", 24, 0, 24, 120, 0, 1.2,
+                        std::min(std::get<0>(aCluster[0]), std::get<0>(aCluster[1])), aratio, "hGMPC");
+      }
+      // Raw-vs-raw sanity check for the "A2 total = 2x A1 total" question, BEFORE
+      // any energy-calibration math touches it: if the raw (gain-matched-only)
+      // 2-wire sum already peaks at ~2x the raw single-wire distribution, the
+      // effect is upstream of the calibration code (neighbor-wire crosstalk on
+      // the raw signal, or a genuine double-hit) rather than introduced by how
+      // the calibration sums wires.
+      plotter->Fill1D("Raw_A2_AnodeSum", 800, 0, 40000, ae0 + ae1, "hGMPC");
+    }
+    else if (aCluster.size() == 1)
+    {
+      plotter->Fill1D("Raw_A1_AnodeSum", 800, 0, 40000, std::get<1>(aCluster[0]), "hGMPC");
+    }
     for (const auto &cCluster : cClusters)
     {
       if (aCluster.size() == 0)
@@ -1419,69 +1528,51 @@ Bool_t TrackRecon::Process(Long64_t entry)
         if (pcEnergyCalibLoaded)
         {
           Event PCEventCalibrated = PCEvent;
-
+          // Calibrate EACH anode wire's charge with its own slope and sum the
+          // results, rather than applying the first wire's factor to the raw
+          // cluster sum: the alpha's total gas dE is the charge summed over the
+          // wires it straddles, and whether it lands on one wire or two is
+          // phi-correlated, so a single-factor-on-sum made the calibrated anode
+          // energy depend on phi (unphysical: dE depends only on theta,z).
+          //
+          // The intercept, though, is a per-EVENT baseline (pedestal/threshold
+          // offset from the single-wire fit, where one wire's ADC stood for the
+          // whole dE_gas) -- not a per-wire quantity. Adding intercept[wi] once
+          // per fired wire double-counts that baseline for multi-wire clusters,
+          // reintroducing a multiplicity- (hence indirectly phi-) dependent bias
+          // through the back door. Apply it exactly once, from the primary
+          // (max-energy, i.e. aCluster[0]) wire.
           double anodeCalibSum = 0.0;
-          double primaryIntercept = 0.0;
+          std::vector<double> calibWireEnergies; // per-wire slope*ADC (no intercept), for the A2 ratio below
-
-          // Vector to track individual calibrated fractional energies for A2 charge sharing
-          std::vector<double> calibWireEnergies;
           calibWireEnergies.reserve(aCluster.size());
-
-          if (!aCluster.empty())
-          {
-            // Use the intercept of the primary (maximum energy) wire exactly ONCE.
-            // The first element of aCluster from Make_Clusters is the maximum energy wire.
-            int primaryWire = std::get<0>(aCluster[0]);
-            if (primaryWire >= 0 && primaryWire < 24)
-            {
-              primaryIntercept = pcEnergyIntercept[primaryWire];
-            }
-          }
-
-          // Apply per-wire gain slopes to the fractional charge, but DO NOT add the intercept yet
           for (const auto &w : aCluster)
           {
             int wi = std::get<0>(w);
-            if (wi >= 0 && wi < 24)
+            double wCalibE = (wi >= 0 && wi < 24) ? pcEnergySlope[wi] * std::get<1>(w) : 0.0;
-            {
+            anodeCalibSum += wCalibE;
-              // Note: std::get<1>(w) already has the relative matching 'pcSlope' applied
+            calibWireEnergies.push_back(wCalibE);
-              double wCalibE = pcEnergySlope[wi] * std::get<1>(w);
-              anodeCalibSum += wCalibE;
-              calibWireEnergies.push_back(wCalibE);
-            }
           }
-
+          int primaryAnodeWire = std::get<0>(aCluster[0]);
-          // Apply the baseline offset (intercept) exactly once for the whole event
+          double primaryIntercept = (primaryAnodeWire >= 0 && primaryAnodeWire < 24) ? pcEnergyIntercept[primaryAnodeWire] : 0.0;
           PCEventCalibrated.Energy1 = anodeCalibSum + primaryIntercept;
-
+          // Cathode uses the single max wire (cpMaxE) -- indexed by z, so it's
-          // Cathode uses the single max wire (cpMaxE) -- indexed by z, so it's already phi-consistent
+          // already phi-consistent; leave it as-is.
           PCEventCalibrated.Energy2 = pcEnergySlope[24 + PCEvent.Cathodech] * cpMaxE + pcEnergyIntercept[24 + PCEvent.Cathodech];
-
           PC_Events_calibrated.push_back(PCEventCalibrated);
 
-          // ---------------------------------------------------------
+          // Calibrated-energy A2 charge-sharing ratio: same diagnostic as the raw
-          // A2 Anode Charge-Sharing Diagnostics
+          // version above (A2_anode_ratio), but on the per-wire CALIBRATED shares,
-          // ---------------------------------------------------------
+          // so miscalibration between the two wires shows up as a ratio pulled
-          if (aCluster.size() == 2 && calibWireEnergies.size() == 2)
+          // away from what the raw-ADC ratio would give. Checks phi/energy
+          // dependence directly on the quantity that actually feeds Energy1.
+          if (calibWireEnergies.size() == 2)
           {
-            double e1 = calibWireEnergies[0];
+            double eSmaller = std::min(calibWireEnergies[0], calibWireEnergies[1]);
-            double e2 = calibWireEnergies[1];
+            double eLarger = std::max(calibWireEnergies[0], calibWireEnergies[1]);
-
-            double eSmaller = std::min(e1, e2);
-            double eLarger = std::max(e1, e2);
-
-            // Ratio goes from 0.0 (almost all charge on one wire) to 1.0 (perfectly shared)
             double ratio = (eLarger > 0.0) ? (eSmaller / eLarger) : 0.0;
-
-            // Plotting the 1D ratio distribution
             plotter->Fill1D("Calib_A2_AnodeRatio", 200, 0.0, 1.0, ratio, "hCalibPC");
-
-            // Checking phi-dependence: Should peak at ratio=1.0 strictly between wires
-            // and approach 0.0 when passing directly over a wire
             plotter->Fill2D("Calib_A2_AnodeRatio_vs_Phi", 360, -180, 180, 200, 0.0, 1.0,
                             PCEvent.pos.Phi() * 180.0 / M_PI, ratio, "hCalibPC");
-
-            // Energy dependence: Verifies if charge induction geometry changes via delta-rays/track-length
             plotter->Fill2D("Calib_A2_AnodeRatio_vs_TotalE", 400, 0, 10, 200, 0.0, 1.0,
                             PCEventCalibrated.Energy1, ratio, "hCalibPC");
           }
@@ -1541,14 +1632,18 @@ Bool_t TrackRecon::Process(Long64_t entry)
 
       if (pcEnergyCalibLoaded)
       {
-        // Per-wire-then-sum, phi-independent (same reasoning as the crossover branch).
+        // Per-wire-then-sum, intercept applied once from the primary wire --
+        // same reasoning as the crossover branch above.
         double anodeCalibSum = 0.0;
         for (const auto &w : aCl)
         {
           int wi = std::get<0>(w);
           if (wi >= 0 && wi < 24)
-            anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w) + pcEnergyIntercept[wi];
+            anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w);
         }
+        int primaryAnodeWireA1C0 = std::get<0>(aCl[0]);
+        double primaryInterceptA1C0 = (primaryAnodeWireA1C0 >= 0 && primaryAnodeWireA1C0 < 24) ? pcEnergyIntercept[primaryAnodeWireA1C0] : 0.0;
+        anodeCalibSum += primaryInterceptA1C0;
         Event ev(pc, anodeCalibSum, -1.0, apTSMaxE, -1.0);
         ev.multi1 = static_cast<int>(aCl.size());
         ev.multi2 = 0; // no cathode -> a{n}c0 topology in pcCalibratedHistograms
@@ -1560,20 +1655,28 @@ Bool_t TrackRecon::Process(Long64_t entry)
       // Anode-wire calibration point -- source runs only. The fixed alpha-source
       // energy is only valid there; proton-run A1C0 has no elastic tag to predict
       // its energy, so it contributes to the display but not the fit.
-      if (doPCEnergyCalibration && !ta_foil_run)
+      if (doPCEnergyCalibration && source_run)
       {
         double path = pathLengthCm(source_pos_a1c0, pc);
         double e_rem = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, pc_calib_alpha_source_mev, path);
         double dE_gas = pc_calib_alpha_source_mev - e_rem;
         if (std::isfinite(dE_gas) && dE_gas > 0.0)
-          pcCalibData[anodeIdx].push_back({apSumE, dE_gas});
+          pcCalibWritePoint(anodeIdx, apSumE, dE_gas);
       }
     }
   }
 
   if (doPCEnergyCalibration)
   {
-    pcEnergyCalibrationAccumulate(PC_Events);
+    // pcEnergyCalibrationAccumulate assumes source_pos = the FIXED alpha-source
+    // position -- only true for the source runs. On a proton-scattering run,
+    // source_vertex is the beam-entrance placeholder, not a real source, so
+    // running this there silently wrote bogus points (huge path length, usually
+    // clamped to the full 5.486 MeV) into the same pool as the real
+    // kinematics-derived proton points. Gate it to non-proton-campaign runs;
+    // pcEnergyCalibrationAccumulateProton already self-gates the other way.
+    if (source_run)
+      pcEnergyCalibrationAccumulate(PC_Events, SX3_Events, QQQ_Events);
     pcEnergyCalibrationAccumulateProton(PC_Events, QQQ_Events, SX3_Events);
   }
 
@@ -1745,26 +1848,14 @@ void TrackRecon::Terminate()
 {
   plotter->FlushToDisk(10);
 
-  if (doPCEnergyCalibration)
+  if (doPCEnergyCalibration && pcCalibOutFile.is_open())
   {
-    gSystem->mkdir("pc_calib_raw", kTRUE); // kTRUE = create parents, no-op if it exists
+    // Points were streamed directly to disk as they were generated (see
-    std::string tag = getenv("RUN_NUMBER") ? std::string("run") + getenv("RUN_NUMBER")
+    // pcCalibWritePoint / Begin()) instead of buffered in memory for the
-                                           : dataset + "_pid" + std::to_string(getpid());
+    // whole run, so there's nothing left to write here -- just close it.
-    std::string outname = "pc_calib_raw/points_" + tag + ".dat";
+    pcCalibOutFile.close();
-    std::ofstream outfile(outname);
+    std::cout << "PC energy calibration: closed raw points file -- run "
-    outfile << std::scientific << std::setprecision(6);
+              << "pccal/fit_pc_energy_calibration.C once all calibration runs are done"
-    long long nPoints = 0;
-    for (int wire = 0; wire < 48; ++wire)
-    {
-      for (const auto &p : pcCalibData[wire])
-      {
-        outfile << wire << " " << p.first << " " << p.second << "\n";
-        ++nPoints;
-      }
-    }
-    outfile.close();
-    std::cout << "PC energy calibration: appended " << nPoints << " raw points to " << outname
-              << " -- run pccal/fit_pc_energy_calibration.C once all calibration runs are done"
               << " to (re)produce pc_energy_calibration.dat" << std::endl;
   }
 }
@@ -1890,7 +1981,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
     const bool hasCathode = (pcevent.Cathodech >= 0);
     const double totalE = hasCathode ? (pcevent.Energy1 + pcevent.Energy2) : pcevent.Energy1;
     if (hasCathode)
-      plotter->Fill2D("Calib_AnodeE_vs_CathodeE_a1c1andup", 800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC");
+      plotter->Fill2D("Calib_AnodeE_vs_CathodeE_a1c1andup",800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC");
     for (const std::string &t : {std::string(""), topo})
     {
       plotter->Fill2D("Calib_AnodeE_vs_AnodeIndex" + t, 24, 0, 24, 800, 0, 3, pcevent.Anodech, pcevent.Energy1, "hCalibPC");
@@ -1901,7 +1992,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
       {
         plotter->Fill2D("Calib_CathodeE_vs_CathodeIndex" + t, 24, 0, 24, 800, 0, 3, pcevent.Cathodech, pcevent.Energy2, "hCalibPC");
         plotter->Fill1D("Calib_CathodeE" + t, 800, 0, 3, pcevent.Energy2, "hCalibPC");
-        plotter->Fill2D("Calib_AnodeE_vs_CathodeE" + t, 800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC");
+        plotter->Fill2D("Calib_AnodeE_vs_CathodeE" + t,800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC");
       }
 
       for (const auto &qqqevent : QQQ_Events)

pc_energy_calibration.dat

@@ -0,0 +1,48 @@
+0 8.024604e-05 0.000000e+00
+1 7.374085e-05 0.000000e+00
+2 9.583027e-05 0.000000e+00
+3 9.008263e-05 0.000000e+00
+4 8.452456e-05 0.000000e+00
+5 1.835247e-04 0.000000e+00
+6 5.125204e-05 0.000000e+00
+7 1.081530e-04 0.000000e+00
+8 1.269708e-04 0.000000e+00
+9 2.018195e-04 0.000000e+00
+10 1.089829e-04 0.000000e+00
+11 8.956875e-05 0.000000e+00
+12 1.000000e+00 0.000000e+00
+13 7.299374e-05 0.000000e+00
+14 7.270788e-05 0.000000e+00
+15 8.406974e-05 0.000000e+00
+16 8.997978e-05 0.000000e+00
+17 4.542671e-05 0.000000e+00
+18 4.755692e-05 0.000000e+00
+19 2.138104e-04 0.000000e+00
+20 9.264217e-05 0.000000e+00
+21 1.701897e-04 0.000000e+00
+22 2.162763e-04 0.000000e+00
+23 3.219591e-04 0.000000e+00
+24 7.364520e-05 0.000000e+00
+25 9.343009e-05 0.000000e+00
+26 1.013807e-04 0.000000e+00
+27 1.495574e-04 0.000000e+00
+28 1.166817e-04 0.000000e+00
+29 9.196472e-05 0.000000e+00
+30 1.030198e-04 0.000000e+00
+31 1.007519e-04 0.000000e+00
+32 1.026175e-04 0.000000e+00
+33 8.109499e-05 0.000000e+00
+34 8.137172e-05 0.000000e+00
+35 8.502646e-05 0.000000e+00
+36 1.042813e-04 0.000000e+00
+37 3.213148e-04 0.000000e+00
+38 9.069578e-05 0.000000e+00
+39 1.600775e-04 0.000000e+00
+40 1.357355e-04 0.000000e+00
+41 1.812490e-04 0.000000e+00
+42 1.186314e-04 0.000000e+00
+43 1.196526e-04 0.000000e+00
+44 2.486180e-04 0.000000e+00
+45 8.347962e-05 0.000000e+00
+46 7.408129e-05 0.000000e+00
+47 8.388388e-05 0.000000e+00

pccal/fit_pc_energy_calibration.C

@@ -1,29 +1,26 @@
-// Aggregates every run's raw PC energy-calibration points (written by
-// TrackRecon.C's Terminate() into pc_calib_raw/points_*.dat, one file per run
-// so parallel jobs never collide) into a single per-wire linear fit, pooling
-// BOTH the 17F and 27Al datasets together -- the PC's wire-level gain isn't
-// expected to differ meaningfully between the two campaigns, so calibration
-// statistics from either dataset's alpha-source or proton-scattering runs are
-// combined into one fit rather than kept separate.
-//
-// Run non-interactively from the top-level ANASEN-Analysis directory once all
-// desired calibration runs have completed:
-//   root -q -l -b -e '.L pccal/fit_pc_energy_calibration.C' -e 'fit_pc_energy_calibration()'
-//
-// Writes pc_energy_calibration.dat (wire slope intercept, 48 rows), which
-// TrackRecon.C's Begin() loads automatically on the next run to populate
-// PC_Events_calibrated. Re-run this any time pc_calib_raw/ has new points --
-// it always re-reads and re-fits everything currently on disk, so it's safe
-// to call repeatedly as more calibration runs accumulate.
 #include <TSystemDirectory.h>
 #include <TSystemFile.h>
 #include <TList.h>
 #include <TString.h>
+#include <TCanvas.h>
+#include <TH2D.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <TFile.h>
+#include <TStyle.h>
+#include <TSystem.h>
+#include <TROOT.h>
+#include <TPaveText.h>
 #include <fstream>
 #include <iostream>
 #include <vector>
+#include <iomanip>
+#include <cmath>
 
-void fit_pc_energy_calibration()
+// Optional dataset_filter (e.g., "17F" or "27Al") prevents mixing different
+// gas pressure/temperature environments which causes gain smearing.
+// Leave blank ("") to pool all files.
+void fit_pc_energy_calibration(const std::string& dataset_filter = "")
 {
   std::vector<std::pair<double, double>> pts[48]; // [wire] -> (ADC, dE_gas MeV)
 
@@ -37,55 +34,180 @@ void fit_pc_energy_calibration()
   }
 
   int nFiles = 0;
+  long long nOverflowCut = 0;
   TIter next(files);
   TSystemFile *f;
+  
+  // Read all points from all files into the single pooled array
   while ((f = (TSystemFile *)next()))
   {
     TString name = f->GetName();
     if (f->IsDirectory() || !name.BeginsWith("points_") || !name.EndsWith(".dat"))
       continue;
+      
+    // Apply dataset safeguard if requested
+    if (!dataset_filter.empty() && !name.Contains(dataset_filter.c_str()))
+      continue;
+
     std::ifstream infile(std::string("pc_calib_raw/") + name.Data());
     if (!infile.is_open())
       continue;
+      
     int wire;
     double adc, dE_gas;
     while (infile >> wire >> adc >> dE_gas)
+    {
       if (wire >= 0 && wire < 48)
-        pts[wire].push_back({adc, dE_gas});
+      {
+        if (adc < 64000.0) {
+          pts[wire].push_back({adc, dE_gas});
+        } else {
+          nOverflowCut++;
+        }
+      }
+    }
     ++nFiles;
   }
-  std::cout << "fit_pc_energy_calibration: read " << nFiles << " run file(s) from pc_calib_raw/" << std::endl;
+  
+  std::cout << "fit_pc_energy_calibration: read " << nFiles 
+            << " run file(s) from pc_calib_raw/ (Filter: '" << dataset_filter << "')" << std::endl;
+  std::cout << "fit_pc_energy_calibration: cut " << nOverflowCut 
+            << " points due to ADC >= 64k overflow." << std::endl;
+
+  // --- Setup ROOT Diagnostic Graphics ---
+  gROOT->SetBatch(kTRUE); // Run silently without popping up windows
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(kBird); 
+  
+  // Create output directory for individual PNGs
+  gSystem->mkdir("pc_calib_plots", kTRUE);
+  
+  TFile *fOut = new TFile("pc_calib_diagnostics.root", "RECREATE");
+  
+  TCanvas *cAnodes = new TCanvas("cAnodes", "Anode Calibrations", 1800, 1200);
+  cAnodes->Divide(6, 4, 0.01, 0.01);
+  
+  TCanvas *cCathodes = new TCanvas("cCathodes", "Cathode Calibrations", 1800, 1200);
+  cCathodes->Divide(6, 4, 0.01, 0.01);
 
   std::ofstream outfile("pc_energy_calibration.dat");
   outfile << std::scientific << std::setprecision(6);
+  
+  // --- Fit and Plot each wire ---
   for (int wire = 0; wire < 48; ++wire)
   {
     double slope = 1.0, intercept = 0.0;
-    bool ok = pts[wire].size() >= 2;
+    double n = static_cast<double>(pts[wire].size());
+    bool ok = (n >= 2);
+    
+    // 1. Calculate Standard Linear Least Squares Fit (y = mx) where intercept is forced to 0
+    double maxX = 0.0, maxY = 0.0;
+    double sxx = 0, sxy = 0; // Only need sum(x^2) and sum(x*y) for fixed-0 intercept
+
+    for (const auto &p : pts[wire])
+    {
+      if (p.first > maxX) maxX = p.first;
+      if (p.second > maxY) maxY = p.second;
+      sxx += p.first * p.first;
+      sxy += p.first * p.second;
+    }
+
     if (ok)
     {
-      double sx = 0, sy = 0, sxx = 0, sxy = 0;
+      if (std::isfinite(sxx) && std::abs(sxx) > 1e-12)
-      for (const auto &p : pts[wire])
       {
-        sx += p.first;
+        slope = sxy / sxx;
-        sy += p.second;
+        intercept = 0.0; // Forced mathematically
-        sxx += p.first * p.first;
-        sxy += p.first * p.second;
       }
-      double n = static_cast<double>(pts[wire].size());
+      else
-      double denom = n * sxx - sx * sx;
-      ok = std::isfinite(denom) && std::abs(denom) > 1e-12;
-      if (ok)
       {
-        slope = (n * sxy - sx * sy) / denom;
+        ok = false;
-        intercept = (sy - slope * sx) / n;
       }
     }
-    if (!ok)
+
-      std::cerr << "fit_pc_energy_calibration: wire " << wire << " has too few points (" << pts[wire].size()
+    if (!ok) {
+      std::cerr << "fit_pc_energy_calibration: wire " << wire << " has too few points (" << n
                 << ") to fit -- writing identity (slope=1, intercept=0)" << std::endl;
+    }
+
     outfile << wire << " " << slope << " " << intercept << "\n";
+
+    // 2. Generate and Fill 2D Density Histogram
+    if (maxX <= 0) maxX = 64000.0;
+    if (maxY <= 0) maxY = 10.0;
+
+    TString wName = Form("%s %02d", wire < 24 ? "Anode" : "Cathode", wire < 24 ? wire : wire - 24);
+    TH2D *h2 = new TH2D(Form("h2_wire_%d", wire), 
+                        wName + "; ADC; dE_{gas} (MeV)", 
+                        150, 0, maxX * 1.05, 
+                        150, 0, maxY * 1.05);
+
+    for (const auto &p : pts[wire]) {
+      h2->Fill(p.first, p.second);
+    }
+
+    // 3. Setup Fit Line and Stats Box
+    TF1 *fitLine = nullptr;
+    TPaveText *pt = nullptr;
+    
+    if (n > 0) {
+        fitLine = new TF1(Form("fit_%d", wire), "[0]*x", 0, maxX * 1.05); // Formula is strictly y = m*x
+        fitLine->SetParameter(0, slope);
+        fitLine->SetLineColor(kRed);
+        fitLine->SetLineWidth(2);
+
+        pt = new TPaveText(0.15, 0.75, 0.55, 0.88, "NDC");
+        pt->SetFillColor(kWhite);
+        pt->SetBorderSize(1);
+        pt->AddText(Form("N = %.0f", n));
+        pt->AddText(Form("m = %.2e", slope));
+        pt->AddText("b = 0 (Fixed)");
+    }
+
+    // ------------------------------------------------------------
+    // 4. Save High-Res Individual PNG
+    // ------------------------------------------------------------
+    TCanvas cTemp("cTemp", "cTemp", 800, 600);
+    cTemp.SetGridx();
+    cTemp.SetGridy();
+    if (n > 0) {
+        h2->Draw("COLZ");
+        fitLine->Draw("SAME");
+        pt->Draw();
+    } else {
+        h2->SetTitle(wName + " (DEAD/NO DATA)");
+        h2->Draw();
+    }
+    cTemp.SaveAs(Form("pc_calib_plots/wire_%02d.png", wire));
+
+    // ------------------------------------------------------------
+    // 5. Draw onto Global PDF Canvas
+    // ------------------------------------------------------------
+    TVirtualPad* pad = (wire < 24) ? cAnodes->cd(wire + 1) : cCathodes->cd(wire - 24 + 1);
+    pad->SetGridx();
+    pad->SetGridy();
+
+    if (n > 0) 
+    {
+      h2->DrawClone("COLZ");
+      fitLine->DrawClone("SAME");
+      pt->DrawClone();
+    }
+    else 
+    {
+      h2->DrawClone(); 
+    }
+    
+    // Write histogram to ROOT file
+    fOut->cd();
+    h2->Write();
   }
+  
   outfile.close();
+  
+  fOut->Close();
+
   std::cout << "fit_pc_energy_calibration: wrote pc_energy_calibration.dat" << std::endl;
-}
+  std::cout << "fit_pc_energy_calibration: individual high-res PNGs saved to pc_calib_plots/" << std::endl;
+  std::cout << "fit_pc_energy_calibration: global contact sheets saved to pc_calib_anodes.pdf and pc_calib_cathodes.pdf" << std::endl;
+}

run_27Al.sh

@@ -6,12 +6,13 @@ export OUT_DIR="Output_27Al"
 export reactiondata=1
 export CO2percent=3
 export pressure_in_torr=250
-export CATHODE_GAIN=1
+export CATHODE_GAIN=3.0
+export source_vertex=-200.0
 export BEAM_AXIS_X=-15
 export BEAM_AXIS_Y=-5
 
 # Clean up previous runs
-rm -f ${OUT_DIR}/results_run*.root output_27Al.root
+rm -f ${OUT_DIR}/*.root
 
 echo "Pre-compiling TrackRecon.C safely on a single core..."
 root -q -l -b -e '.L TrackRecon.C++O'
@@ -23,7 +24,6 @@ process_run() {
     local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root"
     local out="${outdir}/results_run${wrun}.root"
 
-    # Ensure the directory exists so ROOT doesn't fail silently
     mkdir -p "$outdir"
 
     root -q -l -b -x "$infile" \
@@ -39,19 +39,28 @@ process_run() {
 export -f process_run
 
 echo "Starting parallel processing..."
-parallel --bar -j 4 process_run ::: {50..52}
+time parallel --bar -j 8 process_run ::: {50..59}
+time parallel --bar -j 8 process_run ::: 62 63 66 67 73 74
+# time parallel --bar -j 1 run_once {1} ::: 68
+# time parallel --bar -j 6 run_once {1} ::: {78..89}
 
 echo "Merging files..."
 hadd -k -j 4 ${OUT_DIR}/output_27Al.root ${OUT_DIR}/results_run*.root
 
+rootbrowse ${OUT_DIR}/output_27Al.root
+
 unset DATASET
+unset PREFIX
+unset OUT_DIR
 unset reactiondata
 unset CO2percent
 unset pressure_in_torr
 unset CATHODE_GAIN
+unset source_vertex
 unset A1C1_LOWBAND_RFACTOR
 unset A1C1_Z_SCALE_QQQ
 unset A1C1_Z_OFF_QQQ
 unset A1C1_Z_OFF_SX3
 unset BEAM_AXIS_X
-unset BEAM_AXIS_Y
+unset BEAM_AXIS_Y
+echo "Script execution finished."

run_tr.sh

@@ -63,7 +63,7 @@ if [[ 1 -eq 1 ]]; then
     export pressure_in_torr=350
     rm -f ${OUT_DIR}/all.root
     echo "Processing 27Al alpha+gas runs..."
-    export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope"
+    # export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope"
     unset timecut_high
     export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope"
     unset Gain
@@ -89,7 +89,7 @@ if [[ 1 -eq 1 ]]; then
 fi
 
 # --- Block 5: 27Al Protons+Gas Runs (15, 17-22) ---
-if [[ 1 -eq 0 ]]; then
+if [[ 1 -eq 1 ]]; then
 
     # export CO2percent=4
     export DATASET="27Al"
@@ -101,15 +101,15 @@ if [[ 1 -eq 0 ]]; then
     echo "Starting parallel processing for 27Al proton runs..."
 
     # process_run 18
-    parallel --bar -j 8 process_run ::: 15 {17..22}
+    # parallel --bar -j 8 process_run ::: 15 {17..22}
-    # parallel --bar -j 8 process_run :::  {17..22}
+    parallel --bar -j 8 process_run :::  {17..22}
     hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{15..22}.root
     unset CATHODE_GAIN
     # exit
 fi
 
 # --- Block 6: 17F Proton Data  ---
-if [[ 1 -eq 0 ]]; then
+if [[ 1 -eq 1 ]]; then
     export DATASET="17F"
     export PREFIX="ProtonRun_"
     export OUT_DIR="Output_p"