StripSumScatter.cpp

Go to the documentation of this file.

#include "StripSumScatter.hpp"
#include "RegionCuts.hpp"
#include <TParameter.h>
 
StripSumScatter::StripSumScatter() {
  for (Int_t i = 0; i < 64; i++) {
    m_yLo[i] = 0.0;
    m_yHi[i] = 0.0;
  }
  m_nSeen = 0;
  m_nNormed = 0;
  m_normedAt.assign(
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX -
          Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN + 1,
      0);
  m_tagged.assign(
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX -
          Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN + 1,
      0);
}
 
StripSumScatter::~StripSumScatter() {
  std::map<Int_t, TH2F *>::iterator it;
  for (it = m_scatter.begin(); it != m_scatter.end(); ++it)
    delete it->second;
  m_scatter.clear();
}
 
Int_t StripSumScatter::ReacIndex(Int_t reac) {
  return reac - Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
}
 
Int_t StripSumScatter::YLoOf(Int_t reac) { return reac + 1; }
 
Int_t StripSumScatter::YHiOf(Int_t reac) {
  const StripSumScatterConfig &C = Constants::cfg.STRIP_SUM_SCATTER_CONFIG;
  std::map<Int_t, Int_t>::const_iterator it = C.POST_WINDOW_STRIPS.find(reac);
  if (it != C.POST_WINDOW_STRIPS.end())
    return TMath::Min(reac + TMath::Max(1, it->second), 17);
  // The window shrinks where it would otherwise run past the last strip it
  // may reach, and is never shorter than one strip.
  const Int_t hi = TMath::Min(reac + C.POST_TRIGGER_SUM_STRIPS,
                              TMath::Min(C.POST_WINDOW_LAST_STRIP, 17));
  return TMath::Max(hi, reac + 1);
}
 
Double_t StripSumScatter::s_jumpSigma[18] = {0.0};
Double_t StripSumScatter::s_stripSigma[18] = {0.0};
 
Double_t StripSumScatter::JumpSigma(Int_t strip) {
  return (strip >= 0 && strip < 18) ? s_jumpSigma[strip] : 0.0;
}
 
Double_t StripSumScatter::StripSigma(Int_t strip) {
  return (strip >= 0 && strip < 18) ? s_stripSigma[strip] : 0.0;
}
 
void StripSumScatter::SetJumpSigma(const Double_t *sigma) {
  for (Int_t s = 0; s < 18; s++)
    s_jumpSigma[s] = sigma[s];
}
 
void StripSumScatter::SetStripSigma(const Double_t *sigma) {
  for (Int_t s = 0; s < 18; s++)
    s_stripSigma[s] = sigma[s];
}
 
Double_t StripSumScatter::JumpMin(Int_t reac) {
  return Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REAC_JUMP_NSIGMA *
         JumpSigma(reac);
}
 
Bool_t StripSumScatter::BeamUpstreamOf(const EnergyView &ev, Int_t reac) {
  const StripSumScatterConfig &C = Constants::cfg.STRIP_SUM_SCATTER_CONFIG;
  if (!C.REQUIRE_BEAM_UPSTREAM_OF_REAC)
    return kTRUE;
  for (Int_t s = 1; s < reac; s++)
    if (TMath::Abs(ev.total[s] - 1.0) > C.BEAM_UPSTREAM_NSIGMA * StripSigma(s))
      return kFALSE;
  return kTRUE;
}
 
// Sigma-clipped width of a sample.
static Double_t ClippedWidth(const std::vector<Double_t> &values) {
  const Int_t kClipPasses = 3;
  const Double_t kClipNSigma = 3.0;
  Double_t mean = 0.0, width = 0.0;
  for (Int_t pass = 0; pass < kClipPasses; pass++) {
    Double_t sum = 0.0, sum2 = 0.0;
    Long64_t kept = 0;
    for (Int_t k = 0; k < Int_t(values.size()); k++) {
      const Double_t v = values[k];
      if (pass > 0 && TMath::Abs(v - mean) > kClipNSigma * width)
        continue;
      sum += v;
      sum2 += v * v;
      kept++;
    }
    if (kept < 2)
      break;
    mean = sum / Double_t(kept);
    width = TMath::Sqrt(TMath::Max(0.0, sum2 / Double_t(kept) - mean * mean));
  }
  return width;
}
 
Bool_t StripSumScatter::MeasureBeamNoise(TChain *chain, Double_t *jump_sigma,
                                         Double_t *strip_sigma) {
  const Long64_t kMaxEvents = 200000;
  const Long64_t kMinEvents = 1000;
  for (Int_t s = 0; s < 18; s++) {
    jump_sigma[s] = 0.0;
    strip_sigma[s] = 0.0;
  }
  if (!chain)
    return kFALSE;
 
  EnergyView ev;
  ev.Attach(chain);
  EnableEventBranches(chain);
  const Long64_t n = TMath::Min(chain->GetEntries(), kMaxEvents);
  std::vector<std::vector<Double_t>> diff(18), deposit(18);
  for (Long64_t j = 0; j < n; j++) {
    chain->GetEntry(j);
    ev.Decode();
    if (!AllStripsFired(ev) || IsPileup(ev) || IsNoise(ev))
      continue;
    if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REJECT_OFFBEAM && IsOffbeam(ev))
      continue;
    if (IsParityAsymmetric(ev))
      continue;
    deposit[0].push_back(ev.total[0]);
    for (Int_t s = 1; s < 18; s++) {
      diff[s].push_back(ev.total[s] - ev.total[s - 1]);
      deposit[s].push_back(ev.total[s]);
    }
  }
  chain->ResetBranchAddresses();
  if (Long64_t(deposit[0].size()) < kMinEvents)
    return kFALSE;
 
  // Reactions are a few 1e-4 of the sample, so the clipped widths are the
  // beam's noise.
  for (Int_t s = 0; s < 18; s++)
    strip_sigma[s] = ClippedWidth(deposit[s]);
  for (Int_t s = 1; s < 18; s++)
    jump_sigma[s] = ClippedWidth(diff[s]);
  for (Int_t s = 1; s <= 16; s++)
    if (!(jump_sigma[s] > 0.0) || !(strip_sigma[s] > 0.0))
      return kFALSE;
  return kTRUE;
}
 
void StripSumScatter::EnableEventBranches(TChain *chain) {
  chain->SetBranchStatus("*", 0);
  chain->SetBranchStatus("Left_0_17_dE", 1);
  chain->SetBranchStatus("RightdE", 1);
  chain->SetBranchStatus("Cathode", 1);
  // Absent in files built before the branch existed; harmless to enable.
  if (chain->GetBranch("SeedTs"))
    chain->SetBranchStatus("SeedTs", 1);
  // The parity-rejected grid diagnostic reads the Grid trigger channel. Enable
  // it only when that plot is requested so the common update path stays
  // branch-light (the grid is a Short_t; reading it otherwise is wasted I/O).
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PLOT_PARITY_REJECTED_GRID &&
      chain->GetBranch("Grid"))
    chain->SetBranchStatus("Grid", 1);
}
 
Bool_t StripSumScatter::AllStripsFired(const EnergyView &ev) {
  if (!Constants::cfg.IGNORE_STRIP_0 && !(ev.total[0] > 0.0))
    return kFALSE;
  if (!Constants::cfg.IGNORE_STRIP_17 && !(ev.total[17] > 0.0))
    return kFALSE;
  for (Int_t s = 1; s <= 16; s++)
    if (!(ev.total[s] > 0.0))
      return kFALSE;
  return kTRUE;
}
 
Bool_t StripSumScatter::PassesReaction(const EnergyView &ev, Int_t reac) {
  const Double_t kReacJumpMin = JumpMin(reac);
  const Double_t kReacJumpMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REAC_JUMP_MAX;
  const Double_t kSmoothMaxStep =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS_MAX_STEP;
  const Int_t kSmoothHiStrip =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS_END_STRIP;
  const Double_t kEndStripMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.END_STRIP_MAX;
 
  if (!AllStripsFired(ev))
    return kFALSE;
  // Otherwise a reaction at an earlier strip can pass this strip's jump gate
  // on a noise fluctuation and be counted here as well.
  if (!BeamUpstreamOf(ev, reac))
    return kFALSE;
  Double_t reac_jump = ev.total[reac] - ev.total[reac - 1];
  if (!(reac_jump > kReacJumpMin && reac_jump < kReacJumpMax))
    return kFALSE;
  if (!(ev.total[reac] > 1.0 + kReacJumpMin &&
        ev.total[reac] < 1.0 + kReacJumpMax))
    return kFALSE;
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS)
    for (Int_t s = reac + 1; s <= kSmoothHiStrip; s++)
      if (TMath::Abs(ev.total[s] - ev.total[s - 1]) > kSmoothMaxStep)
        return kFALSE;
  if (Constants::cfg.IGNORE_STRIP_17)
    return ev.total[16] < kEndStripMax;
  Int_t end_strip =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_STRIP_16_BELOW_BEAM ? 16
                                                                          : 17;
  return ev.total[end_strip] < kEndStripMax;
}
 
Bool_t StripSumScatter::IsPureBeam(const EnergyView &ev,
                                   const BeamEllipses &be) {
  if (!be.ok)
    return kFALSE;
  if (!AllStripsFired(ev))
    return kFALSE;
  // Must pass BOTH the entrance AND exit ellipses.
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PURE_BEAM_GATE ==
      StripSumScatterConfig::PURE_BEAM_GATE_S1_S2) {
    if (!PassesGate(be.s1_s2, ev, 1, 2))
      return kFALSE;
  } else {
    if (!PassesGate(be.s0_s1, ev, 0, 1))
      return kFALSE;
  }
  if (be.use_s15_s16) {
    if (!PassesGate(be.s15_s16, ev, 15, 16))
      return kFALSE;
  } else {
    if (!PassesGate(be.s16_s17, ev, 16, 17))
      return kFALSE;
  }
  return kTRUE;
}
 
Bool_t StripSumScatter::IsPileup(const EnergyView &ev) {
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REJECT_PILEUP) {
    const Double_t kThresh =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PILEUP_THRESH_PY;
    const Int_t kMinStrips =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PILEUP_MIN_STRIPS;
    Int_t n = 0;
    for (Int_t s = 1; s <= 16; s++)
      if (ev.total[s] >= kThresh && ++n >= kMinStrips)
        return kTRUE;
    return kFALSE;
  }
  const Double_t kThresh =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PILEUP_THRESHOLD;
  const Int_t kMinStrips =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PILEUP_MIN_STRIPS;
  Int_t n = 0;
  for (Int_t s = 1; s <= 16; s++)
    if (ev.total[s] >= kThresh && ++n >= kMinStrips)
      return kTRUE;
  return kFALSE;
}
 
Bool_t StripSumScatter::IsNoise(const EnergyView &ev) {
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REJECT_NOISE) {
    const Double_t kThresh =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.NOISE_THRESH_PY;
    const Int_t kMinStrips =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.NOISE_MIN_STRIPS;
    Int_t n = 0;
    for (Int_t s = 1; s <= 16; s++)
      if (ev.total[s] <= kThresh && ++n >= kMinStrips)
        return kTRUE;
    return kFALSE;
  }
  const Double_t kThresh =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.NOISE_THRESHOLD;
  const Int_t kMinStrips =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.NOISE_MIN_STRIPS;
  Int_t n = 0;
  for (Int_t s = 1; s <= 16; s++)
    if (ev.total[s] <= kThresh && ++n >= kMinStrips)
      return kTRUE;
  return kFALSE;
}
 
Bool_t StripSumScatter::IsOffbeam(const EnergyView &ev) {
  const Double_t kDist = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.OFFBEAM_DIST;
  const Int_t kMinStrips =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.OFFBEAM_MIN_STRIPS;
  // Use flat normed beam level (1.0) as reference, matching Python which
  // computes beam reference after these pre-filters.
  const Double_t kBeamRef = 1.0;
  Int_t n = 0;
  for (Int_t s = 1; s <= 16; s++)
    if (TMath::Abs(ev.total[s] - kBeamRef) >= kDist && ++n >= kMinStrips)
      return kTRUE;
  return kFALSE;
}
 
// Even strips against odd strips over the whole trace, before the tag; see
// PARITY_ASYM_MAX. Strip 16 is even and strip 1 odd, so both means run over
// eight strips; a residue's excess spreads over both parities alike and moves
// the ratio by a few percent at most.
Bool_t StripSumScatter::IsParityAsymmetric(const EnergyView &ev) {
  const Double_t kMax = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PARITY_ASYM_MAX;
  if (kMax <= 0.0)
    return kFALSE;
  Double_t odd = 0.0, even = 0.0;
  for (Int_t s = 1; s <= 16; s++)
    ((s % 2) ? odd : even) += ev.total[s];
  if (odd <= 0.0)
    return kTRUE;
  return TMath::Abs(even / odd - 1.0) > kMax;
}
 
Double_t StripSumScatter::SumRange(const Double_t *total, Int_t lo, Int_t hi) {
  Double_t sum = 0.0;
  for (Int_t s = lo; s <= hi; s++)
    sum += total[s];
  return sum;
}
 
void StripSumScatter::PlaneXY(const Double_t *total, Int_t reac, Double_t &x,
                              Double_t &y) {
  x = SumRange(total, Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO,
               Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI);
  y = SumRange(total, YLoOf(reac), YHiOf(reac));
}
 
std::vector<GateSpec> StripSumScatter::ActiveGates() {
  std::vector<GateSpec> gates;
  GateSpec g;
  g.sx = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_STRIP_X;
  g.sy = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_STRIP_Y;
  gates.push_back(g);
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_GATE_S3_S4) {
    g.sx = 3;
    g.sy = 4;
    gates.push_back(g);
  }
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_GATE_S5_S6) {
    g.sx = 5;
    g.sy = 6;
    gates.push_back(g);
  }
  return gates;
}
 
TString StripSumScatter::CacheName() {
  TString name = "StripSumScatter_cache";
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_GATE_S3_S4)
    name += "_g34";
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_GATE_S5_S6)
    name += "_g56";
  name += ".root";
  return name;
}
 
Bool_t StripSumScatter::PassesGate(const BeamFit2D &gate, const EnergyView &ev,
                                   Int_t sx, Int_t sy) {
  Double_t g0 = ev.total[sx];
  Double_t g1 = ev.total[sy];
  if (!(g0 > 0.0 && g1 > 0.0))
    return kFALSE;
  const Double_t kGateNSigmaX =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_X;
  const Double_t kGateNSigmaY =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_Y;
  return BeamFitUtils::InEllipseXY(gate, g0, g1, kGateNSigmaX, kGateNSigmaY);
}
 
BeamFit2D
StripSumScatter::FindBeamGate(TChain *chain, Int_t sx, Int_t sy,
                              const std::vector<GateSpec> &prior_specs,
                              const std::vector<BeamFit2D> &prior_gates,
                              const TString &tag, const TString &subdir) {
  const Int_t kGateBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_BINS;
  const Double_t kGateMin = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_MIN;
  const Double_t kGateMax = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_MAX;
  const Int_t kSeedHalfBins =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.SEED_HALF_BINS;
  const Double_t kSeedFrac = 0.3;
  const Long64_t kSampleMaxPoints =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.SAMPLE_MAX_POINTS;
 
  BeamFit2D out;
  EnergyView ev;
  ev.Attach(chain);
  EnableEventBranches(chain);
  TH2F *h =
      new TH2F(Form("h2_beamgate_s%d_s%d_%s", sx, sy, tag.Data()),
               Form(";#DeltaE strip %d [a.u.];#DeltaE strip %d [a.u.]", sx, sy),
               kGateBins, kGateMin, kGateMax, kGateBins, kGateMin, kGateMax);
  h->SetDirectory(nullptr);
  Long64_t n = chain->GetEntries();
  Long64_t stride = FileSet::SampleStride(n, kSampleMaxPoints);
  for (Long64_t j = 0; j < n; j += stride) {
    chain->GetEntry(j);
    ev.Decode();
    // Series gating: only events passing every prior gate feed this fit.
    Bool_t prior_ok = kTRUE;
    for (Int_t gi = 0; gi < Int_t(prior_specs.size()); gi++)
      if (!PassesGate(prior_gates[gi], ev, prior_specs[gi].sx,
                      prior_specs[gi].sy)) {
        prior_ok = kFALSE;
        break;
      }
    if (!prior_ok)
      continue;
    Double_t x = ev.total[sx];
    Double_t y = ev.total[sy];
    if (x > 0.0 && y > 0.0)
      h->Fill(x, y);
  }
  if (h->GetEntries() < 100) {
    delete h;
    return out;
  }
  Double_t bw_x = h->GetXaxis()->GetBinWidth(1);
  Double_t bw_y = h->GetYaxis()->GetBinWidth(1);
  Int_t bx = 0, by = 0, bz = 0;
  h->GetMaximumBin(bx, by, bz);
  Double_t peak_val = h->GetBinContent(bx, by);
  Int_t lo_bx = std::max(1, bx - kSeedHalfBins);
  Int_t hi_bx = std::min(h->GetNbinsX(), bx + kSeedHalfBins);
  Int_t lo_by = std::max(1, by - kSeedHalfBins);
  Int_t hi_by = std::min(h->GetNbinsY(), by + kSeedHalfBins);
  Moments2D m = BeamFitUtils::ComputeMoments(h, lo_bx, hi_bx, lo_by, hi_by,
                                             kSeedFrac * peak_val, bw_x, bw_y);
  if (m.weight <= 0) {
    delete h;
    return out;
  }
  out.amp = peak_val;
  out.mu_x = m.mu_x;
  out.mu_y = m.mu_y;
  out.sigma_x = m.sigma_x;
  out.sigma_y = m.sigma_y;
  out.rho = m.rho;
  out.ok = kTRUE;
 
  const Double_t kGateNSigmaX =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_X;
  const Double_t kGateNSigmaY =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_Y;
 
  std::lock_guard<std::mutex> lock(g_plot_mutex);
  TCanvas *c = PlottingUtils::GetConfiguredCanvas(kFALSE);
  PlottingUtils::ConfigureAndDraw2DHistogram(h, c);
  // Correlated 2D Gaussian ellipse matching the InEllipseXY gate.
  Double_t sxx = out.sigma_x * out.sigma_x;
  Double_t syy = out.sigma_y * out.sigma_y;
  Double_t sxy = out.rho * out.sigma_x * out.sigma_y;
  Double_t sum = sxx + syy;
  Double_t diff = sxx - syy;
  Double_t det = TMath::Sqrt(diff * diff + 4.0 * sxy * sxy);
  Double_t lambda1 = 0.5 * (sum + det);
  Double_t lambda2 = 0.5 * (sum - det);
  Double_t theta = 0.5 * TMath::ATan2(2.0 * sxy, diff) * 180.0 / TMath::Pi();
  TEllipse *e =
      new TEllipse(out.mu_x, out.mu_y, kGateNSigmaX * TMath::Sqrt(lambda1),
                   kGateNSigmaX * TMath::Sqrt(lambda2), 0, 360, theta);
  e->SetFillStyle(0);
  e->SetLineColor(kRed + 1);
  e->SetLineWidth(2);
  e->Draw();
  PlottingUtils::SaveFigure(c, Form("beam_gate_s%d_s%d", sx, sy), subdir,
                            PlotSaveOptions::kLINEAR);
  delete c;
 
  delete h;
  return out;
}
 
void StripSumScatter::DrawTraceSet(const std::vector<TGraph *> &traces,
                                   Int_t color) {
  for (Int_t i = 0; i < Int_t(traces.size()); i++) {
    traces[i]->SetLineColor(color);
    traces[i]->SetLineWidth(1);
    traces[i]->Draw("L SAME");
  }
}
 
// Re-render the SAME selected events under the other decode.
// IGNORE_SHORT_STRIPS is a decode-time switch, so both renderings come from one
// calibration: the long-only trace is long_au, and the summed trace is (long_au
// + short_au) rescaled per strip so that the summed BEAM peak sits at 1.0 a.u.
// Without that rescale the two are not on a common scale -- pass 1 anchors the
// long side's own beam peak at 1.0, so the sum reads 1+f with f the short-side
// fraction.
void StripSumScatter::DrawAltDecodeRegionTraces(Int_t reac, TCutG *cutAn,
                                                TCutG *cutAa) {
  const Bool_t long_only_is_current = Constants::cfg.IGNORE_SHORT_STRIPS;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kTracesPerRegion =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRACES_PER_CLASS;
  UInt_t bit = (1u << ReacIndex(reac));
 
  std::vector<const TraceEvt *> ev_beam, ev_aa, ev_an;
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    const TraceEvt &e = m_reservoir[k];
    if (e.beam_flat) {
      if (Int_t(ev_beam.size()) < kTracesPerRegion)
        ev_beam.push_back(&e);
      continue;
    }
    if (!(e.reac_mask & bit))
      continue;
    Double_t td[18];
    for (Int_t s = 0; s < 18; s++)
      td[s] = Double_t(e.total[s]);
    Double_t x = 0.0, y = 0.0;
    PlaneXY(td, reac, x, y);
    if (cutAn && Int_t(ev_an.size()) < kTracesPerRegion &&
        cutAn->IsInside(x, y))
      ev_an.push_back(&e);
    else if (cutAa && Int_t(ev_aa.size()) < kTracesPerRegion &&
             cutAa->IsInside(x, y))
      ev_aa.push_back(&e);
  }
 
  Double_t norm[18];
  for (Int_t s = 0; s < 18; s++)
    norm[s] = 1.0;
  for (Int_t s = 1; s <= 16; s++) {
    std::vector<Double_t> v;
    for (Int_t k = 0; k < Int_t(ev_beam.size()); k++) {
      Double_t sum =
          Double_t(ev_beam[k]->long_au[s]) + Double_t(ev_beam[k]->short_au[s]);
      if (sum > 0)
        v.push_back(sum);
    }
    if (v.size() < 20)
      continue;
    std::sort(v.begin(), v.end());
    Double_t med = v[v.size() / 2];
    if (med > 0)
      norm[s] = 1.0 / med;
  }
 
  std::vector<TGraph *> g_beam, g_aa, g_an;
  std::vector<const std::vector<const TraceEvt *> *> srcs;
  srcs.push_back(&ev_beam);
  srcs.push_back(&ev_aa);
  srcs.push_back(&ev_an);
  std::vector<std::vector<TGraph *> *> dsts;
  dsts.push_back(&g_beam);
  dsts.push_back(&g_aa);
  dsts.push_back(&g_an);
  for (Int_t c = 0; c < 3; c++) {
    for (Int_t k = 0; k < Int_t(srcs[c]->size()); k++) {
      const TraceEvt &e = *(*srcs[c])[k];
      Float_t alt[18];
      for (Int_t s = 0; s < 18; s++)
        alt[s] = e.long_au[s];
      for (Int_t s = 1; s <= 16; s++)
        alt[s] =
            long_only_is_current
                ? Float_t((Double_t(e.long_au[s]) + Double_t(e.short_au[s])) *
                          norm[s])
                : e.long_au[s];
      dsts[c]->push_back(TraceFromTotal(alt));
    }
  }
 
  const char *tag = long_only_is_current ? "sum" : "longonly";
  DrawRegionTraces(Form("region_traces_reac%d_altdecode_%s", reac, tag),
                   "strip_sum_scatter", g_beam, g_aa, g_an, 0.6, 1.6,
                   "#DeltaE [a.u.]");
  DrawRegionMeanTraces(
      Form("region_mean_traces_reac%d_altdecode_%s", reac, tag),
      "strip_sum_scatter", g_beam, g_aa, g_an, 0.6, 1.6, "#DeltaE [a.u.]");
 
  TString txt =
      Paths::ResultsDir() +
      Form("/plots/strip_sum_scatter/region_traces_reac%d_altdecode_%s.txt",
           reac, tag);
  std::ofstream out(txt.Data());
  if (out) {
    out << "# selected under "
        << (long_only_is_current ? "long-side-only" : "L+R sum")
        << " decode; alternate rendering is " << tag << std::endl;
    out << "# per-strip normalisation applied to the summed decode so that the"
           " summed BEAM median is 1.0 a.u."
        << std::endl;
    out << "# strip  norm  then per class: <long_only>  <sum>  ratio  "
           "short_frac"
        << std::endl;
    const char *cls[3] = {"beam", "aa", "an"};
    out << "strip norm";
    for (Int_t c = 0; c < 3; c++)
      out << " " << cls[c] << "_long " << cls[c] << "_sum " << cls[c]
          << "_ratio " << cls[c] << "_shortfrac";
    out << std::endl;
    for (Int_t s = 1; s <= 16; s++) {
      out << s << " " << Form("%.6f", norm[s]);
      for (Int_t c = 0; c < 3; c++) {
        Double_t sl = 0, ss = 0;
        Int_t nn = Int_t(srcs[c]->size());
        for (Int_t k = 0; k < nn; k++) {
          sl += Double_t((*srcs[c])[k]->long_au[s]);
          ss += Double_t((*srcs[c])[k]->short_au[s]);
        }
        if (nn > 0) {
          sl /= nn;
          ss /= nn;
        }
        Double_t sum = (sl + ss) * norm[s];
        out << " " << Form("%.5f", sl) << " " << Form("%.5f", sum) << " "
            << Form("%.5f", sl > 0 ? sum / sl : 0.0) << " "
            << Form("%.5f", (sl + ss) > 0 ? ss / (sl + ss) : 0.0);
      }
      out << std::endl;
    }
    out.close();
    std::cout << "Wrote " << txt << std::endl;
  }
 
  for (Int_t c = 0; c < 3; c++)
    for (Int_t k = 0; k < Int_t(dsts[c]->size()); k++)
      delete (*dsts[c])[k];
}
 
TGraph *StripSumScatter::TraceFromTotal(const Float_t *total) {
 
  Double_t td[18];
  for (Int_t s = 0; s < 18; s++)
    td[s] = Double_t(total[s]);
  return EventsSummary::BuildTraceFromTotals(td);
}
 
void StripSumScatter::DrawRegionTraces(const TString &save_name,
                                       const TString &subdir,
                                       const std::vector<TGraph *> &beam,
                                       const std::vector<TGraph *> &aa,
                                       const std::vector<TGraph *> &an,
                                       Double_t y_min, Double_t y_max,
                                       const char *y_title) {
  std::lock_guard<std::mutex> lock(g_plot_mutex);
  Int_t s_lo = Constants::cfg.IGNORE_STRIP_0 ? 1 : 0;
  Int_t s_hi = Constants::cfg.IGNORE_STRIP_17 ? 16 : 17;
  TH2F *frame =
      new TH2F("h_region_trace_frame", Form(";Strip;%s", y_title),
               s_hi - s_lo + 1, s_lo - 0.5, s_hi + 0.5, 100, y_min, y_max);
  frame->SetStats(0);
  TCanvas *c = PlottingUtils::GetConfiguredCanvas(kFALSE);
  frame->Draw();
  DrawTraceSet(beam, kGray + 2);
  DrawTraceSet(aa, kAzure + 2);
  DrawTraceSet(an, kRed + 1);
 
  TGraph *p_beam = new TGraph(1);
  TGraph *p_aa = new TGraph(1);
  TGraph *p_an = new TGraph(1);
  TGraph *proxies[3] = {p_beam, p_aa, p_an};
  Int_t pcol[3] = {kGray + 2, kAzure + 2, kRed + 1};
  for (Int_t i = 0; i < 3; i++) {
    proxies[i]->SetPoint(0, -1e9, -1e9);
    proxies[i]->SetLineColor(pcol[i]);
    proxies[i]->SetLineWidth(3);
  }
  TLegend *leg = PlottingUtils::AddLegend(0.725, 0.875, 0.70, 0.86);
  leg->AddEntry(p_beam, "Beam", "l");
  leg->AddEntry(p_aa, "(#alpha,#alpha')", "l");
  leg->AddEntry(p_an, "(#alpha,n)", "l");
  leg->Draw();
 
  PlottingUtils::SaveFigure(c, save_name, subdir, PlotSaveOptions::kLINEAR);
  delete c;
}
 
void StripSumScatter::DrawRegionMeanTraces(const TString &save_name,
                                           const TString &subdir,
                                           const std::vector<TGraph *> &beam,
                                           const std::vector<TGraph *> &aa,
                                           const std::vector<TGraph *> &an,
                                           Double_t y_min, Double_t y_max,
                                           const char *y_title) {
  std::lock_guard<std::mutex> lock(g_plot_mutex);
  TH2F *frame = new TH2F("h_region_mean_frame", Form(";Strip;%s", y_title), 18,
                         -0.5, 17.5, 100, y_min, y_max);
  frame->SetStats(0);
  frame->SetDirectory(nullptr);
  TCanvas *c = PlottingUtils::GetConfiguredCanvas(kFALSE);
  frame->Draw();
 
  const std::vector<TGraph *> *regions[3] = {&beam, &aa, &an};
  Int_t colors[3] = {kGray + 2, kAzure + 2, kRed + 1};
  const char *labels[3] = {"Beam", "(#alpha,#alpha')", "(#alpha,n)"};
  std::vector<TGraphErrors *> means;
  TLegend *leg = PlottingUtils::AddLegend(0.725, 0.875, 0.70, 0.86);
 
  for (Int_t r = 0; r < 3; r++) {
    const std::vector<TGraph *> &tr = *regions[r];
    if (tr.empty())
      continue;
    Int_t npts = tr[0]->GetN();
    std::vector<Double_t> mean(npts, 0.0), m2(npts, 0.0);
    for (Int_t t = 0; t < Int_t(tr.size()); t++) {
      Double_t *yv = tr[t]->GetY();
      for (Int_t p = 0; p < npts; p++) {
        mean[p] += yv[p];
        m2[p] += yv[p] * yv[p];
      }
    }
    Double_t nt = Double_t(tr.size());
    TGraphErrors *ge = new TGraphErrors(npts);
    Double_t *xv = tr[0]->GetX();
    for (Int_t p = 0; p < npts; p++) {
      mean[p] /= nt;
      Double_t var = m2[p] / nt - mean[p] * mean[p];
      ge->SetPoint(p, xv[p], mean[p]);
      ge->SetPointError(p, 0.0, var > 0.0 ? TMath::Sqrt(var) : 0.0);
    }
    ge->SetLineColor(colors[r]);
    ge->SetLineWidth(3);
    ge->SetFillColorAlpha(colors[r], 0.15);
    ge->Draw("3 SAME"); // +-1 RMS band (all bands first, behind the lines)
    means.push_back(ge);
    leg->AddEntry(ge, labels[r], "l");
  }
  // Mean lines on top of every band.
  for (Int_t i = 0; i < Int_t(means.size()); i++)
    means[i]->Draw("LX SAME"); // mean line, no end caps
  leg->Draw();
 
  PlottingUtils::SaveFigure(c, save_name, subdir, PlotSaveOptions::kLINEAR);
  for (Int_t i = 0; i < Int_t(means.size()); i++)
    delete means[i];
  delete leg;
  delete c;
  delete frame;
}
 
void StripSumScatter::TraceYRange(const std::vector<TGraph *> &beam,
                                  const std::vector<TGraph *> &aa,
                                  const std::vector<TGraph *> &an,
                                  Double_t &y_min, Double_t &y_max) {
  y_min = std::numeric_limits<Double_t>::max();
  y_max = -std::numeric_limits<Double_t>::max();
  const std::vector<TGraph *> *sets[3] = {&beam, &aa, &an};
  for (Int_t si = 0; si < 3; si++) {
    const std::vector<TGraph *> &v = *sets[si];
    for (Int_t i = 0; i < Int_t(v.size()); i++) {
      Double_t x = 0.0, y = 0.0;
      for (Int_t k = 0; k < v[i]->GetN(); k++) {
        v[i]->GetPoint(k, x, y);
        if (x < 0.5 || x > 16.5)
          continue;
        if (y < y_min)
          y_min = y;
        if (y > y_max)
          y_max = y;
      }
    }
  }
  if (y_min > y_max) { // no in-range points sampled
    y_min = 0.0;
    y_max = 1.0;
  }
  Double_t pad = 0.05 * (y_max - y_min);
  if (pad <= 0.0)
    pad = 1.0;
  y_min -= pad;
  y_max += pad;
}
 
TCutG *StripSumScatter::PromptCut(TCanvas *c, const char *name,
                                  const char *label) {
  std::cout << "  >>> draw the " << label
            << " region: left-click vertices, double-click to close"
            << std::endl;
  c->cd();
  TCutG *cut = static_cast<TCutG *>(c->WaitPrimitive("CUTG", "CutG"));
  if (!cut) {
    std::cerr << "  no " << label << " cut drawn" << std::endl;
    return nullptr;
  }
  cut->SetName(name);
  cut->SetLineColor(kBlack);
  cut->SetLineWidth(2);
  return cut;
}
 
// Saved region cuts.
//
// Drawing the regions by hand is the only way to place them when the reaction
// population's location is not yet known, but a hand-drawn cut that is not
// stored makes the run unreproducible: the next pass gets a different polygon
// and no two results are comparable. Persisting them separates "decide where
// the region is", which needs a person once, from "apply it", which should be
// automatic from then on -- and without a DISPLAY.
void StripSumScatter::SaveRegionCuts(Int_t reac, TCutG *cut_an, TCutG *cut_aa) {
  RegionCutStore::Save(reac, cut_an, cut_aa);
}
 
TCutG *StripSumScatter::LoadRegionCut(const char *name, Int_t reac) {
  return RegionCutStore::Load(name, reac);
}
 
void StripSumScatter::SmoothTrace(const Double_t *in, Double_t *out,
                                  Int_t width) {
  Int_t half = width / 2;
  for (Int_t s = 0; s < 18; s++) {
    Int_t lo = TMath::Max(0, s - half);
    Int_t hi = TMath::Min(17, s + half);
    Double_t sum = 0.0;
    for (Int_t t = lo; t <= hi; t++)
      sum += in[t];
    out[s] = sum / Double_t(hi - lo + 1);
  }
}
 
// Savitzky-Golay smoothing: 3rd-degree polynomial, half-window of 2
// (5-point convolution). Uses standard SG coefficients that sum to 1.
// For a 5-point window with 3rd-degree polynomial, the smoothed value at
// the center uses coefficients: [-3, 12, 17, 12, -3] / 35.
// At edges, the window shrinks and coefficients are renormalised.
void StripSumScatter::SavitzkyGolay(const Double_t *in, Double_t *out) {
  static const Int_t K = 2; // half-width (5-point window)
 
  // Standard SG coefficients for 5-point, 3rd-degree polynomial (smoothed
  // value): These are translation-invariant - same for all center positions.
  static const Double_t sg_coeff[2 * K + 1] = {
      -3.0 / 35.0, // coefficient for t = s - 2
      12.0 / 35.0, // coefficient for t = s - 1
      17.0 / 35.0, // coefficient for t = s (center)
      12.0 / 35.0, // coefficient for t = s + 1
      -3.0 / 35.0  // coefficient for t = s + 2
  };
 
  for (Int_t s = 0; s < 18; s++) {
    Int_t lo = TMath::Max(0, s - K);
    Int_t hi = TMath::Min(17, s + K);
    Double_t val = 0.0;
    Double_t wsum = 0.0;
 
    // Apply SG coefficients for the positions within the clipped window
    for (Int_t t = lo; t <= hi; t++) {
      Int_t offset = t - s + K; // 0..4, position within 5-point window
      val += sg_coeff[offset] * in[t];
      wsum += sg_coeff[offset];
    }
 
    // Renormalise at edges where window shrinks
    out[s] = (wsum != 0.0) ? val / wsum : in[s];
  }
}
 
// CFD-style trigger finder: scan left-to-right for the first strip whose
// beam-subtracted signal exceeds both a fraction of the trace peak and a
// multiple of the beam sigma. Returns the strip index, or -1 if none fires.
 
Int_t StripSumScatter::FindTrigger(const Double_t *td, const Double_t *base,
                                   Double_t beam_sigma) {
  const Int_t s_lo = 2;
  const Int_t s_hi = 16;
 
  const Double_t frac =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRIGGER_CFD_FRAC;
 
  const Double_t nsigma =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRIGGER_NSIGMA * beam_sigma;
 
  Double_t peak_signal = -1.0e30;
 
  for (Int_t s = s_lo; s <= s_hi; s++) {
    Double_t signal = td[s] - base[s];
    if (signal > peak_signal)
      peak_signal = signal;
  }
 
  if (peak_signal < nsigma)
    return -1;
 
  Double_t thresh = frac * peak_signal;
 
  for (Int_t s = s_lo; s <= s_hi; s++) {
    Double_t signal = td[s] - base[s];
 
    if (signal >= thresh && signal >= nsigma)
      return s;
  }
 
  return -1;
}
 
// Build a TGraph trace from a Savitzky-Golay-smoothed set of per-strip totals.
TGraph *StripSumScatter::SmoothedTraceFromTotal(const Float_t *total) {
  Double_t td[18], sgd[18];
  for (Int_t s = 0; s < 18; s++)
    td[s] = Double_t(total[s]);
  SavitzkyGolay(td, sgd);
  return EventsSummary::BuildTraceFromTotals(sgd);
}
 
void StripSumScatter::ClusterVarHists(Int_t reac, TCutG *cut_aa, TCutG *cut_an,
                                      const TString &subdir) {
  const Int_t NV = 9;
  const Int_t NC = 3;
  const char *vkey[NV] = {"energy",      "peak3",      "plateau",
                          "tail",        "reacstrip",  "mult",
                          "trigtaildev", "reacslope3", "beamdev"};
  const char *vtitle[NV] = {
      "#Sigma_{all strips}(#DeltaE#minus1) [a.u.]",
      "#Sigma_{trig#pm1}#DeltaE (0 if no trigger) [a.u.]",
      "Plateau Excess  #Sigma_{trig+1..trig+POST}(#DeltaE#minus1) [a.u.]",
      "#DeltaE(s17) [a.u.]",
      "Trigger Strip",
      "Both-side Multiplicity (strips 1-16)",
      "|#DeltaE#minusbeam| at trigger + at s17 [a.u.]",
      "#DeltaE(reac+3) #minus #DeltaE(reac#minus3) [a.u.]",
      "RMS_{8-17}(#DeltaE#minusbeam) [a.u.]"};
  const char *clabel[NC] = {"beam", "(a,a')", "(a,n)"};
 
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kClusterSmoothWindow =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.CLUSTER_SMOOTH_WINDOW;
 
  // Per (class, variable) value lists for raw traces.
  std::vector<Double_t> vals_raw[NC][NV];
  // Same structure, but cluster variables computed on Savitzky-Golay-smoothed
  // traces (SG kernel removes the L_odd/R_even sawtooth so peak/onset features
  // land on the real physical peak rather than an odd strip).
  std::vector<Double_t> vals_sg[NC][NV];
  UInt_t bit = (1u << ReacIndex(reac));
 
  // Per-strip beam baseline = mean over the beam-flat reservoir events. It
  // carries the L_odd/R_even sawtooth, so subtracting it removes that
  // systematic exactly before the reaction-onset jump search (better than
  // blurring it with smoothing).
  Double_t base[18];
  for (Int_t s = 0; s < 18; s++)
    base[s] = 0.0;
  Long64_t nbeam = 0;
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++)
    if (m_reservoir[k].beam_flat) {
      for (Int_t s = 0; s < 18; s++)
        base[s] += Double_t(m_reservoir[k].total[s]);
      nbeam++;
    }
  for (Int_t s = 0; s < 18; s++)
    base[s] = (nbeam > 0) ? base[s] / Double_t(nbeam) : 1.0;
 
  // Pooled beam-noise RMS = sqrt(mean over beam-flat events and all strips of
  // (total - base)^2). The reaction-onset threshold is this many sigma (an
  // N-sigma discriminator), matching the Python pipeline -- so on flat beam
  // the excess does NOT cross on average and the event gets NO trigger.
  Double_t beam_sumsq = 0.0;
  Long64_t beam_npt = 0;
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++)
    if (m_reservoir[k].beam_flat)
      for (Int_t s = 0; s < 18; s++) {
        Double_t d = Double_t(m_reservoir[k].total[s]) - base[s];
        beam_sumsq += d * d;
        beam_npt++;
      }
  Double_t beam_sigma =
      (beam_npt > 0) ? TMath::Sqrt(beam_sumsq / beam_npt) : 0.0;
 
  // Count triggers over ALL reservoir events (not just classified subset) so
  // the numbers are directly comparable to the Python pipeline.
  Long64_t triggered = 0, no_trigger = 0;
  Double_t td_all[18];
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    for (Int_t s = 0; s < 18; s++)
      td_all[s] = Double_t(m_reservoir[k].total[s]);
    if (FindTrigger(td_all, base, beam_sigma) >= 0)
      triggered++;
    else
      no_trigger++;
  }
  const Double_t reac_onset_gate =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRIGGER_NSIGMA * beam_sigma;
  const Double_t cf_frac =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRIGGER_CFD_FRAC;
  std::cout << "  beam reference: mean+RMS of " << nbeam
            << " pure-beam events (fitted s0,s1 & s16,s17 ellipses); "
            << Form("noise sigma=%.4f", beam_sigma) << std::endl;
  std::cout << Form("  reaction onset: gate %g-sigma = %.4f, CF fraction %g; "
                    "triggered %lld of %lld (no trigger: %lld)",
                    Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRIGGER_NSIGMA,
                    reac_onset_gate, cf_frac, Long64_t(triggered),
                    Long64_t(m_reservoir.size()), Long64_t(no_trigger))
            << std::endl;
 
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    const TraceEvt &e = m_reservoir[k];
    Double_t td[18];
    for (Int_t s = 0; s < 18; s++)
      td[s] = Double_t(e.total[s]);
    Int_t cls = -1;
    if (e.beam_flat)
      cls = 0;
    else if (e.reac_mask & bit) {
      Double_t x = 0.0, y = 0.0;
      PlaneXY(td, reac, x, y);
      if (cut_aa && cut_aa->IsInside(x, y))
        cls = 1;
      else if (cut_an && cut_an->IsInside(x, y))
        cls = 2;
    }
    if (cls < 0)
      continue;
    // The same five variables the blind clustering uses (normed total, beam
    // at 1 per strip; guards 0/17 included).
    Double_t energy = 0.0;
    for (Int_t s = 0; s < 18; s++) {
      energy += td[s] - 1.0;
    }
 
    Int_t trigger_strip = FindTrigger(td, base, beam_sigma);
    Bool_t has_trig = (trigger_strip >= 0);
    // Plateau excess: beam-subtracted sum over a sliding window that tracks
    // the trigger (reacstrip+1 .. reacstrip+PLATEAU_POST). 0 when no trigger.
    // Out-of-range strips are dropped from the sum.
    const Int_t kPlateauPost =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PLATEAU_POST;
    Double_t plateau = 0.0;
    if (has_trig)
      for (Int_t d = 1; d <= kPlateauPost; d++) {
        Int_t s = trigger_strip + d;
        if (s >= 0 && s < 18)
          plateau += td[s] - 1.0;
      }
    // peak-3 sum centered on the TRIGGER strip (not the argmax). 0 when there
    // is no trigger, which is the discriminator: flat beam scores 0.
    Double_t peak3 = 0.0;
    if (has_trig)
      for (Int_t s = trigger_strip - 1; s <= trigger_strip + 1; s++)
        if (s >= 0 && s < 18)
          peak3 += td[s];
    // |deviation from beam| at the TRIGGER strip plus at the end strip s17 --
    // the (a,n) signature is a rise at the trigger AND a collapse at s17.
    Double_t trigtaildev =
        has_trig ? TMath::Abs(td[trigger_strip] - base[trigger_strip]) +
                       TMath::Abs(td[17] - base[17])
                 : 0.0;
    // Slopes ACROSS the trigger: dE(reac+n) - dE(reac-n). reac+-n share parity
    // so the L_odd/R_even sawtooth cancels. Only valid (and filled) when there
    // is a trigger AND both endpoints are in range (symmetric, no clamping).
    Bool_t ok3 =
        has_trig && (trigger_strip - 3 >= 0) && (trigger_strip + 3 <= 17);
    Double_t reacslope3 =
        ok3 ? td[trigger_strip + 3] - td[trigger_strip - 3] : 0.0;
 
    // How beam-like the BACK HALF (strips 8-17) is: RMS deviation of the trace
    // from the beam baseline over those strips. Subtracting base[] removes the
    // L/R sawtooth, so this is clean; LOW = beam-like (flat at beam level),
    // high for a reaction's plateau/collapse or the elevation of pileup.
    // Amplitude-aware (NOT max-normalized like the template-prune residual):
    // the beam has a fixed level, so a flat-but-elevated trace must not read as
    // beam. No trigger needed (fixed window), so always filled.
    Double_t beamdev = 0.0;
    Int_t n_bl = 0;
    for (Int_t s = 8; s <= 17; s++) {
      Double_t d = td[s] - base[s];
      beamdev += d * d;
      n_bl++;
    }
    beamdev = TMath::Sqrt(beamdev / Double_t(n_bl));
    Double_t v[NV] = {energy,
                      peak3,
                      plateau,
                      td[17], // raw end strip (was td[17] - 1.0)
                      Double_t(trigger_strip),
                      Double_t(e.both_mult),
                      trigtaildev,
                      reacslope3,
                      beamdev};
    // peak3 is filled even with no trigger (it scores 0 -- the discriminator).
    // The other trigger-centered vars are skipped when there is no trigger
    // (reacstrip) or the symmetric window runs off an edge (slopes /
    // jaggedness).
    Bool_t vok[NV] = {kTRUE, kTRUE,    kTRUE, kTRUE, kTRUE,
                      kTRUE, has_trig, ok3,   kTRUE};
    for (Int_t iv = 0; iv < NV; iv++)
      if (vok[iv])
        vals_raw[cls][iv].push_back(v[iv]);
 
    // Savitzky-Golay-smoothed trace: recompute trigger and cluster variables on
    // the SG-filtered copy. This removes the L_odd/R_even sawtooth so onset /
    // peak / slope features land on the real physical structure.
    Double_t sgd[18];
    SavitzkyGolay(td, sgd);
    Double_t energy_sg = 0.0;
    for (Int_t s = 0; s < 18; s++)
      energy_sg += sgd[s] - 1.0;
    Double_t ex_sg[18], sm_ex_sg[18];
    for (Int_t s = 0; s < 18; s++)
      ex_sg[s] = sgd[s] - base[s];
    SmoothTrace(ex_sg, sm_ex_sg, kClusterSmoothWindow);
    Int_t reacstrip_sg = FindTrigger(sgd, base, beam_sigma);
    Bool_t has_trig_sg = (reacstrip_sg >= 0);
    // Plateau excess on SG trace: sliding window over reacstrip_sg+1 ..
    // reacstrip_sg+PLATEAU_POST. 0 when no trigger; out-of-range strips
    // dropped.
    Double_t plateau_sg = 0.0;
    if (has_trig_sg)
      for (Int_t d = 1; d <= kPlateauPost; d++) {
        Int_t s = reacstrip_sg + d;
        if (s >= 0 && s < 18)
          plateau_sg += sgd[s] - 1.0;
      }
    Double_t peak3_sg = 0.0;
    if (has_trig_sg)
      for (Int_t s = reacstrip_sg - 1; s <= reacstrip_sg + 1; s++)
        if (s >= 0 && s < 18)
          peak3_sg += sgd[s];
    Double_t trigtaildev_sg =
        has_trig_sg ? TMath::Abs(sgd[reacstrip_sg] - base[reacstrip_sg]) +
                          TMath::Abs(sgd[17] - base[17])
                    : 0.0;
    Bool_t ok3_sg =
        has_trig_sg && (reacstrip_sg - 3 >= 0) && (reacstrip_sg + 3 <= 17);
    Double_t reacslope3_sg =
        ok3_sg ? sgd[reacstrip_sg + 3] - sgd[reacstrip_sg - 3] : 0.0;
    Double_t beamdev_sg = 0.0;
    Int_t n_bl_sg = 0;
    for (Int_t s = 8; s <= 17; s++) {
      Double_t d = sgd[s] - base[s];
      beamdev_sg += d * d;
      n_bl_sg++;
    }
    beamdev_sg = TMath::Sqrt(beamdev_sg / Double_t(n_bl_sg));
    Double_t v_sg[NV] = {energy_sg,
                         peak3_sg,
                         plateau_sg,
                         sgd[17],
                         Double_t(reacstrip_sg),
                         Double_t(e.both_mult),
                         trigtaildev_sg,
                         reacslope3_sg,
                         beamdev_sg};
    Bool_t vok_sg[NV] = {kTRUE, kTRUE,       kTRUE,  kTRUE, kTRUE,
                         kTRUE, has_trig_sg, ok3_sg, kTRUE};
    for (Int_t iv = 0; iv < NV; iv++)
      if (vok_sg[iv])
        vals_sg[cls][iv].push_back(v_sg[iv]);
  }
 
  std::cout << "cluster-var hists (reac " << reac
            << "): beam=" << vals_raw[0][0].size()
            << " (a,a')=" << vals_raw[1][0].size()
            << " (a,n)=" << vals_raw[2][0].size() << std::endl;
 
  std::vector<Int_t> colors = PlottingUtils::GetDefaultColors();
  // Two smoothing passes: raw traces, then Savitzky-Golay smoothed.
  const Int_t kNP = 2;
  const char *pass_label[kNP] = {"raw", "sg"};
 
  for (Int_t ip = 0; ip < kNP; ip++) {
    if (ip == 1 && Constants::cfg.STRIP_SUM_SCATTER_CONFIG.SKIP_SAVGOL_PLOTS)
      continue;
    std::vector<Double_t>(*vals)[NC][NV] = (ip == 0) ? &vals_raw : &vals_sg;
    for (Int_t iv = 0; iv < NV; iv++) {
      Double_t lo = 1.0e30, hi = -1.0e30;
      Double_t mean[NC] = {0.0, 0.0, 0.0};
      for (Int_t ic = 0; ic < NC; ic++) {
        for (Int_t m = 0; m < Int_t((*vals)[ic][iv].size()); m++) {
          lo = TMath::Min(lo, (*vals)[ic][iv][m]);
          hi = TMath::Max(hi, (*vals)[ic][iv][m]);
          mean[ic] += (*vals)[ic][iv][m];
        }
        if (!(*vals)[ic][iv].empty())
          mean[ic] /= Double_t((*vals)[ic][iv].size());
      }
      if (ip == 0) {
        std::cout << "  [" << pass_label[ip] << "] " << vkey[iv]
                  << ": mean beam=" << mean[0] << " (a,a')=" << mean[1]
                  << " (a,n)=" << mean[2] << " [range " << lo << ".." << hi
                  << "]" << std::endl;
      }
      Int_t nbins = 80;
      if (iv == 4) { // reaction strip: integer bins 0..17
        lo = -1.5;
        hi = 17.5;
        nbins = 19;
      } else if (iv == 5) { // both-channel multiplicity: integer bins 0..16
        lo = -0.5;
        hi = 16.5;
        nbins = 17;
      } else {
        if (!(hi > lo)) { // constant -> give it a drawable range
          lo -= 0.5;
          hi += 0.5;
        }
        Double_t pad = 0.05 * (hi - lo);
        lo -= pad;
        hi += pad;
      }
 
      TCanvas *c = PlottingUtils::GetConfiguredCanvas(kTRUE);
      TString axis = Form(";%s;Counts", vtitle[iv]);
      std::vector<TH1F *> hs;
      Double_t ymax = 0.0;
      for (Int_t ic = 0; ic < NC; ic++) {
        TH1F *h = new TH1F(
            Form("h_cv_%s_%s_c%d_r%d", vkey[iv], pass_label[ip], ic, reac),
            axis, nbins, lo, hi);
        h->SetDirectory(nullptr);
        for (Int_t m = 0; m < Int_t((*vals)[ic][iv].size()); m++)
          h->Fill((*vals)[ic][iv][m]);
 
        PlottingUtils::ConfigureHistogram(h, colors[ic % Int_t(colors.size())],
                                          axis);
        h->SetStats(0);
        ymax = TMath::Max(ymax, h->GetMaximum());
        hs.push_back(h);
      }
      if (ymax <= 0.0)
        ymax = 1.0;
      TLegend *leg = PlottingUtils::AddLegend(0.775, 0.875, 0.70, 0.86);
 
      for (Int_t ic = 0; ic < NC; ic++) {
        if (ic == 0) {
          hs[0]->SetMaximum(3.0 * ymax);
          hs[0]->SetMinimum(1.0e-1);
          hs[0]->Draw("HIST");
        } else {
          hs[ic]->Draw("HIST SAME");
        }
        leg->AddEntry(hs[ic], clabel[ic], "l");
      }
      leg->Draw();
 
      TString sub_subdir = subdir + "/clusters_" + pass_label[ip];
 
      PlottingUtils::SaveFigure(
          c, Form("cluster_var_%s_%s_reac%d", vkey[iv], pass_label[ip], reac),
          sub_subdir, PlotSaveOptions::kLOG);
      for (Int_t m = 0; m < Int_t(hs.size()); m++)
        delete hs[m];
      delete c;
    }
  }
}
 
TString StripSumScatter::BuildFingerprint(const std::vector<Int_t> &run_order,
                                          std::map<Int_t, TChain *> &chains) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Double_t kReacJumpNSigma =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REAC_JUMP_NSIGMA;
  const Double_t kReacJumpMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REAC_JUMP_MAX;
  const Int_t kSmoothHiStrip =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS_END_STRIP;
  const Double_t kSmoothMaxStep =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS_MAX_STEP;
  const Double_t kEndStripMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.END_STRIP_MAX;
 
  const Int_t kGateStripX =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_STRIP_X;
  const Int_t kGateStripY =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_STRIP_Y;
  const Double_t kGateNSigmaX =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_X;
  const Double_t kGateNSigmaY =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_NSIGMA_Y;
  const Int_t kGateBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_BINS;
  const Double_t kGateMin = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_MIN;
  const Double_t kGateMax = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.GATE_MAX;
 
  const Int_t kXBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.XBINS;
  const Int_t kYBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.YBINS;
 
  // Two parts. Before the bar: everything that decides which events are
  // tagged and kept, so a change there means a pass over the events files.
  // After it: only what is built from the tagged events, all of which the
  // reservoir keeps, so a cache whose tagging matches but whose plane does
  // not is re-projected from its reservoir in a minute rather than refilled.
  TString s = Form(
      "v19 reac[%d,%d] bmult[%d,%d] jump[%.2fsig,%.3f] smooth=%d,%d "
      "step=%.3f s17=%.3f gate[s%d,s%d,%.2f,%.2f,%d,%.3f,%.3f] par=%.3f",
      kReacMin, kReacMax, Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BOTH_MULT_MAX,
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BOTH_MULT_COUNT_TO,
      kReacJumpNSigma, kReacJumpMax,
      Int_t(Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REQUIRE_SMOOTHNESS),
      kSmoothHiStrip, kSmoothMaxStep, kEndStripMax, kGateStripX, kGateStripY,
      kGateNSigmaX, kGateNSigmaY, kGateBins, kGateMin, kGateMax,
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PARITY_ASYM_MAX);
  // The gate resolves through the measured noise, so the thresholds actually
  // applied are stamped too: a drift in the noise refills.
  s += " jmin[";
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
    s += Form("%s%.4f", reac == kReacMin ? "" : ",", JumpMin(reac));
  s += "]";
  // The upstream-beam precondition changes which events are tagged; like the
  // jump gate it resolves through measured noise, so the tolerances are
  // stamped.
  s += Form(" up=%d,%.2fsig[",
            Int_t(Constants::cfg.STRIP_SUM_SCATTER_CONFIG
                      .REQUIRE_BEAM_UPSTREAM_OF_REAC),
            Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BEAM_UPSTREAM_NSIGMA);
  for (Int_t strip = 1; strip < kReacMax; strip++)
    s += Form("%s%.4f", strip == 1 ? "" : ",",
              Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BEAM_UPSTREAM_NSIGMA *
                  StripSigma(strip));
  s += "]";
  // Active beam gates (also keyed by cache filename, but folded in here too so
  // a mismatch never silently reuses a stale same-named cache).
  std::vector<GateSpec> gates = ActiveGates();
  for (Int_t i = 0; i < Int_t(gates.size()); i++)
    s += Form(" g[s%d,s%d]", gates[i].sx, gates[i].sy);
 
  // Display windows are deliberately absent: they no longer change what is
  // built, so retuning them must not invalidate the cache.
  for (Int_t i = 0; i < Int_t(run_order.size()); i++) {
    Int_t run = run_order[i];
    s += Form(" r%d:%lld", run, chains[run]->GetEntries());
  }
 
  // The built quantity: build range, binning, the x range and each strip's
  // y window, so a change to the window rule re-projects.
  TString plane =
      Form("buildx[%.3f,%.3f] buildy[%.3f,%.3f] bins[%d,%d] x[%d,%d]",
           ScatterBuildRange::kXMin, ScatterBuildRange::kXMax,
           ScatterBuildRange::kYMin, ScatterBuildRange::kYMax, kXBins, kYBins,
           Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO,
           Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI);
  plane += " y";
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
    plane +=
        Form("%s%d-%d", reac == kReacMin ? "[" : ",", YLoOf(reac), YHiOf(reac));
  plane += "]";
  return s + " | " + plane;
}
 
// The tagging half of a fingerprint; a cache from before the split has no
// bar and never matches.
static TString TagPart(const TString &fingerprint) {
  const Ssiz_t bar = fingerprint.Index(" | ");
  return bar < 0 ? TString("") : TString(fingerprint(0, bar));
}
 
// Per-reaction-strip y-axis bounds straight from
// Constants::cfg.STRIP_SUM_SCATTER_CONFIG.Y_DISPLAY_RANGE (tunable per dataset,
// per strip); strips absent from the map fall back to
// Y_DISPLAY_MIN/Y_DISPLAY_MAX. x stays fixed (strip-independent).
void StripSumScatter::YBounds(Double_t *y_lo, Double_t *y_hi) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
    Int_t ri = reac - kReacMin;
    std::map<Int_t, std::pair<Double_t, Double_t>>::const_iterator it =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.Y_DISPLAY_RANGE.find(reac);
    if (it != Constants::cfg.STRIP_SUM_SCATTER_CONFIG.Y_DISPLAY_RANGE.end()) {
      y_lo[ri] = it->second.first;
      y_hi[ri] = it->second.second;
    } else {
      y_lo[ri] = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.Y_DISPLAY_MIN;
      y_hi[ri] = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.Y_DISPLAY_MAX;
    }
  }
}
 
TString StripSumScatter::PrettyLabel(const TString &tag) {
  TString base = RemixSim::TagWithoutStrip(tag);
  base.ReplaceAll("_eres", "");
  if (base == "aa")
    return "(#alpha,#alpha')";
  if (base == "an")
    return "(#alpha,n)";
  if (base == "beam")
    return "Beam";
  return base;
}
 
// Per-strip sim normalization gains: read the sim beam file and, for each
// strip, average the (unit-gain) per-strip beam total, then set gain[s] = 1 /
// mean[s] so every strip's sim beam lands on 1 a.u. This flattens the sim beam
// the SAME way the per-channel data normalization flattens the experimental
// beam
// -- a single global factor would not, since it preserves the sim's per-strip
// structure. Strips with no beam signal keep gain 0 (drop out like an
// uncalibrated channel).
Bool_t StripSumScatter::SimBeamGains(Double_t *gain) {
  const Long64_t kSampleMaxPoints =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.SAMPLE_MAX_POINTS;
 
  for (Int_t s = 0; s < 18; s++)
    gain[s] = 0.0;
  // Reference the ERES beam file -- the same file type as the eres populations
  // plotted in SimOverlay/SimTraceOverlay -- so the normalized eres beam lands
  // exactly on 1. Falls back to the non-eres SIM_BEAM_FILE if no eres beam
  // control file is present.
  TString file;
  std::vector<RemixSim::SimFileSpec> specs = RemixSim::BuildFileSpecs();
  for (Int_t i = 0; i < Int_t(specs.size()); i++) {
    if (!RemixSim::IsEresTag(specs[i].tag))
      std::cout << "No eres sim file for " << specs[i].tag << ", using standard"
                << std::endl;
 
    TString base = RemixSim::TagWithoutStrip(specs[i].tag);
    base.ReplaceAll("_eres", "");
    if (base == "beam") {
      file = RemixSim::SimRootPath(specs[i]);
      break;
    }
  }
  if (file.Length() == 0)
    file =
        Paths::DatasetDir() + "/sim_root_files/" + Constants::cfg.SIM_BEAM_FILE;
  TFile *f = IO::OpenForReading(file);
  if (!f || f->IsZombie()) {
    std::cerr << "strip-sum-scatter: cannot open sim beam file " << file
              << "; sim overlay stays in raw sim units." << std::endl;
    if (f)
      delete f;
    return kFALSE;
  }
  TTree *t = static_cast<TTree *>(f->Get("events_MeV"));
  if (!t) {
    std::cerr << "strip-sum-scatter: no events_MeV tree in sim beam file "
              << file << "; sim overlay stays in raw sim units." << std::endl;
    f->Close();
    delete f;
    return kFALSE;
  }
  Float_t left[18] = {0}, right[18] = {0};
  t->SetBranchAddress("Left_0_17_dE", left);
  t->SetBranchAddress("RightdE", right);
  Long64_t n = t->GetEntries();
  Long64_t stride = FileSet::SampleStride(n, kSampleMaxPoints);
  Double_t sum[18] = {0};
  Long64_t cnt[18] = {0};
  // Unit gains so SimTotal yields the raw per-strip beam total (IGNORE_SHORT
  // aware), which is exactly the quantity these gains will later normalize.
  Double_t unit[18];
  for (Int_t s = 0; s < 18; s++)
    unit[s] = 1.0;
  for (Long64_t j = 0; j < n; j += stride) {
    t->GetEntry(j);
    Double_t total[18];
    SimTotal(left, right, unit, total);
    for (Int_t s = 0; s < 18; s++)
      if (total[s] > 0.0) {
        sum[s] += total[s];
        cnt[s]++;
      }
  }
  f->Close();
  delete f;
  Int_t n_set = 0;
  for (Int_t s = 0; s < 18; s++) {
    if (cnt[s] > 0 && sum[s] > 0.0) {
      gain[s] = 1.0 / (sum[s] / Double_t(cnt[s]));
      n_set++;
    }
  }
  if (n_set == 0)
    return kFALSE;
  std::cout << "strip-sum-scatter: sim per-strip beam normalization to 1 a.u. ("
            << n_set << " strips)." << std::endl;
  return kTRUE;
}
 
void StripSumScatter::SimTotal(const Float_t *left, const Float_t *right,
                               const Double_t *gain, Double_t *total) {
  for (Int_t s = 0; s < 18; s++)
    total[s] = gain[s] * (Double_t(left[s]) + Double_t(right[s]));
  if (Constants::cfg.IGNORE_SHORT_STRIPS)
    for (Int_t s = 1; s <= 16; s++)
      total[s] =
          gain[s] * ((s % 2) != 0 ? Double_t(left[s]) : Double_t(right[s]));
}
 
TGraph *StripSumScatter::SimPopScatter(const TString &file, Int_t reac,
                                       const Double_t *gain,
                                       Long64_t max_points) {
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
 
  TFile *f = IO::OpenForReading(file);
  if (!f || f->IsZombie()) {
    if (f)
      delete f;
    return nullptr;
  }
  TTree *t = static_cast<TTree *>(f->Get("events_MeV"));
  if (!t) {
    std::cerr << "  no events_MeV tree in " << file << std::endl;
    f->Close();
    delete f;
    return nullptr;
  }
  Int_t y_lo = reac + 1;
  Int_t y_hi = TMath::Min(reac + 6, 17);
  Float_t left[18] = {0}, right[18] = {0};
  t->SetBranchAddress("Left_0_17_dE", left);
  t->SetBranchAddress("RightdE", right);
  Long64_t n = t->GetEntries();
  Long64_t stride = FileSet::SampleStride(n, max_points);
  TGraph *g = new TGraph();
  Long64_t k = 0;
  for (Long64_t j = 0; j < n; j += stride) {
    t->GetEntry(j);
    Double_t total[18];
    SimTotal(left, right, gain, total);
    Double_t x = SumRange(total, kXLo, kXHi);
    Double_t y = SumRange(total, y_lo, y_hi);
    if (x > 0.0)
      g->SetPoint(k++, x, y);
  }
  g->Set(k);
  f->Close();
  delete f;
  return g;
}
 
// Up to max_traces per-strip dE-profile traces, stride-sampled across one sim
// file. Sim energies are arbitrary-unit floats; total[s] is the per-strip
// normalized gain[s]*(left[s] + right[s]) so the traces share the data's axis
// (and the sim beam is flat at NORM, like the data).
std::vector<TGraph *> StripSumScatter::SimPopTraces(const TString &file,
                                                    const Double_t *gain,
                                                    Long64_t max_traces) {
  std::vector<TGraph *> traces;
  TFile *f = IO::OpenForReading(file);
  if (!f || f->IsZombie()) {
    if (f)
      delete f;
    return traces;
  }
  TTree *t = static_cast<TTree *>(f->Get("events_MeV"));
  if (!t) {
    f->Close();
    delete f;
    return traces;
  }
  Float_t left[18] = {0}, right[18] = {0};
  t->SetBranchAddress("Left_0_17_dE", left);
  t->SetBranchAddress("RightdE", right);
  Long64_t n = t->GetEntries();
  Long64_t stride = FileSet::SampleStride(n, max_traces);
  for (Long64_t j = 0; j < n && Int_t(traces.size()) < max_traces;
       j += stride) {
    t->GetEntry(j);
    Double_t total[18];
    SimTotal(left, right, gain, total);
    traces.push_back(EventsSummary::BuildTraceFromTotals(total));
  }
  f->Close();
  delete f;
  return traces;
}
 
// Per reaction strip, overlay TRACES_PER_CLASS sampled per-strip traces of each
// sim population in the experimental DrawRegionTraces style (beam grey, (a,a')
// azure, (a,n) red). The beam reference is the same for every strip. Sampled
// fresh each run (40 traces/file is trivial).
void StripSumScatter::SimTraceOverlay() {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kTracesPerRegion =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRACES_PER_CLASS;
 
  std::vector<RemixSim::SimFileSpec> specs = RemixSim::BuildFileSpecs();
  if (specs.empty())
    return;
  // One file per class and strip, the _eres twin preferred where it exists
  // (alphabetical order would otherwise leave the plain file in the map).
  std::map<Int_t, TString> aa_file, an_file; // reaction strip -> sim file
  std::map<Int_t, Bool_t> aa_eres, an_eres;
  std::vector<TString> beam_files;
  Bool_t beam_eres = kFALSE;
  for (Int_t i = 0; i < Int_t(specs.size()); i++) {
    TString base = RemixSim::TagWithoutStrip(specs[i].tag);
    base.ReplaceAll("_eres", "");
    TString file = RemixSim::SimRootPath(specs[i]);
    Int_t strip = RemixSim::ReactionStripOf(specs[i].tag);
    const Bool_t eres = RemixSim::IsEresTag(specs[i].tag);
    if (base == "beam") {
      if (beam_eres && !eres)
        continue;
      if (eres && !beam_eres)
        beam_files.clear();
      beam_files.push_back(file);
      beam_eres = beam_eres || eres;
    } else if (strip >= kReacMin && strip <= kReacMax) {
      std::map<Int_t, TString> &files = base == "aa" ? aa_file : an_file;
      std::map<Int_t, Bool_t> &have_eres = base == "aa" ? aa_eres : an_eres;
      if (base != "aa" && base != "an")
        continue;
      if (files.count(strip) && (have_eres[strip] || !eres))
        continue;
      files[strip] = file;
      have_eres[strip] = eres;
    }
  }
 
  Double_t gain[18];
  if (!SimBeamGains(gain))
    for (Int_t s = 0; s < 18; s++)
      gain[s] = 1.0;
 
  std::vector<TGraph *> beam_traces;
  for (Int_t i = 0; i < Int_t(beam_files.size()) &&
                    Int_t(beam_traces.size()) < kTracesPerRegion;
       i++) {
    std::vector<TGraph *> t = SimPopTraces(
        beam_files[i], gain, kTracesPerRegion - Int_t(beam_traces.size()));
    for (Int_t k = 0; k < Int_t(t.size()); k++)
      beam_traces.push_back(t[k]);
  }
 
  for (Int_t r = kReacMin; r <= kReacMax; r++) {
    std::vector<TGraph *> aa_traces, an_traces;
    if (aa_file.find(r) != aa_file.end())
      aa_traces = SimPopTraces(aa_file[r], gain, kTracesPerRegion);
    if (an_file.find(r) != an_file.end())
      an_traces = SimPopTraces(an_file[r], gain, kTracesPerRegion);
    if (aa_traces.empty() && an_traces.empty())
      continue;
    DrawRegionTraces(Form("sim_region_traces_reac%d", r), "sim_scatter",
                     beam_traces, aa_traces, an_traces, 0.6, 1.6,
                     "#DeltaE [a.u.]");
    for (Int_t i = 0; i < Int_t(aa_traces.size()); i++)
      delete aa_traces[i];
    for (Int_t i = 0; i < Int_t(an_traces.size()); i++)
      delete an_traces[i];
  }
  for (Int_t i = 0; i < Int_t(beam_traces.size()); i++)
    delete beam_traces[i];
}
 
// Fingerprint of the sim inputs + window geometry: each eres file's size+mtime
// (cheap, no open) plus the reaction-strip range and x window. Regenerating the
// sim (new mtimes) or changing the windows invalidates the cached overlay.
TString StripSumScatter::SimFingerprint(
    const std::vector<RemixSim::SimFileSpec> &specs) {
  // v2: sim is per-strip normalized via SimBeamGains(); the gain source (the
  // eres beam file) is stamped by the per-spec loop below (it covers every
  // eres file, including beam_eres).
  // v3: normalization hardcoded to 1 a.u. (NORM_MUSIC_MEV removed); bump so
  // overlays cached at other norms rebuild.
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
 
  // v4: one population per class and strip (the eres twin preferred), so
  // caches that held every population twice are rebuilt.
  TString s = Form("v4 reac[%d,%d] x[%d,%d]", kReacMin, kReacMax, kXLo, kXHi);
  for (Int_t i = 0; i < Int_t(specs.size()); i++) {
    TString f = RemixSim::SimRootPath(specs[i]);
    Long_t id = 0, flags = 0, mtime = 0;
    Long64_t size = -1;
    if (gSystem->GetPathInfo(f, &id, &size, &flags, &mtime) != 0) {
      size = -1;
      mtime = 0;
    }
    s += Form(" %s:%lld:%ld", specs[i].tag.Data(), size, mtime);
  }
  return s;
}
 
// Reload cached sim scatter graphs (grouped by reaction strip; each graph's
// title holds its population label) if the fingerprint matches. Caller owns the
// returned graphs.
Bool_t StripSumScatter::LoadSimCache(
    const TString &fp, std::map<Int_t, std::vector<TGraph *>> &by_strip) {
  TString full = IO::GetRootFilesBaseDir() + TString("/") +
                 "StripSumScatter_simcache.root";
  if (gSystem->AccessPathName(full))
    return kFALSE;
  TFile *f = IO::OpenForReading("StripSumScatter_simcache.root");
  if (!f || f->IsZombie()) {
    if (f)
      delete f;
    return kFALSE;
  }
  TNamed *cfp = static_cast<TNamed *>(f->Get("sim_fingerprint"));
  if (!cfp || fp != cfp->GetTitle()) {
    f->Close();
    delete f;
    return kFALSE;
  }
  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = static_cast<TKey *>(next()))) {
    TString name = key->GetName();
    if (!name.BeginsWith("simg_r"))
      continue;
    TString rest = name(6, name.Length() - 6); // after "simg_r": <strip>_p<idx>
    Int_t us = rest.Index("_p");
    if (us < 0)
      continue;
    Int_t r = TString(rest(0, us)).Atoi();
    TGraph *g = static_cast<TGraph *>(f->Get(name));
    if (!g)
      continue;
    by_strip[r].push_back(static_cast<TGraph *>(g->Clone()));
  }
  f->Close();
  delete f;
  return kTRUE;
}
 
void StripSumScatter::WriteSimCache(
    const TString &fp, const std::map<Int_t, std::vector<TGraph *>> &by_strip) {
  TFile *out = IO::OpenForWriting("StripSumScatter_simcache.root", "RECREATE");
  if (!out || out->IsZombie()) {
    if (out)
      delete out;
    return;
  }
  out->cd();
  TNamed cfp("sim_fingerprint", fp.Data());
  cfp.Write();
  std::map<Int_t, std::vector<TGraph *>>::const_iterator it;
  for (it = by_strip.begin(); it != by_strip.end(); ++it)
    for (Int_t i = 0; i < Int_t(it->second.size()); i++)
      it->second[i]->Write(Form("simg_r%d_p%d", it->first, Int_t(i)));
  out->Close();
  delete out;
}
 
// Sim-only comparison plots: one per reaction strip, each sim population a
// coloured+labelled point cloud on the same axes as that strip's data scatter,
// for side-by-side comparison with the data PID scatters. Beam (no reaction
// strip) overlays on every strip. The scatter graphs are fingerprint-cached
// (sim file sizes/mtimes + window geometry), so re-runs reload them instead of
// rescanning the sim files.
void StripSumScatter::SimOverlay() {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
 
  std::vector<RemixSim::SimFileSpec> specs = RemixSim::BuildFileSpecs();
  if (specs.empty()) {
    std::cerr
        << "strip-sum-scatter: no sim control files; skipping sim overlay."
        << std::endl;
    return;
  }
  TString fp = SimFingerprint(specs);
 
  std::map<Int_t, std::vector<TGraph *>>
      by_strip; // strip -> graphs (title=label)
  Bool_t loaded = LoadSimCache(fp, by_strip);
  if (!loaded) {
    // One population per class and strip. Every population has two control
    // files, the plain one and its _eres twin with the measured widths, and
    // both carry the same label; taking each file as a population drew every
    // class twice. The eres file wins where it exists.
    std::map<std::pair<TString, Int_t>, std::pair<SimPop, Bool_t>> chosen;
    for (Int_t i = 0; i < Int_t(specs.size()); i++) {
      TString base = RemixSim::TagWithoutStrip(specs[i].tag);
      base.ReplaceAll("_eres", "");
      const Int_t strip = RemixSim::ReactionStripOf(specs[i].tag);
      const Bool_t eres = RemixSim::IsEresTag(specs[i].tag);
      std::pair<TString, Int_t> key(base, strip);
      if (chosen.count(key) && (chosen[key].second || !eres))
        continue;
      SimPop p;
      p.file = RemixSim::SimRootPath(specs[i]);
      p.label = PrettyLabel(specs[i].tag);
      chosen[key] = std::make_pair(p, eres);
    }
    std::map<Int_t, std::vector<SimPop>> reacted;
    std::vector<SimPop> refs;
    for (std::map<std::pair<TString, Int_t>,
                  std::pair<SimPop, Bool_t>>::const_iterator it =
             chosen.begin();
         it != chosen.end(); ++it) {
      if (!it->second.second)
        std::cout << "No eres sim file for " << it->first.first
                  << (it->first.second >= 0 ? Form("_s%d", it->first.second)
                                            : "")
                  << ", using standard" << std::endl;
      if (it->first.second < 0)
        refs.push_back(it->second.first);
      else
        reacted[it->first.second].push_back(it->second.first);
    }
    Double_t gain[18];
    if (!SimBeamGains(gain))
      for (Int_t s = 0; s < 18; s++)
        gain[s] = 1.0;
    const Long64_t kSimMaxPoints = 25000;
    for (Int_t r = kReacMin; r <= kReacMax; r++) {
      std::vector<SimPop> group = reacted[r];
      for (Int_t i = 0; i < Int_t(refs.size()); i++)
        group.push_back(refs[i]);
      for (Int_t i = 0; i < Int_t(group.size()); i++) {
        TGraph *g = SimPopScatter(group[i].file, r, gain, kSimMaxPoints);
        if (!g || g->GetN() == 0) {
          if (g)
            delete g;
          continue;
        }
        g->SetTitle(group[i].label);
        by_strip[r].push_back(g);
      }
    }
    Int_t n_graphs = 0;
    std::map<Int_t, std::vector<TGraph *>>::const_iterator cit;
    for (cit = by_strip.begin(); cit != by_strip.end(); ++cit)
      n_graphs += Int_t(cit->second.size());
    if (n_graphs == 0) {
      std::cerr << "strip-sum-scatter: no sim data found (regenerate "
                   "sim_root_files); skipping sim overlay."
                << std::endl;
      return;
    }
    WriteSimCache(fp, by_strip);
    std::cout << "strip-sum-scatter: built + cached sim overlay (" << n_graphs
              << " population graphs)." << std::endl;
  } else {
    std::cout
        << "strip-sum-scatter: loaded cached sim overlay (fingerprint match)."
        << std::endl;
  }
 
  std::map<Int_t, std::vector<TGraph *>>::iterator it;
  for (it = by_strip.begin(); it != by_strip.end(); ++it) {
    Int_t r = it->first;
    std::map<Int_t, TH2F *>::const_iterator sit = m_scatter.find(r);
    if (sit == m_scatter.end())
      continue;
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    TH2F *ref = sit->second;
    TH2F *frame =
        new TH2F(Form("sim_frame_r%d", r), "", 10, ref->GetXaxis()->GetXmin(),
                 ref->GetXaxis()->GetXmax(), 10, ref->GetYaxis()->GetXmin(),
                 ref->GetYaxis()->GetXmax());
    frame->SetStats(0);
    frame->GetXaxis()->SetTitle(ref->GetXaxis()->GetTitle());
    frame->GetYaxis()->SetTitle(ref->GetYaxis()->GetTitle());
    TCanvas *c = PlottingUtils::GetConfiguredCanvas(kFALSE);
    c->SetLeftMargin(0.18);
    frame->Draw();
    // Match the experimental region-traces legend placement (top-right).
    TLegend *leg = PlottingUtils::AddLegend(0.725, 0.875, 0.70, 0.86);
    for (Int_t i = 0; i < Int_t(it->second.size()); i++) {
      TGraph *g = it->second[i];
      // Match the experimental region-trace colours (DrawRegionTraces): beam
      // grey, (a,a') azure, (a,n) red -- keyed off the population label.
      TString lab = g->GetTitle();
      Int_t color = kBlack;
      if (lab == "Beam")
        color = kGray + 2;
      else if (lab == "(#alpha,#alpha')")
        color = kAzure + 2;
      else if (lab == "(#alpha,n)")
        color = kRed + 1;
      g->SetMarkerStyle(20);
      g->SetMarkerSize(0.3);
      g->SetMarkerColorAlpha(color, 0.35);
      g->SetLineColor(color);
      g->Draw("P SAME");
      leg->AddEntry(g, g->GetTitle(), "p");
    }
    leg->Draw();
    PlottingUtils::SaveFigure(c,
                              Form("sim_normsumE_reac%d_s%d_%d_vs_s%d_%d", r,
                                   YLoOf(r), YHiOf(r), kXLo, kXHi),
                              "sim_scatter", PlotSaveOptions::kLINEAR);
    delete leg;
    delete c;
    delete frame;
  }
 
  std::map<Int_t, std::vector<TGraph *>>::iterator dit;
  for (dit = by_strip.begin(); dit != by_strip.end(); ++dit)
    for (Int_t i = 0; i < Int_t(dit->second.size()); i++)
      delete dit->second[i];
}
 
Bool_t StripSumScatter::TryLoadCache(const TString &cacheName,
                                     const TString &fingerprint) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
 
  TString cache_full = IO::GetRootFilesBaseDir() + TString("/") + cacheName;
  if (gSystem->AccessPathName(cache_full)) {
    std::cout << "strip-sum-scatter: no cache file found; will rebuild."
              << std::endl;
    return kFALSE;
  }
 
  TFile *cf = IO::OpenForReading(cacheName);
  if (!cf || cf->IsZombie()) {
    if (cf)
      delete cf;
    std::cout << "strip-sum-scatter: cache file unreadable; rebuilding."
              << std::endl;
    return kFALSE;
  }
 
  TNamed *fp = static_cast<TNamed *>(cf->Get("fingerprint"));
  const Bool_t exact = fp && fingerprint == fp->GetTitle();
  // Same tagging, different plane: the reservoir has every tagged event, so
  // the scatters are rebuilt from it below and the cache rewritten.
  const Bool_t reproject = !exact && fp && !TagPart(fingerprint).IsNull() &&
                           TagPart(fingerprint) == TagPart(fp->GetTitle());
  if (!exact && !reproject) {
    std::cout << "strip-sum-scatter: cache present but stale; rebuilding."
              << std::endl;
    std::cout << "  cached: " << (fp ? fp->GetTitle() : "(none)") << std::endl;
    std::cout << "  wanted: " << fingerprint << std::endl;
    cf->Close();
    delete cf;
    return kFALSE;
  }
 
  Bool_t ok = kTRUE;
  for (Int_t reac = kReacMin; reac <= kReacMax && ok && exact; reac++) {
    TH2F *h = static_cast<TH2F *>(cf->Get(Form("scatter_r%d", reac)));
    if (!h) {
      ok = kFALSE;
      break;
    }
    TH2F *hc = static_cast<TH2F *>(h->Clone());
    hc->SetDirectory(nullptr);
    m_scatter[reac] = hc;
  }
 
  // The normalization counts ride along, so a re-projected cache keeps them.
  if (TParameter<Long64_t> *p =
          dynamic_cast<TParameter<Long64_t> *>(cf->Get("n_seen")))
    m_nSeen = p->GetVal();
  if (TParameter<Long64_t> *p =
          dynamic_cast<TParameter<Long64_t> *>(cf->Get("n_normed")))
    m_nNormed = p->GetVal();
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
    if (TParameter<Long64_t> *p = dynamic_cast<TParameter<Long64_t> *>(
            cf->Get(Form("n_normed_r%d", reac))))
      m_normedAt[ReacIndex(reac)] = p->GetVal();
    if (TParameter<Long64_t> *p = dynamic_cast<TParameter<Long64_t> *>(
            cf->Get(Form("n_tagged_r%d", reac))))
      m_tagged[ReacIndex(reac)] = p->GetVal();
  }
 
  TTree *tt = static_cast<TTree *>(cf->Get("traces"));
  if (reproject && !tt)
    ok = kFALSE;
  if (ok && tt) {
    TraceEvt e;
    tt->SetBranchAddress("total", e.total);
    tt->SetBranchAddress("total_adc", e.total_adc);
    tt->SetBranchAddress("reac_mask", &e.reac_mask);
    tt->SetBranchAddress("beam_flat", &e.beam_flat);
    tt->SetBranchAddress("both_mult", &e.both_mult);
    e.seed_ts = 0;
    if (tt->GetBranch("seed_ts"))
      tt->SetBranchAddress("seed_ts", &e.seed_ts);
    Long64_t nt = tt->GetEntries();
    m_reservoir.reserve(nt);
    for (Long64_t j = 0; j < nt; j++) {
      tt->GetEntry(j);
      m_reservoir.push_back(e);
    }
  }
 
  cf->Close();
  delete cf;
 
  if (!ok) {
    std::cout << "strip-sum-scatter: cache partially corrupt; rebuilding."
              << std::endl;
    m_reservoir.clear();
    return kFALSE;
  }
  if (reproject) {
    std::cout << "strip-sum-scatter: cache tagging matches but the plane "
                 "changed; re-projecting "
              << m_reservoir.size() << " reservoir events." << std::endl;
    AllocateScatters();
    ReprojectFromReservoir();
    WriteCache(cacheName, fingerprint);
  } else {
    std::cout << "strip-sum-scatter: loaded cached scatters + "
              << m_reservoir.size() << " reservoir events (fingerprint match)."
              << std::endl;
  }
  return kTRUE;
}
 
void StripSumScatter::WriteCache(const TString &cacheName,
                                 const TString &fingerprint) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
 
  TFile *out = IO::OpenForWriting(cacheName, "RECREATE");
  if (!out || out->IsZombie()) {
    if (out)
      delete out;
    return;
  }
  out->cd();
  TNamed fp("fingerprint", fingerprint.Data());
  fp.Write();
  // Normalization counts travel with the scatters they describe, so a cross
  // section never has to re-derive them from a different pass over the data.
  TParameter<Long64_t>("n_seen", m_nSeen).Write();
  TParameter<Long64_t>("n_normed", m_nNormed).Write();
  // The noise the jump gate and the upstream tolerance were resolved against.
  for (Int_t s = 1; s < 18; s++)
    TParameter<Double_t>(Form("jump_sigma_s%d", s), s_jumpSigma[s]).Write();
  for (Int_t s = 0; s < 18; s++)
    TParameter<Double_t>(Form("strip_sigma_s%d", s), s_stripSigma[s]).Write();
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
    TParameter<Long64_t>(Form("n_normed_r%d", reac),
                         m_normedAt[ReacIndex(reac)])
        .Write();
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
    TParameter<Long64_t>(Form("n_tagged_r%d", reac), m_tagged[ReacIndex(reac)])
        .Write();
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
    m_scatter[reac]->Write(Form("scatter_r%d", reac));
 
  TTree *tt = new TTree("traces", "strip-sum trace reservoir");
  TraceEvt e;
  tt->Branch("total", e.total, "total[18]/F");
  tt->Branch("total_adc", e.total_adc, "total_adc[18]/F");
  tt->Branch("reac_mask", &e.reac_mask, "reac_mask/i");
  tt->Branch("beam_flat", &e.beam_flat, "beam_flat/O");
  tt->Branch("both_mult", &e.both_mult, "both_mult/I");
  tt->Branch("seed_ts", &e.seed_ts, "seed_ts/l");
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    e = m_reservoir[k];
    tt->Fill();
  }
  tt->Write();
  out->Close();
  delete out;
  std::cout << "strip-sum-scatter: wrote cache " << cacheName << std::endl;
}
 
// One run's beam ellipses and series gates. Split out so runs can be fitted in
// parallel: each call touches only its own chain and returns its own result,
// with no shared state to guard.
SingleRunFitResult
StripSumScatter::FitRunGates(Int_t run, TChain *chain,
                             const std::vector<GateSpec> &activeGates) {
  SingleRunFitResult res;
  if (!chain || chain->GetEntries() == 0)
    return res;
 
  // --- Beam classification ellipses (no gating) ---
  BeamEllipses be;
  be.ok = kFALSE;
  {
    std::vector<GateSpec> emptyPrior;
    std::vector<BeamFit2D> emptyGates;
    const TString tag = Form("run%d", run);
    const TString subdir = Form("strip_sum_scatter/run%d", run);
    Int_t ent_sx = 0, ent_sy = 1;
    const Char_t *ent_tag = "s0/s1";
    if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.PURE_BEAM_GATE ==
        StripSumScatterConfig::PURE_BEAM_GATE_S1_S2) {
      ent_sx = 1;
      ent_sy = 2;
      ent_tag = "s1/s2";
    }
    BeamFit2D ent_ell = FindBeamGate(chain, ent_sx, ent_sy, emptyPrior,
                                     emptyGates, tag, subdir);
    if (ent_sx == 0)
      be.s0_s1 = ent_ell;
    else
      be.s1_s2 = ent_ell;
    if (ent_ell.ok) {
      std::lock_guard<std::mutex> lk(g_log_mutex);
      std::cout << "  run " << run << " beam ellipse " << ent_tag << ": mu=("
                << ent_ell.mu_x << "," << ent_ell.mu_y << ")" << std::endl;
    } else {
      std::lock_guard<std::mutex> lk(g_log_mutex);
      std::cerr << "  run " << run << " beam ellipse " << ent_tag
                << " failed; skipping run" << std::endl;
      return res;
    }
    if (Constants::cfg.IGNORE_STRIP_17) {
      be.use_s15_s16 = kTRUE;
      be.s15_s16 =
          FindBeamGate(chain, 15, 16, emptyPrior, emptyGates, tag, subdir);
      if (be.s15_s16.ok) {
        std::lock_guard<std::mutex> lk(g_log_mutex);
        std::cout << "  run " << run << " beam ellipse s15/s16: mu=("
                  << be.s15_s16.mu_x << "," << be.s15_s16.mu_y << ")"
                  << std::endl;
      } else {
        std::lock_guard<std::mutex> lk(g_log_mutex);
        std::cerr << "  run " << run
                  << " beam ellipse s15/s16 failed; skipping run" << std::endl;
        return res;
      }
    } else {
      be.use_s15_s16 = kFALSE;
      be.s16_s17 =
          FindBeamGate(chain, 16, 17, emptyPrior, emptyGates, tag, subdir);
      if (be.s16_s17.ok) {
        std::lock_guard<std::mutex> lk(g_log_mutex);
        std::cout << "  run " << run << " beam ellipse s16/s17: mu=("
                  << be.s16_s17.mu_x << "," << be.s16_s17.mu_y << ")"
                  << std::endl;
      } else {
        std::lock_guard<std::mutex> lk(g_log_mutex);
        std::cerr << "  run " << run
                  << " beam ellipse s16/s17 failed; skipping run" << std::endl;
        return res;
      }
    }
    be.ok = kTRUE;
  }
  res.pure_beam = be;
 
  // --- Scatter filter gates (series gating) ---
  std::vector<BeamFit2D> runGates;
  std::vector<GateSpec> priorSpecs;
  Bool_t allOk = kTRUE;
  for (Int_t gi = 0; gi < Int_t(activeGates.size()); gi++) {
    BeamFit2D g = FindBeamGate(chain, activeGates[gi].sx, activeGates[gi].sy,
                               priorSpecs, runGates, Form("run%d", run),
                               Form("strip_sum_scatter/run%d", run));
    if (g.ok) {
      std::lock_guard<std::mutex> lk(g_log_mutex);
      std::cout << "  run " << run << " beam gate s" << activeGates[gi].sx
                << "/s" << activeGates[gi].sy << ": mu=(" << g.mu_x << ","
                << g.mu_y << ")" << std::endl;
    } else {
      std::lock_guard<std::mutex> lk(g_log_mutex);
      std::cerr << "  run " << run << " beam gate s" << activeGates[gi].sx
                << "/s" << activeGates[gi].sy << " failed; skipping run"
                << std::endl;
      allOk = kFALSE;
    }
    runGates.push_back(g);
    priorSpecs.push_back(activeGates[gi]);
  }
  res.series_gates = runGates;
  res.ok = allOk;
  return res;
}
 
// One run's scatter fill. Each call builds PRIVATE scatter histograms and its
// own reservoir slice, so the runs never touch shared state; the caller merges
// them in run order, which makes the threaded result identical to sequential.
SingleRunFillResult StripSumScatter::FillRunScatters(
    Int_t run, TChain *chain, const std::vector<GateSpec> &activeGates,
    const std::vector<BeamFit2D> &runGates, const BeamEllipses &runBeam) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kXBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.XBINS;
  const Int_t kYBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.YBINS;
  const Int_t kBeamReservoirCap =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRACES_PER_CLASS * 10;
  const Int_t nReacStrips = kReacMax - kReacMin + 1;
 
  SingleRunFillResult res;
  res.scatters.assign(nReacStrips, nullptr);
  if (!chain)
    return res;
  // Private, directory-less clones over the same fixed build range as the
  // merged ones.
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
    TH2F *h =
        new TH2F(Form("scatter_r%d_run%d", reac, run), "", kXBins,
                 ScatterBuildRange::kXMin, ScatterBuildRange::kXMax, kYBins,
                 ScatterBuildRange::kYMin, ScatterBuildRange::kYMax);
    h->SetDirectory(nullptr);
    res.scatters[ReacIndex(reac)] = h;
  }
  Long64_t totalGated = 0, totalSeen = 0, totalNormed = 0;
  res.tagged.assign(nReacStrips, 0);
  res.normed_at.assign(nReacStrips, 0);
  Int_t nBeamKept = 0;
  EnergyView ev;
  ev.Attach(chain);
  EnableEventBranches(chain);
  ULong64_t seed_ts_in = 0;
  if (chain->GetBranch("SeedTs"))
    chain->SetBranchAddress("SeedTs", &seed_ts_in);
  Long64_t n = chain->GetEntries();
  Int_t nReac = kReacMax - kReacMin + 1;
  {
    std::lock_guard<std::mutex> lk(g_log_mutex);
    std::cout << "Run " << run << ": filling " << nReac
              << " reaction-strip scatters over " << n << " events..."
              << std::endl;
  }
 
  for (Long64_t j = 0; j < n; j++) {
    chain->GetEntry(j);
    ev.Decode();
    totalSeen++;
 
    Bool_t passesAll = kTRUE;
    for (Int_t gi = 0; gi < Int_t(activeGates.size()); gi++)
      if (!PassesGate(runGates[gi], ev, activeGates[gi].sx,
                      activeGates[gi].sy)) {
        passesAll = kFALSE;
        break;
      }
    if (!passesAll)
      continue;
    if (IsPileup(ev)) // reject overlapping-beam pileup
      continue;
    if (IsNoise(ev))
      continue;
    if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REJECT_OFFBEAM && IsOffbeam(ev))
      continue;
    if (IsParityAsymmetric(ev))
      continue;
    // Both-ends multiplicity: counted on raw ADC so it sees the short end
    // even when IGNORE_SHORT_STRIPS zeroes it in the decode.
    if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BOTH_MULT_MAX >= 0) {
      const Int_t hi = TMath::Min(
          16, Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BOTH_MULT_COUNT_TO);
      Int_t nboth = 0;
      for (Int_t s = 1; s <= hi; s++)
        if (ev.left_0_17_adc[s] > 0.0 && ev.rightdE_adc[s] > 0.0)
          nboth++;
      if (nboth > Constants::cfg.STRIP_SUM_SCATTER_CONFIG.BOTH_MULT_MAX)
        continue;
    }
    // Last point at which nothing about a reaction has been asked. Counting
    // here, rather than at `seen`, is what makes the ratio to a tag count a
    // cross section: both sides carry the same gate and quality efficiencies.
    totalNormed++;
    // Per-strip denominator. A beam particle counts toward strip `reac` only
    // if it met the conditions a reaction there would also have had to meet,
    // so those efficiencies cancel in the ratio rather than tilting the
    // excitation function. Everything after this in PassesReaction is about
    // the reaction itself and belongs only to the numerator.
    if (AllStripsFired(ev)) {
      for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
        if (BeamUpstreamOf(ev, reac))
          res.normed_at[ReacIndex(reac)]++;
    }
 
    UInt_t mask = 0;
    for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
      if (!PassesReaction(ev, reac))
        continue;
      mask |= (1u << ReacIndex(reac));
      res.tagged[ReacIndex(reac)]++;
      Double_t x = 0.0, y = 0.0;
      PlaneXY(ev.total, reac, x, y);
      res.scatters[ReacIndex(reac)]->Fill(x, y);
    }
 
    // Keep reaction-passing events for traces; cap pure-beam events (only
    // ~TRACES_PER_CLASS are ever drawn). The two are mutually exclusive -- a
    // pure-beam event has no reaction jump. This cap is only a per-task
    // memory bound; FillScatters re-applies it across all tasks when it
    // merges, which is what actually fixes the kept beam population.
    Bool_t beam = (mask == 0) && IsPureBeam(ev, runBeam);
    if (mask == 0 && !(beam && nBeamKept < kBeamReservoirCap))
      continue;
    if (beam)
      nBeamKept++;
 
    TraceEvt e;
    for (Int_t s = 0; s < 18; s++) {
      e.total[s] = Float_t(ev.total[s]);
      e.total_adc[s] =
          Float_t(ev.left_0_17_adc[s]) + Float_t(ev.rightdE_adc[s]);
      e.long_au[s] = Float_t(ev.total[s]);
      e.short_au[s] = 0.0f;
    }
    for (Int_t s = 1; s <= 16; s++) {
      Double_t lv = Double_t(ev.gain_left[s]) * Double_t(ev.left_0_17_adc[s]) *
                    Double_t(ev.strip_factor[s]);
      Double_t rv = Double_t(ev.gain_right[s]) * Double_t(ev.rightdE_adc[s]) *
                    Double_t(ev.strip_factor[s]);
      Bool_t l_is_long = ((s % 2) != 0);
      e.long_au[s] = Float_t(l_is_long ? lv : rv);
      e.short_au[s] = Float_t(l_is_long ? rv : lv);
    }
    // Mirror IGNORE_SHORT_STRIPS: the normed total keeps only the long side
    // of a split strip, so the raw trace must drop the same side to stay
    // comparable.
    if (Constants::cfg.IGNORE_SHORT_STRIPS)
      for (Int_t s = 1; s <= 16; s++)
        e.total_adc[s] = ((s % 2) != 0) ? Float_t(ev.left_0_17_adc[s])
                                        : Float_t(ev.rightdE_adc[s]);
    // Both-channel multiplicity: split strips (1-16) where both ends
    // FIRED. Read off the RAW ADC, not the calibrated ends -- the
    // short-end gains are 0 (uncalibrated, no sim anchor), so the
    // calibrated short ends are always zero; the raw ADC still carries
    // whether the channel fired.
    Int_t both = 0;
    for (Int_t s = 1; s <= 16; s++)
      if (ev.left_0_17_adc[s] > 0.0 && ev.rightdE_adc[s] > 0.0)
        both++;
    e.both_mult = both;
    e.seed_ts = seed_ts_in;
    e.reac_mask = mask;
    e.beam_flat = beam;
    res.reservoir.push_back(e);
    if (mask != 0)
      totalGated++;
  }
  res.gated = totalGated;
  res.seen = totalSeen;
  res.normed = totalNormed;
  return res;
}
 
// Run `n` indexed tasks on `workers` threads, pulling from a shared queue.
static void RunIndexedParallel(Int_t n, Int_t workers,
                               const std::function<void(Int_t)> &task) {
  std::queue<Int_t> work;
  for (Int_t i = 0; i < n; i++)
    work.push(i);
  std::mutex work_mutex;
  std::vector<std::thread> pool;
  for (Int_t w = 0; w < workers; w++) {
    pool.emplace_back([&]() {
      while (true) {
        Int_t i;
        {
          std::lock_guard<std::mutex> lk(work_mutex);
          if (work.empty())
            return;
          i = work.front();
          work.pop();
        }
        task(i);
      }
    });
  }
  for (Int_t w = 0; w < Int_t(pool.size()); w++)
    pool[w].join();
}
 
void StripSumScatter::AllocateScatters() {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kXBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.XBINS;
  const Int_t kYBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.YBINS;
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
    TH2F *h =
        new TH2F(Form("scatter_r%d", reac),
                 Form(";norm. #DeltaE strips %d#rightarrow%d [a.u.];norm. "
                      "#DeltaE strips %d#rightarrow%d [a.u.]",
                      kXLo, kXHi, YLoOf(reac), YHiOf(reac)),
                 kXBins, ScatterBuildRange::kXMin, ScatterBuildRange::kXMax,
                 kYBins, ScatterBuildRange::kYMin, ScatterBuildRange::kYMax);
    h->SetDirectory(nullptr);
    h->SetStats(0);
    m_scatter[reac] = h;
  }
}
 
void StripSumScatter::ReprojectFromReservoir() {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  Long64_t nFilled = 0;
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    const TraceEvt &e = m_reservoir[k];
    if (e.reac_mask == 0)
      continue;
    Double_t total[18];
    for (Int_t s = 0; s < 18; s++)
      total[s] = e.total[s];
    for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
      if (!(e.reac_mask & (1u << ReacIndex(reac))))
        continue;
      Double_t x = 0.0, y = 0.0;
      PlaneXY(total, reac, x, y);
      m_scatter[reac]->Fill(x, y);
      nFilled++;
    }
  }
  std::cout << "strip-sum-scatter: re-projected " << nFilled
            << " tagged entries." << std::endl;
}
 
void StripSumScatter::FillScatters(const std::vector<Int_t> &runOrder,
                                   std::map<Int_t, TChain *> &chains) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.XBINS;
  const Int_t kYBins = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.YBINS;
 
  AllocateScatters();
 
  std::vector<GateSpec> activeGates = ActiveGates();
 
  // Pre-index the chains so worker threads only touch their own run.
  const Int_t nRuns = Int_t(runOrder.size());
  std::vector<TChain *> chainVec(nRuns);
  for (Int_t i = 0; i < nRuns; i++)
    chainVec[i] = chains[runOrder[i]];
 
  Int_t n_workers =
      TMath::Min(Int_t(std::thread::hardware_concurrency()), nRuns);
  n_workers = TMath::Min(
      n_workers, Constants::cfg.STRIP_SUM_SCATTER_CONFIG.MAX_STRIP_SUM_WORKERS);
  if (n_workers < 1)
    n_workers = 1;
  std::cout << "strip-sum-scatter: " << nRuns << " runs on " << n_workers
            << " workers" << std::endl;
 
  // Phase 1: beam ellipses + series gates, one task per run. The gates within
  // a run stay sequential -- each only sees events passing the prior ones.
  std::vector<SingleRunFitResult> fits(nRuns);
  RunIndexedParallel(nRuns, n_workers, [&](Int_t i) {
    fits[i] = FitRunGates(runOrder[i], chainVec[i], activeGates);
  });
 
  // Phase 2: event filling. Work is split per EVENTS FILE, not per run: a
  // CoMPASS run is hundreds of subfiles collapsed into one chain by
  // GroupEventsByRun, so per-run tasks would leave 87Rb at two-way parallelism
  // while a SOLARIS dataset with one file per run already gets full width.
  // Gates stay per-run above (they need whole-run statistics); the fill does
  // not, so each file is its own task carrying its run's gates.
  struct FillTask {
    Int_t run_idx;
    TString path;
  };
  std::vector<FillTask> tasks;
  {
    std::map<Int_t, Int_t> idx_of_run;
    for (Int_t i = 0; i < nRuns; i++)
      idx_of_run[runOrder[i]] = i;
    std::vector<FileSpec> specs = FileSet::BuildProcessedFileSpecs();
    for (Int_t k = 0; k < Int_t(specs.size()); k++) {
      std::map<Int_t, Int_t>::const_iterator it = idx_of_run.find(specs[k].run);
      if (it == idx_of_run.end() || !fits[it->second].ok)
        continue;
      TString full = IO::GetRootFilesBaseDir() + "/" +
                     FileSet::EventsName(specs[k]) + ".root";
      if (gSystem->AccessPathName(full))
        continue;
      FillTask t;
      t.run_idx = it->second;
      t.path = full;
      tasks.push_back(t);
    }
  }
  Int_t nTasks = Int_t(tasks.size());
  Int_t fill_workers =
      TMath::Min(Int_t(std::thread::hardware_concurrency()), nTasks);
  fill_workers =
      TMath::Min(fill_workers,
                 Constants::cfg.STRIP_SUM_SCATTER_CONFIG.MAX_STRIP_SUM_WORKERS);
  if (fill_workers < 1)
    fill_workers = 1;
  std::cout << "strip-sum-scatter: filling " << nTasks << " files on "
            << fill_workers << " workers" << std::endl;
 
  // Every task's result is merged into the totals in TASK ORDER as soon as
  // every task before it has finished, and freed on the spot. Holding all
  // results until the end and merging afterwards kept fourteen scatters of
  // XBINS x YBINS bins and a reservoir slice per file alive at once, then
  // copied every slice into the total while the slices still existed: for a
  // SOLARIS dataset of seventy files that is tens of GB at the moment of the
  // merge, and the process was killed there. The order is kept so the pure-
  // beam budget below is spent the same way whatever order the workers
  // finish in, which makes the threaded result identical to sequential.
  std::vector<SingleRunFillResult> fills(nTasks);
  std::vector<Bool_t> filled(nTasks, kFALSE);
  Int_t next_merge = 0;
  std::mutex merge_mutex;
  Long64_t totalGated = 0, totalSeen = 0;
  const Int_t kBeamReservoirCap =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRACES_PER_CLASS * 10;
  Int_t nBeamKept = 0;
  auto merge = [&](Int_t t) {
    for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
      Int_t ri = ReacIndex(reac);
      if (ri < Int_t(fills[t].scatters.size()) && fills[t].scatters[ri])
        m_scatter[reac]->Add(fills[t].scatters[ri]);
    }
    // Pure-beam events are capped GLOBALLY, not per task. Each task applies
    // the same cap to its own slice, so without this the kept beam population
    // would scale with the number of events files -- hundreds of them for a
    // CoMPASS run. Reaction-tagged events are never dropped.
    for (Int_t k = 0; k < Int_t(fills[t].reservoir.size()); k++) {
      const TraceEvt &e = fills[t].reservoir[k];
      if (e.beam_flat) {
        if (nBeamKept >= kBeamReservoirCap)
          continue;
        nBeamKept++;
      }
      m_reservoir.push_back(e);
    }
    totalGated += fills[t].gated;
    totalSeen += fills[t].seen;
    m_nNormed += fills[t].normed;
    for (Int_t k = 0; k < Int_t(fills[t].normed_at.size()); k++)
      m_normedAt[k] += fills[t].normed_at[k];
    for (Int_t k = 0; k < Int_t(fills[t].tagged.size()); k++)
      m_tagged[k] += fills[t].tagged[k];
    for (Int_t k = 0; k < Int_t(fills[t].scatters.size()); k++)
      delete fills[t].scatters[k];
    fills[t].scatters.clear();
    std::vector<TraceEvt>().swap(fills[t].reservoir);
  };
  RunIndexedParallel(nTasks, fill_workers, [&](Int_t t) {
    Int_t i = tasks[t].run_idx;
    {
      TChain ch("events");
      ch.Add(tasks[t].path);
      if (ch.GetEntries() > 0)
        fills[t] = FillRunScatters(runOrder[i], &ch, activeGates,
                                   fits[i].series_gates, fits[i].pure_beam);
    }
    std::lock_guard<std::mutex> lk(merge_mutex);
    filled[t] = kTRUE;
    while (next_merge < nTasks && filled[next_merge]) {
      merge(next_merge);
      next_merge++;
    }
  });
  m_nSeen = totalSeen;
  std::cout << "strip-sum-scatter: " << totalGated << " reaction-tagged of "
            << totalSeen << " events (" << m_nNormed
            << " past every pre-tag cut); reservoir " << m_reservoir.size()
            << std::endl;
}
 
void StripSumScatter::PlotScatters() {
  // Display windows, recomputed here rather than in FillScatters: a cached run
  // skips the fill entirely, and leaving these at their constructor zeros made
  // SetRangeUser(0, 0) throw the zoom away on every run after the first.
  YBounds(m_yLo, m_yHi);
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
 
  std::lock_guard<std::mutex> lock(g_plot_mutex);
  for (Int_t reac = kReacMin; reac <= kReacMax; reac++) {
    TCanvas *c = PlottingUtils::GetConfiguredCanvas(kFALSE);
    // Display-only zoom: the histogram is built over the fixed build range, so
    // retuning these windows never forces a refill.
    Int_t ri = ReacIndex(reac);
    m_scatter[reac]->GetXaxis()->SetRangeUser(
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_DISPLAY_MIN,
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_DISPLAY_MAX);
    m_scatter[reac]->GetYaxis()->SetRangeUser(m_yLo[ri], m_yHi[ri]);
    PlottingUtils::ConfigureAndDraw2DHistogram(m_scatter[reac], c);
    m_scatter[reac]->GetYaxis()->SetTitleOffset(1.3);
    c->SetLeftMargin(0.18);
    PlottingUtils::SaveFigure(c,
                              Form("normsumE_reac%d_s%d_%d_vs_s%d_%d", reac,
                                   YLoOf(reac), YHiOf(reac), kXLo, kXHi),
                              "strip_sum_scatter", PlotSaveOptions::kLINEAR);
    delete c;
  }
}
 
// Diagnostic: the Grid #DeltaE of events that pass the cheap pre-tag cuts and
// are then rejected by the parity cut. Self-contained pass over the chains
// (not the fill's cuts), so it works on a cached dataset too. Two views of the
// same population are saved: the decoded Grid signal in a.u. ([0,1]) and the
// raw Grid branch in ADC. Purely visual: does not change what is tagged or the
// cache. Both are drawn on a log-y axis (PlotSaveOptions::kLOG).
void StripSumScatter::PlotParityRejectedGrid(
    const std::vector<Int_t> &run_order, std::map<Int_t, TChain *> &chains) {
  const StripSumScatterConfig &C = Constants::cfg.STRIP_SUM_SCATTER_CONFIG;
  if (!C.PLOT_PARITY_REJECTED_GRID)
    return;
  if (!(C.PARITY_ASYM_MAX > 0.0)) {
    std::cout << "strip-sum-scatter: PLOT_PARITY_REJECTED_GRID is set but "
                 "PARITY_ASYM_MAX <= 0; no parity cut to reject events, "
                 "skipping the grid diagnostic."
              << std::endl;
    return;
  }
 
  const Double_t grid_max_adc = Constants::ActiveGridMaxAdc();
 
  // a.u. view: the decoded grid (grid_adc / 16384 when calibrated), which is
  // what the rest of the analysis works in. ADC view: the raw trigger channel.
  TH1F *h_au = new TH1F("grid_parity_rejected_au",
                        ";Grid #DeltaE [a.u.];Counts", 400, 0.0, 1.0);
  TH1F *h_adc =
      new TH1F("grid_parity_rejected_adc", ";Grid #DeltaE [ADC];Counts", 400,
               0.0, grid_max_adc > 0.0 ? grid_max_adc : 16384.0);
  h_au->SetDirectory(nullptr);
  h_adc->SetDirectory(nullptr);
 
  Long64_t nRejected = 0;
  for (Int_t i = 0; i < Int_t(run_order.size()); i++) {
    TChain *chain = chains[run_order[i]];
    if (!chain)
      continue;
    EnergyView ev;
    ev.Attach(chain);
    EnableEventBranches(chain);
    const Long64_t n = chain->GetEntries();
    for (Long64_t j = 0; j < n; j++) {
      chain->GetEntry(j);
      ev.Decode();
      if (!AllStripsFired(ev) || IsPileup(ev) || IsNoise(ev))
        continue;
      if (C.REJECT_OFFBEAM && IsOffbeam(ev))
        continue;
      if (IsParityAsymmetric(ev)) {
        h_au->Fill(ev.grid);
        h_adc->Fill(ev.grid_adc);
        nRejected++;
      }
    }
    chain->ResetBranchAddresses();
  }
 
  std::cout << "strip-sum-scatter: parity-rejected grid diagnostic: "
            << nRejected << " events rejected by the parity cut." << std::endl;
  if (nRejected == 0) {
    delete h_au;
    delete h_adc;
    return;
  }
 
  {
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    TCanvas *c_au = PlottingUtils::GetConfiguredCanvas(kFALSE);
    PlottingUtils::ConfigureAndDrawHistogram(h_au, kBlue + 1);
    PlottingUtils::SaveFigure(c_au, "grid_parity_rejected_au",
                              "strip_sum_scatter", PlotSaveOptions::kLOG);
    delete c_au;
    TCanvas *c_adc = PlottingUtils::GetConfiguredCanvas(kFALSE);
    PlottingUtils::ConfigureAndDrawHistogram(h_adc, kBlue + 1);
    PlottingUtils::SaveFigure(c_adc, "grid_parity_rejected_adc",
                              "strip_sum_scatter", PlotSaveOptions::kLOG);
    delete c_adc;
  }
  delete h_au;
  delete h_adc;
}
 
void StripSumScatter::InteractiveOverlay(Int_t reac) {
  const Int_t kReacMin =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
  const Int_t kReacMax =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
  const Int_t kXLo = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_LO;
  const Int_t kXHi = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.X_HI;
  const Int_t kTracesPerRegion =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.TRACES_PER_CLASS;
 
  if (reac < kReacMin || reac > kReacMax) {
    std::cerr << "strip-sum-scatter: candidate reaction strip " << reac
              << " outside [" << kReacMin << "," << kReacMax
              << "]; skipping interactive overlay." << std::endl;
    return;
  }
  // A cut saved from a previous pass wins over prompting, so a run that has
  // already been decided once repeats without a person in the loop -- and
  // without a DISPLAY, since nothing interactive has to open.
  TCutG *cutAn = nullptr;
  TCutG *cutAa = nullptr;
  // Saved cuts -- drawn here earlier, or fitted by compute-regions -- win over
  // prompting, so a decided run repeats without a person in the loop.
  if (!Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REGION_CUT_REDRAW) {
    cutAn = LoadRegionCut("region_an", reac);
    cutAa = LoadRegionCut("region_aa", reac);
  }
  TCanvas *cutCanvas = nullptr;
 
  // The TApplication and its argv must outlive every GUI canvas made below,
  // which includes the trace canvases drawn after this block -- tearing it down
  // while those still exist faults in the ROOT paint path. So it is never
  // deleted: the process exits shortly after this function returns, and letting
  // it leak is the only way to keep it alive past the last canvas. The argv
  // storage is static for the same reason -- TApplication keeps the pointers we
  // hand it, and they must not dangle into a dead stack frame.
  static Int_t app_argc = 1;
  static char app_arg0[] = "strip-sum-scatter";
  static char *app_argv[] = {app_arg0};
 
  if (cutAn && cutAa) {
    std::cout << "  [region] loaded saved cuts for reac " << reac << std::endl;
  } else {
    delete cutAn;
    cutAn = nullptr;
    delete cutAa;
    cutAa = nullptr;
    if (!gSystem->Getenv("DISPLAY")) {
      std::cerr << "strip-sum-scatter: no saved region cuts for reac " << reac
                << " and no DISPLAY to draw them; skipping interactive "
                   "region-trace overlay (scatters already saved)."
                << std::endl;
      return;
    }
    // Intentionally not stored and never deleted -- see the note above.
    new TApplication("strip-sum-scatter", &app_argc, app_argv);
    gROOT->SetBatch(kFALSE);
 
    cutCanvas = new TCanvas("c_strip_sum_regions",
                            "Draw (a,n) then (a,a') regions", 900, 700);
    cutCanvas->SetLogz(kTRUE); // match the saved scatter's z-scale
    m_scatter[reac]->Draw("COLZ");
    cutCanvas->Update();
    cutAn = PromptCut(cutCanvas, "region_an", "(a,n)");
    cutAa = PromptCut(cutCanvas, "region_aa", "(a,a')");
    SaveRegionCuts(reac, cutAn, cutAa);
 
    cutCanvas->GetListOfPrimitives()->Remove(cutAn);
    cutCanvas->GetListOfPrimitives()->Remove(cutAa);
    gROOT->SetEditorMode();
    gSystem->ProcessEvents();
    // Stay in GUI mode for the rest of the run, and keep this canvas alive.
    //
    // PlotScatters() builds batch canvases earlier in this same process without
    // trouble, so batch canvas creation is fine on its own; what faults is
    // going back to batch once a TApplication and a GUI canvas exist.
    // TCanvas::Build() paints its border while constructing, and that paint
    // path is what dies in TPad::PaintBox. Deleting this canvas or nulling gPad
    // does not help -- only not re-entering batch does. The remaining canvases
    // open windows, which is fine for a command that already requires DISPLAY;
    // SaveFigure still writes them to disk, and the process exits shortly
    // after.
    cutCanvas->Clear();
    gSystem->ProcessEvents();
  }
 
  // ClusterVarHists(reac, cutAa, cutAn, "strip_sum_scatter");
 
  std::vector<TGraph *> tr_an, tr_aa, tr_beam;
  // Same selected events, raw (un-normalized) ADC -- one entry per normed
  // trace, kept in lock-step so the two overlays show the identical events.
  std::vector<TGraph *> tr_an_adc, tr_aa_adc, tr_beam_adc;
  // Same selected events again, Savitzky-Golay smoothed (normed a.u. space),
  // for the with-smoothing overlay -- also in lock-step with the raw normed
  // traces, so the two a.u. overlays show the identical events.
  const Bool_t kSkipSg =
      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.SKIP_SAVGOL_PLOTS;
  std::vector<TGraph *> tr_an_sg, tr_aa_sg, tr_beam_sg;
  UInt_t bit = (1u << ReacIndex(reac));
 
  for (Int_t k = 0; k < Int_t(m_reservoir.size()); k++) {
    if (Int_t(tr_an.size()) >= kTracesPerRegion &&
        Int_t(tr_aa.size()) >= kTracesPerRegion &&
        Int_t(tr_beam.size()) >= kTracesPerRegion)
      break;
    const TraceEvt &e = m_reservoir[k];
    if (e.beam_flat && Int_t(tr_beam.size()) < kTracesPerRegion) {
      tr_beam.push_back(TraceFromTotal(e.total));
      tr_beam_adc.push_back(TraceFromTotal(e.total_adc));
      if (!kSkipSg)
        tr_beam_sg.push_back(SmoothedTraceFromTotal(e.total));
      continue;
    }
    if (!(e.reac_mask & bit))
      continue;
    Double_t td[18];
    for (Int_t s = 0; s < 18; s++)
      td[s] = Double_t(e.total[s]);
    Double_t x = 0.0, y = 0.0;
    PlaneXY(td, reac, x, y);
    if (cutAn && Int_t(tr_an.size()) < kTracesPerRegion &&
        cutAn->IsInside(x, y)) {
      tr_an.push_back(TraceFromTotal(e.total));
      tr_an_adc.push_back(TraceFromTotal(e.total_adc));
      if (!kSkipSg)
        tr_an_sg.push_back(SmoothedTraceFromTotal(e.total));
    } else if (cutAa && Int_t(tr_aa.size()) < kTracesPerRegion &&
               cutAa->IsInside(x, y)) {
      tr_aa.push_back(TraceFromTotal(e.total));
      tr_aa_adc.push_back(TraceFromTotal(e.total_adc));
      if (!kSkipSg)
        tr_aa_sg.push_back(SmoothedTraceFromTotal(e.total));
    }
  }
 
  std::cout << "Sampled traces: beam=" << tr_beam.size()
            << " (a,a')=" << tr_aa.size() << " (a,n)=" << tr_an.size()
            << std::endl;
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.ALT_DECODE_REGION_TRACES)
    DrawAltDecodeRegionTraces(reac, cutAn, cutAa);
  DrawRegionTraces(Form("region_traces_reac%d", reac), "strip_sum_scatter",
                   tr_beam, tr_aa, tr_an, 0.6, 1.6, "#DeltaE [a.u.]");
  DrawRegionMeanTraces(Form("region_mean_traces_reac%d", reac),
                       "strip_sum_scatter", tr_beam, tr_aa, tr_an, 0.6, 1.6,
                       "#DeltaE [a.u.]");
  Double_t adc_y_lo = 0.0, adc_y_hi = 0.0;
  TraceYRange(tr_beam_adc, tr_aa_adc, tr_an_adc, adc_y_lo, adc_y_hi);
  DrawRegionTraces(Form("region_traces_reac%d_adc", reac), "strip_sum_scatter",
                   tr_beam_adc, tr_aa_adc, tr_an_adc, adc_y_lo, adc_y_hi,
                   "#DeltaE [ADC]");
  DrawRegionMeanTraces(Form("region_mean_traces_reac%d_adc", reac),
                       "strip_sum_scatter", tr_beam_adc, tr_aa_adc, tr_an_adc,
                       adc_y_lo, adc_y_hi, "#DeltaE [ADC]");
  if (!kSkipSg) {
    DrawRegionTraces(Form("region_traces_reac%d_sg", reac), "strip_sum_scatter",
                     tr_beam_sg, tr_aa_sg, tr_an_sg, 0.6, 1.6,
                     "#DeltaE [a.u.]");
    DrawRegionMeanTraces(Form("region_mean_traces_reac%d_sg", reac),
                         "strip_sum_scatter", tr_beam_sg, tr_aa_sg, tr_an_sg,
                         0.7, 1.3, "#DeltaE [a.u.]");
  }
 
  for (Int_t i = 0; i < Int_t(tr_an.size()); i++)
    delete tr_an[i];
  for (Int_t i = 0; i < Int_t(tr_aa.size()); i++)
    delete tr_aa[i];
  for (Int_t i = 0; i < Int_t(tr_beam.size()); i++)
    delete tr_beam[i];
  for (Int_t i = 0; i < Int_t(tr_an_adc.size()); i++)
    delete tr_an_adc[i];
  for (Int_t i = 0; i < Int_t(tr_aa_adc.size()); i++)
    delete tr_aa_adc[i];
  for (Int_t i = 0; i < Int_t(tr_beam_adc.size()); i++)
    delete tr_beam_adc[i];
  for (Int_t i = 0; i < Int_t(tr_an_sg.size()); i++)
    delete tr_an_sg[i];
  for (Int_t i = 0; i < Int_t(tr_aa_sg.size()); i++)
    delete tr_aa_sg[i];
  for (Int_t i = 0; i < Int_t(tr_beam_sg.size()); i++)
    delete tr_beam_sg[i];
 
  delete cutAn;
  delete cutAa;
}
 
void StripSumScatter::Run() {
  // Required before the threaded fill touches TChains from worker threads.
  ROOT::EnableThreadSafety();
  InitUtils::SetROOTPreferences(PlotSaveFormat::kPNG,
                                Paths::ResultsDir() + "/plots",
                                Paths::ResultsDir() + "/root_files");
  gROOT->SetBatch(kTRUE);
 
  std::vector<Int_t> run_order;
  std::map<Int_t, TChain *> chain_by_run = FileSet::GroupEventsByRun(run_order);
  if (run_order.empty()) {
    std::cerr << "strip-sum-scatter: no runs found" << std::endl;
    return;
  }
 
  // The jump gate and the upstream-beam tolerance are set in sigma of the
  // beam's noise, so the noise comes first: the fingerprint stamps the
  // thresholds they resolve to.
  Double_t jump_sigma[18], strip_sigma[18];
  if (!MeasureBeamNoise(chain_by_run[run_order[0]], jump_sigma, strip_sigma)) {
    std::cerr << "strip-sum-scatter: cannot measure the beam noise on run "
              << run_order[0] << std::endl;
    return;
  }
  SetJumpSigma(jump_sigma);
  SetStripSigma(strip_sigma);
  {
    const Int_t kReacMin =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MIN;
    const Int_t kReacMax =
        Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REACTION_STRIP_MAX;
    TString line =
        Form("strip-sum-scatter: jump noise sigma from run %d:", run_order[0]);
    for (Int_t reac = kReacMin; reac <= kReacMax; reac++)
      line += Form(" s%d %.4f", reac, jump_sigma[reac]);
    std::cout << line << std::endl;
    line =
        Form("strip-sum-scatter: strip noise sigma from run %d:", run_order[0]);
    for (Int_t strip = 1; strip < kReacMax; strip++)
      line += Form(" s%d %.4f", strip, strip_sigma[strip]);
    std::cout << line << std::endl;
    std::cout << Form("strip-sum-scatter: jump gate %.2f sigma -> %.4f at "
                      "strip %d, %.4f at strip %d",
                      Constants::cfg.STRIP_SUM_SCATTER_CONFIG.REAC_JUMP_NSIGMA,
                      JumpMin(kReacMin), kReacMin, JumpMin(kReacMax), kReacMax)
              << std::endl;
  }
 
  // Build fingerprint and try cache.
  TString fingerprint = BuildFingerprint(run_order, chain_by_run);
  TString cache_name = CacheName();
 
  Bool_t loaded = TryLoadCache(cache_name, fingerprint);
 
  if (!loaded) {
    FillScatters(run_order, chain_by_run);
    WriteCache(cache_name, fingerprint);
  }
 
  // Batch plotting (always done).
  PlotScatters();
 
  // Optional diagnostic: grid of events rejected by the parity cut. Reads the
  // Grid branch directly, so it works whether or not the scatter cache was
  // freshly filled. Self-gated on PLOT_PARITY_REJECTED_GRID.
  PlotParityRejectedGrid(run_order, chain_by_run);
 
  // Optional sim overlays.
  if (Constants::cfg.STRIP_SUM_SCATTER_CONFIG.RERUN_SIM) {
    SimOverlay();
    SimTraceOverlay();
  }
 
  // Interactive region-trace overlay (requires DISPLAY).
  Int_t reac = Constants::cfg.STRIP_SUM_SCATTER_CONFIG.CANDIDATE_REAC_STRIP;
  InteractiveOverlay(reac);
 
  // Cleanup chains.
  for (Int_t i = 0; i < Int_t(run_order.size()); i++)
    delete chain_by_run[run_order[i]];
}