CalibrateBeam.cpp

Go to the documentation of this file.

#include "CalibrateBeam.hpp"
#include <Rtypes.h>
#include <TH2D.h>
#include <utility>
 
const Int_t kMaxChannels = 36;
// Raw total of one strip: strips 0 and 17 read once, strips 1-16 the sum
// of their two ends, indexed by strip - 1.
inline Double_t StripTotalAdc(const UShort_t *l, const UShort_t *r, UShort_t s0,
                              UShort_t s17, Int_t s) {
  if (s <= 0)
    return Double_t(s0);
  if (s >= 17)
    return Double_t(s17);
  return Double_t(l[s - 1]) + Double_t(r[s - 1]);
}

inline Int_t GatePartner(Int_t s) { return s <= 1 ? 0 : s - 1; }
const Long64_t kMinSamples = 200;
const Long64_t kSampleCap = 20000;
const Long64_t kPairCap = 100000;
 
// Defined below; forward-declared for the anchor reduction.
void RobustPeakSeed(const std::vector<Float_t> &v, Double_t &mode,
                    Double_t &sigma);
 
std::vector<ChannelCal> CalibrateBeam::BuildChannels() {
  std::vector<ChannelCal> chans;
  ChannelCal c{};
  c.name = "Strip0";
  c.side = 'S';
  c.strip = 0;
  chans.push_back(c);
  c.name = "Strip17";
  c.side = 'S';
  c.strip = 17;
  chans.push_back(c);
  for (Int_t s = 1; s <= 16; s++) {
    c.name = Form("L%d", s);
    c.side = 'L';
    c.strip = s;
    chans.push_back(c);
    c.name = Form("R%d", s);
    c.side = 'R';
    c.strip = s;
    chans.push_back(c);
  }
  c.name = "Cathode";
  c.side = 'C';
  c.strip = -1;
  chans.push_back(c);
  // The grid: anchored on its modal beam peak like a strip, so the analysis
  // can gate on it in beam units (GATE_AXIS_GRID).
  c.name = "Grid";
  c.side = 'G';
  c.strip = -1;
  chans.push_back(c);
  return chans;
}
 
Char_t LongSide(Int_t strip) { return (strip % 2 == 0) ? 'R' : 'L'; }
 
Bool_t IsCalibrated(const ChannelCal &c) { return c.fit_adc > 0; }
 
Double_t Gain(const ChannelCal &c) {
  return c.gain >= 0.0 ? c.gain : 1.0 / c.fit_adc;
}

BeamFit2D FindBeamGateStrips(const FileSpec &spec, Int_t sx, Int_t sy,
                             const TString &run_label,
                             const TString &plot_subdir,
                             Bool_t save_plot = kTRUE) {
  BeamFit2D out;
 
  TString sub = FileSet::EventsName(spec) + ".root";
  TFile *sf = IO::OpenForReading(sub);
  if (!sf || sf->IsZombie()) {
    if (sf)
      delete sf;
    return out;
  }
  TTree *tree = static_cast<TTree *>(sf->Get("events"));
  if (!tree) {
    sf->Close();
    delete sf;
    return out;
  }
  // Raw ADC, read before any calibration exists (this fit produces it), so read
  // the branches directly. strip1 total = L1 + R1; strip2 total = L2 + R2.
  UShort_t leftdE_adc[16], rightdE_adc[16], strip0_adc = 0, strip17_adc = 0;
  tree->SetBranchAddress("LeftdE", leftdE_adc);
  tree->SetBranchAddress("RightdE", rightdE_adc);
  tree->SetBranchAddress("Strip0dE", &strip0_adc);
  tree->SetBranchAddress("Strip17dE", &strip17_adc);

  const Int_t kBeamGateNBins = 1024;
  TH2F *h =
      new TH2F(Form("h2_gate_s%d_s%d_%s", sx, sy, run_label.Data()),
               Form(";Strip%d #DeltaE [ADC];Strip%d #DeltaE [ADC]", sx, sy),
               kBeamGateNBins, 0.0, Constants::ActiveStripEMaxAdc(),
               kBeamGateNBins, 0.0, Constants::ActiveStripEMaxAdc());
  h->SetDirectory(nullptr);
  Long64_t n = tree->GetEntries();
  // The same events unbinned (up to a cap), for the clipped moments below.
  const Long64_t kMaxPoints = 2000000;
  std::vector<std::pair<Float_t, Float_t>> pts;
  pts.reserve(std::min(n, kMaxPoints));
  for (Long64_t j = 0; j < n; j++) {
    tree->GetEntry(j);
    const Double_t tx =
        StripTotalAdc(leftdE_adc, rightdE_adc, strip0_adc, strip17_adc, sx);
    const Double_t ty =
        StripTotalAdc(leftdE_adc, rightdE_adc, strip0_adc, strip17_adc, sy);
    if (tx > 0.0 && ty > 0.0) {
      h->Fill(tx, ty);
      if (Long64_t(pts.size()) < kMaxPoints)
        pts.push_back(std::make_pair(Float_t(tx), Float_t(ty)));
    }
  }
  sf->Close();
  delete sf;
 
  if (h->GetEntries() < 100) {
    Constants::DetailErr() << "  " << run_label
                           << ": too few events for the strip " << sy
                           << " beam gate (strips " << sx << " vs " << sy << ")"
                           << std::endl;
    delete h;
    return out;
  }

  const Int_t kSeedHalfBins = 40;
  Moments2D m = BeamFitUtils::SeedSpotMoments(h, kSeedHalfBins);
  if (m.weight <= 0) {
    Constants::DetailErr() << "  " << run_label
                           << ": no bins above beam seed threshold"
                           << std::endl;
    delete h;
    return out;
  }
  // Re-center iteratively (moments inside ±2.5σ of the centroid); at high
  // rate (run84: 45 kHz) the ~2x-beam pileup blob passes the seed cut.
  const Int_t kMomentRefineIters = 4;
  const Double_t kMomentRefineNSigma = 2.5;
  const Double_t bw_x = h->GetXaxis()->GetBinWidth(1);
  const Double_t bw_y = h->GetYaxis()->GetBinWidth(1);
  Int_t bx, by, bz;
  h->GetMaximumBin(bx, by, bz);
  const Double_t peak_val = h->GetBinContent(bx, by);
  for (Int_t iter = 0; iter < kMomentRefineIters; iter++) {
    Int_t wlo_bx = std::max(
        1, h->GetXaxis()->FindBin(m.mu_x - kMomentRefineNSigma * m.sigma_x));
    Int_t whi_bx = std::min(
        h->GetNbinsX(),
        h->GetXaxis()->FindBin(m.mu_x + kMomentRefineNSigma * m.sigma_x));
    Int_t wlo_by = std::max(
        1, h->GetYaxis()->FindBin(m.mu_y - kMomentRefineNSigma * m.sigma_y));
    Int_t whi_by = std::min(
        h->GetNbinsY(),
        h->GetYaxis()->FindBin(m.mu_y + kMomentRefineNSigma * m.sigma_y));
    Moments2D m_ref = BeamFitUtils::ComputeMoments(
        h, wlo_bx, whi_bx, wlo_by, whi_by, 0.30 * peak_val, bw_x, bw_y);
    if (m_ref.weight <= 0)
      break;
    m = m_ref;
  }
  // The width BEAM_GATE_NSIGMA is in is a fitted sigma, finely binned
  // around the seed; over the full ADC range the spot is a few bins wide.
  BeamFitUtils::SpotFit spot = BeamFitUtils::FitSpotFromPoints(h, pts, m);
  out = spot.fit;
  Constants::Detail() << "  beam gate strip " << sy << " (strips " << sx
                      << " vs " << sy << "): "
                      << (spot.fit_used ? "fit"
                                        : "CLIPPED MOMENTS (fit failed)")
                      << " mu=(" << out.mu_x << "," << out.mu_y << ") sigma=("
                      << out.sigma_x << "," << out.sigma_y
                      << ") rho=" << out.rho << " chi2/ndf=" << spot.chi2_ndf
                      << std::endl;
 
  if (save_plot) {
    TCanvas *cv = PlottingUtils::GetConfiguredCanvas(kFALSE);
    PlottingUtils::ConfigureAndDraw2DHistogram(h, cv);
    BeamFitUtils::DrawEllipse(out, Constants::cfg.BEAM_GATE_NSIGMA,
                              kViolet + 2);
    if (Constants::SavePlots())
      PlottingUtils::SaveFigure(cv, Form("beam_gate_s%02d", sy),
                                plot_subdir + "/beam_gate",
                                PlotSaveOptions::kLINEAR);
    delete cv;
  }
  delete h;
  return out;
}
 
// Median by partial sort; v need not be sorted, but is only partially
// ordered afterwards, so a sort is still needed for other order statistics.
template <class T> inline Double_t Median(std::vector<T> &v) {
  if (v.empty())
    return 0.0;
  Int_t n = Int_t(v.size());
  std::nth_element(v.begin(), v.begin() + n / 2, v.end());
  Double_t med = Double_t(v[n / 2]);
  if (n % 2 == 0) {
    std::nth_element(v.begin(), v.begin() + n / 2 - 1, v.end());
    med = 0.5 * (med + Double_t(v[n / 2 - 1]));
  }
  return med;
}

void RobustPeakSeed(const std::vector<Float_t> &v, Double_t &mode,
                    Double_t &sigma) {
  mode = 0.0;
  sigma = 0.0;
  Int_t n = Int_t(v.size());
  if (n < 4)
    return;
  std::vector<Float_t> s(v);
  std::sort(s.begin(), s.end());
  Int_t i_lo = Int_t(0.05 * n);
  Int_t i_hi = Int_t(0.95 * n);
  Int_t i_q1 = Int_t(0.25 * n);
  Int_t i_q3 = Int_t(0.75 * n);
  if (i_hi >= n)
    i_hi = n - 1;
  Double_t p_lo = Double_t(s[i_lo]);
  Double_t p_hi = Double_t(s[i_hi]);
  Double_t med = Double_t(s[n / 2]);
  sigma = (Double_t(s[i_q3]) - Double_t(s[i_q1])) / 1.349;
  if (!(sigma > 0.0))
    sigma = 0.05 * (med > 0.0 ? med : 1.0);
  if (!(p_hi > p_lo)) {
    mode = med;
    return;
  }
  const Int_t nbins = 100;
  TH1F h("robust_peak_seed_h", "", nbins, p_lo, p_hi);
  h.SetDirectory(nullptr);
  for (Int_t j = 0; j < n; j++)
    if (Double_t(s[j]) >= p_lo && Double_t(s[j]) <= p_hi)
      h.Fill(Double_t(s[j]));
  mode = h.GetBinCenter(h.GetMaximumBin());
  if (!(mode > 0.0))
    mode = med;
}

struct StripPairSamples {
  std::vector<Float_t> l;
  std::vector<Float_t> r;
  // Beam-gated, long-triggered pairs: keep events with silent short end
  // (recorded 0) -- the population the decode sums; `l`/`r` need both ends.
  std::vector<Float_t> gated_long;
  std::vector<Float_t> gated_short;
};
 
// check_LR gain-match slice constants (37Cl_an_check_LR.ipynb); defined
// before CollectAnchorSamplesOneSubfile, which picks the shoulder slice.

void CollectAnchorSamplesOneSubfile(const FileSpec &spec,
                                    const std::vector<ChannelCal> &chans,
                                    const BeamFit2D gate[18],
                                    std::vector<std::vector<Float_t>> &samples,
                                    StripPairSamples pairs[18]) {
  Int_t n_chans = Int_t(chans.size());
  samples.assign(n_chans, std::vector<Float_t>());
  Bool_t any = kFALSE;
  for (Int_t s = 0; s <= 17; s++)
    any = any || gate[s].ok;
  if (!any)
    return;
 
  TString sub = FileSet::EventsName(spec) + ".root";
  TFile *sf = IO::OpenForReading(sub);
  if (!sf || sf->IsZombie()) {
    if (sf)
      sf->Close();
    return;
  }
  TTree *tree = static_cast<TTree *>(sf->Get("events"));
  if (!tree) {
    sf->Close();
    delete sf;
    return;
  }
  // Raw ADC, pre-calibration: the two ends of strips 1-16 and the
  // unsegmented strips; strip totals are L+R.
  UShort_t leftdE_adc[16], rightdE_adc[16], strip0_adc = 0, strip17_adc = 0;
  Short_t cathode_adc = 0, grid_adc = 0;
  tree->SetBranchAddress("LeftdE", leftdE_adc);
  tree->SetBranchAddress("RightdE", rightdE_adc);
  tree->SetBranchAddress("Strip0dE", &strip0_adc);
  tree->SetBranchAddress("Strip17dE", &strip17_adc);
  tree->SetBranchAddress("Cathode", &cathode_adc);
  tree->SetBranchAddress("Grid", &grid_adc);
 
  Long64_t n = tree->GetEntries();
  for (Long64_t j = 0; j < n; j++) {
    tree->GetEntry(j);
    // (L, R) pairs for strips 1–16, UNGATED (check_LR runs on all events):
    // a beam gate would remove the shoulder events (long side ~300-350 ADC).
    if (pairs) {
      for (Int_t s = 1; s <= 16; s++) {
        Int_t lv = Int_t(leftdE_adc[s - 1]);
        Int_t rv = Int_t(rightdE_adc[s - 1]);
        if (lv > 0 && rv > 0 && Long64_t(pairs[s].l.size()) < kPairCap) {
          pairs[s].l.push_back(Float_t(lv));
          pairs[s].r.push_back(Float_t(rv));
        }
      }
    }
    // Which strips this event is beam-like in, each judged by its own gate.
    Bool_t pass[18] = {kFALSE};
    for (Int_t s = 1; s <= 17; s++) {
      if (!gate[s].ok)
        continue;
      const Double_t x = StripTotalAdc(leftdE_adc, rightdE_adc, strip0_adc,
                                       strip17_adc, GatePartner(s));
      const Double_t y =
          StripTotalAdc(leftdE_adc, rightdE_adc, strip0_adc, strip17_adc, s);
      if (x <= 0.0 || y <= 0.0)
        continue;
      pass[s] = BeamFitUtils::InEllipseXY(gate[s], x, y,
                                          Constants::cfg.BEAM_GATE_NSIGMA);
    }
    // Strip 0 shares strip 1's gate, so share its verdict: GatePartner(0) is
    // 0, and a self-test lands every event outside, zeroing its gain.
    pass[0] = pass[1];
    if (pairs) {
      for (Int_t s = 1; s <= 16; s++) {
        if (!pass[s] || Long64_t(pairs[s].gated_long.size()) >= kPairCap)
          continue;
        Bool_t l_is_long = (LongSide(s) == 'L');
        Int_t long_v =
            l_is_long ? Int_t(leftdE_adc[s - 1]) : Int_t(rightdE_adc[s - 1]);
        Int_t short_v =
            l_is_long ? Int_t(rightdE_adc[s - 1]) : Int_t(leftdE_adc[s - 1]);
        if (long_v > 0) {
          pairs[s].gated_long.push_back(Float_t(long_v));
          pairs[s].gated_short.push_back(Float_t(short_v > 0 ? short_v : 0));
        }
      }
    }
    for (Int_t i = 0; i < n_chans; i++) {
      if (Long64_t(samples[i].size()) >= kSampleCap)
        continue;
      const ChannelCal &c = chans[i];
      // The cathode and the grid have no strip of their own, so they ride on
      // strip 1's gate.
      const Int_t s = (c.side == 'C' || c.side == 'G') ? 1 : c.strip;
      if (s < 0 || s > 17 || !pass[s])
        continue;
      Int_t v = 0;
      if (c.side == 'S')
        v = Int_t(c.strip == 0 ? strip0_adc : strip17_adc);
      else if (c.side == 'L')
        v = Int_t(leftdE_adc[c.strip - 1]);
      else if (c.side == 'R')
        v = Int_t(rightdE_adc[c.strip - 1]);
      else if (c.side == 'C')
        v = Int_t(cathode_adc);
      else if (c.side == 'G')
        v = Int_t(grid_adc);
      if (v > 0)
        samples[i].push_back(Float_t(v));
    }
  }
  sf->Close();
  delete sf;
}

const Int_t kGmBins = 512;
const Double_t kRidgeShortMinFrac =
    0.15; // start of the short window, as a fraction of C_long
const Double_t kRidgeShortMaxFrac =
    0.5;                            // cap on short, as a fraction of C_long
const Double_t kRidgeBandLo = 0.30; // single-particle band, x C_long
const Double_t kRidgeBandHi = 1.45;
const Int_t kRidgeSlices = 60;   // at most this many slices of the window
const Int_t kRidgeMinSlices = 6; // at least this many, for sparse ridges
const Long64_t kRidgeMinPerSlice = 200;
// Absolute noise floor per slice when the per-slice minimum is relaxed for
// thinly-spread ridges: below this a slice median is dominated by shot noise.
const Long64_t kRidgeNoiseFloor = 20;
const Int_t kRidgeMinPts = 6;
// C_short/C_long is a preamp gain ratio, order unity; outside this range
// the ridge fit failed in a way the intercept test cannot see.
const Double_t kRidgeRatioLo = 0.50;
const Double_t kRidgeRatioHi = 2.00;
const Double_t kGmEsumLo = 0.8; // a.u. eSum peak search window (pass 2)
const Double_t kGmEsumHi = 2.5;

struct RidgeFit {
  std::vector<Double_t> x, y, ey; // slice centres, medians, median errors
  Double_t lo = 0.0, hi = 0.0;    // short-axis fit window
  Double_t slope = 0.0, intercept = 0.0;
  Double_t c_long = 0.0, c_short = 0.0;
  Bool_t fitted = kFALSE; // a line was fitted (it may still be rejected)
  // Diagnostic: how far a strip got before the fit, so a "not measurable"
  // log line can say which check rejected it instead of only slope/intercept.
  Long64_t n_gated = 0;        // pairs handed to RidgeShortAnchor
  Long64_t n_short_window = 0; // stayed inside the short window
  Long64_t n_long_band = 0;    // also inside the long band
  Int_t n_slices = 0;          // slices with >= kRidgeMinPerSlice
  Long64_t min_per_slice = 0;  // adaptive per-slice minimum used this strip
  const char *fail = "";       // the check that returned 0, if any
};

Bool_t TheilSenLine(const std::vector<Double_t> &x,
                    const std::vector<Double_t> &y, Double_t &slope,
                    Double_t &intercept) {
  const Int_t n = Int_t(x.size());
  if (n < 2)
    return kFALSE;
  std::vector<Double_t> slopes;
  slopes.reserve(n * (n - 1) / 2);
  for (Int_t i = 0; i < n; i++)
    for (Int_t j = i + 1; j < n; j++) {
      Double_t dx = x[j] - x[i];
      if (dx == 0.0)
        continue;
      slopes.push_back((y[j] - y[i]) / dx);
    }
  if (slopes.empty())
    return kFALSE;
  slope = Median(slopes);
  std::vector<Double_t> resid;
  resid.reserve(n);
  for (Int_t i = 0; i < n; i++)
    resid.push_back(y[i] - slope * x[i]);
  intercept = Median(resid);
  return kTRUE;
}

Double_t RidgeShortAnchor(const std::vector<Float_t> &v_short,
                          const std::vector<Float_t> &v_long, Double_t c_long,
                          Double_t &slope_out, Double_t &intercept_out,
                          RidgeFit &dbg) {
  slope_out = 0.0;
  intercept_out = 0.0;
  dbg.c_long = c_long;
  dbg.n_gated = Long64_t(v_short.size());
  if (c_long <= 0 || v_short.size() != v_long.size() || v_short.size() < 500) {
    dbg.fail = c_long <= 0                       ? "c_long<=0"
               : v_short.size() != v_long.size() ? "size mismatch"
                                                 : "fewer than 500 gated pairs";
    return 0.0;
  }
  const Double_t hi = kRidgeShortMaxFrac * c_long;
  const Double_t lo = kRidgeShortMinFrac * c_long;
  // Pairs inside the short window and the long band, kept so the slice count
  // can be chosen from how many there are before they are binned.
  std::vector<std::pair<Double_t, Double_t>> in_window;
  for (Int_t j = 0; j < Int_t(v_short.size()); j++) {
    Double_t sh = Double_t(v_short[j]), lg = Double_t(v_long[j]);
    if (sh < lo || sh >= hi)
      continue;
    dbg.n_short_window++;
    if (lg <= kRidgeBandLo * c_long || lg >= kRidgeBandHi * c_long)
      continue;
    dbg.n_long_band++;
    in_window.push_back(std::make_pair(sh, lg));
  }
  // Slice count from the population: kRidgeMinPerSlice per slice on average,
  // between kRidgeMinSlices (sparse censored ends) and kRidgeSlices.
  Int_t n_slices = Int_t(Long64_t(in_window.size()) / kRidgeMinPerSlice);
  n_slices = TMath::Max(kRidgeMinSlices, TMath::Min(kRidgeSlices, n_slices));
  std::vector<std::vector<Double_t>> slice(n_slices);
  for (Int_t j = 0; j < Int_t(in_window.size()); j++) {
    Int_t b = Int_t((in_window[j].first - lo) / (hi - lo) * n_slices);
    if (b >= 0 && b < n_slices)
      slice[b].push_back(in_window[j].second);
  }
  std::vector<Double_t> x, y, ey;
  {
    // Per-slice min: strict floor first; if kRidgeMinPts slices don't fill,
    // halve & retry down to the noise floor (the censored ends are sparse).
    Long64_t min_per_slice = kRidgeMinPerSlice;
    for (;;) {
      x.clear();
      y.clear();
      ey.clear();
      for (Int_t b = 0; b < n_slices; b++) {
        if (Long64_t(slice[b].size()) < min_per_slice)
          continue;
        std::sort(slice[b].begin(), slice[b].end());
        const Int_t nb = Int_t(slice[b].size());
        // The IQR reads the sorted slice, so it comes before the median,
        // which only partially orders it.
        const Double_t iqr = slice[b][nb * 3 / 4] - slice[b][nb / 4];
        const Double_t med = Median(slice[b]);
        x.push_back(lo + (b + 0.5) * (hi - lo) / n_slices);
        y.push_back(med);
        ey.push_back(1.253 * (iqr / 1.349) / TMath::Sqrt(Double_t(nb)));
      }
      if (Int_t(x.size()) >= kRidgeMinPts || min_per_slice <= kRidgeNoiseFloor)
        break;
      // Halve the floor and retry.
      min_per_slice = TMath::Max(kRidgeNoiseFloor, min_per_slice / 2);
    }
    dbg.min_per_slice = min_per_slice;
  }
  dbg.x = x;
  dbg.y = y;
  dbg.ey = ey;
  dbg.lo = lo;
  dbg.hi = hi;
  dbg.n_slices = Int_t(x.size());
  if (Int_t(x.size()) < kRidgeMinPts) {
    dbg.fail = "fewer than kRidgeMinPts filled slices";
    return 0.0;
  }
  // TheilSen: median of pairwise slopes; outlying slice medians (2/3-particle
  // bands, low-long background) can't bend the slope as a chi-square pol1.
  Double_t slope = 0.0, inter = 0.0;
  if (!TheilSenLine(x, y, slope, inter)) {
    dbg.fail = "degenerate slice set (TheilSenLine)";
    return 0.0;
  }
  slope_out = slope;
  intercept_out = inter;
  dbg.slope = slope;
  dbg.intercept = inter;
  dbg.fitted = kTRUE;
  if (slope < 0)
    dbg.c_short = -inter / slope;
  if (slope >= 0) {
    dbg.fail = "slope >= 0";
    return 0.0;
  }
  // The line must cross the long axis near the modal seed (the short-end
  // offset adds up to ~20% on the R-chain); a large departure means a miss.
  if (inter < 0.80 * c_long || inter > 1.35 * c_long) {
    dbg.fail = "intercept outside [0.80, 1.35]*C_long";
    return 0.0;
  }
  // The line's short-axis crossing, as the "measurable" flag; the caller
  // works from slope and intercept directly.
  return -inter / slope;
}

void SaveRidgeFitPlots(const StripPairSamples pairs[18], const RidgeFit dbg[18],
                       const TString &plot_subdir) {
  TString subdir = plot_subdir + "/ridge";
  for (Int_t s = 1; s <= 16; s++) {
    const RidgeFit &d = dbg[s];
    const StripPairSamples &p = pairs[s];
    if (d.c_long <= 0 || p.gated_short.empty())
      continue;
    Bool_t l_is_long = (LongSide(s) == 'L');
    // A little past the fit window, so the cap the fit stops at is visible.
    Double_t xhi = (d.hi > 0 ? d.hi : kRidgeShortMaxFrac * d.c_long) * 1.35;
    TH2F *h = new TH2F(Form("h_ridge_s%d", s),
                       Form(";Strip %d Short (%c) #DeltaE [ADC];Strip %d Long "
                            "(%c) #DeltaE [ADC]",
                            s, l_is_long ? 'R' : 'L', s, l_is_long ? 'L' : 'R'),
                       200, 0.0, xhi, 200, 0.0, 1.60 * d.c_long);
    for (Int_t j = 0; j < Int_t(p.gated_short.size()); j++)
      h->Fill(Double_t(p.gated_short[j]), Double_t(p.gated_long[j]));
    TCanvas *cv = PlottingUtils::GetConfiguredCanvas(kFALSE);
    PlottingUtils::Configure2DHistogram(h, cv);
    h->Draw("COLZ");
    TGraphErrors *g = nullptr;
    if (!d.x.empty()) {
      g = new TGraphErrors(Int_t(d.x.size()), &d.x[0], &d.y[0], nullptr,
                           &d.ey[0]);
      g->SetMarkerStyle(20);
      g->SetMarkerSize(0.8);
      g->SetMarkerColor(kBlack);
      g->SetLineColor(kBlack);
      g->Draw("P SAME");
    }
    // Extended to short = 0 on purpose: where the line lands above the
    // long-only peak is the short-end offset times the slope.
    TLine *lf = nullptr;
    if (d.fitted) {
      lf = new TLine(0.0, d.intercept, xhi, d.intercept + d.slope * xhi);
      lf->SetLineColor(kViolet + 2);
      lf->SetLineWidth(2);
      lf->Draw();
    }
    // The short window the slices were taken from, so the threshold blob left
    // of it can be seen to be excluded.
    TLine *lw[2] = {nullptr, nullptr};
    if (d.hi > 0) {
      const Double_t ymax = 1.60 * d.c_long;
      lw[0] = new TLine(d.lo, 0.0, d.lo, ymax);
      lw[1] = new TLine(d.hi, 0.0, d.hi, ymax);
      for (Int_t k = 0; k < 2; k++) {
        lw[k]->SetLineColor(kGray + 2);
        lw[k]->SetLineStyle(2);
        lw[k]->Draw();
      }
    }
    if (Constants::SavePlots())
      PlottingUtils::SaveFigure(cv, Form("ridge_s%02d", s), subdir,
                                PlotSaveOptions::kLINEAR);
    delete cv;
    delete g;
    delete lf;
    delete lw[0];
    delete lw[1];
    delete h;
  }
}

Double_t GmFindPeak(const std::vector<Float_t> &v, Double_t lo, Double_t hi,
                    Double_t skip_frac) {
  if (v.empty() || hi <= lo)
    return 0.0;
  std::vector<Long64_t> h(kGmBins, 0);
  Double_t bw = (hi - lo) / kGmBins;
  for (Int_t j = 0; j < Int_t(v.size()); j++) {
    Double_t x = Double_t(v[j]);
    if (x < lo || x >= hi)
      continue;
    Int_t b = Int_t((x - lo) / bw);
    if (b >= 0 && b < kGmBins)
      h[b]++;
  }
  Int_t skip = Int_t(kGmBins * skip_frac);
  Int_t bmax = -1;
  Long64_t vmax = 0;
  for (Int_t b = skip; b < kGmBins; b++) {
    if (h[b] > vmax) {
      vmax = h[b];
      bmax = b;
    }
  }
  if (bmax < 0 || vmax <= 0)
    return 0.0;
  return lo + (bmax + 0.5) * bw;
}
 
// The L/R gain match state per strip: the anchors and short offset, the
// measured flags, and the found and median ratios and offsets.
struct GainMatchState {
  Bool_t matched[18] = {kFALSE};
  Bool_t ridge_ok[18] = {kFALSE};
  Double_t long_anchor_adc[18] = {0};
  Double_t short_anchor_adc[18] = {0};
  Double_t short_offset_adc[18] = {0};
  std::vector<Double_t> ratios_found[2];
  std::vector<Double_t> offsets_found[2];
  Double_t median_ratio[2] = {0.0, 0.0};
  Double_t median_offset[2] = {0.0, 0.0};
  Double_t global_ratio = 0.0;
};
 
// Each strip's L/R channel indices in the cal vector, -1 where absent.
static void IndexStripEnds(const std::vector<ChannelCal> &chans,
                           Int_t idx_l[18], Int_t idx_r[18]) {
  for (Int_t s = 0; s < 18; s++) {
    idx_l[s] = -1;
    idx_r[s] = -1;
  }
  for (Int_t i = 0; i < Int_t(chans.size()); i++) {
    if (chans[i].strip >= 1 && chans[i].strip <= 16) {
      if (chans[i].side == 'L')
        idx_l[chans[i].strip] = i;
      else if (chans[i].side == 'R')
        idx_r[chans[i].strip] = i;
    }
  }
}
 
// One strip's beam-gated ridge fit: the gate removes pile-up/junk while
// keeping off-centre crossings; the line is long = a - k short.
static void FitStripRidge(Int_t s, const StripPairSamples &p,
                          const ChannelCal &c_long, RidgeFit &dbg,
                          Double_t &peak_short, Double_t &peak_long,
                          Double_t &offset) {
  const std::vector<Float_t> &v_short = p.gated_short;
  const std::vector<Float_t> &v_long = p.gated_long;
  Double_t slope = 0.0, inter = 0.0;
  const Double_t crossing =
      RidgeShortAnchor(v_short, v_long, c_long.fit_adc, slope, inter, dbg);
  if (crossing > 0) {
    const Double_t k = -slope;
    peak_long = c_long.fit_adc;
    peak_short = peak_long / k;
    offset = (inter - peak_long) / k;
    Constants::Detail() << "  strip " << s
                        << " ridge slope=" << Form("%.3f", slope)
                        << " intercept=" << Form("%.1f", inter)
                        << " ADC; long anchor=" << Form("%.1f", peak_long)
                        << " (modal)  short anchor=" << Form("%.1f", peak_short)
                        << " ADC  short offset=" << Form("%.1f", offset)
                        << " ADC (" << Form("%.3f", offset / peak_long)
                        << " x C_long)" << std::endl;
  } else
    Constants::DetailErr() << "  strip " << s << ": ridge not measurable "
                           << "(slope=" << Form("%.3f", slope)
                           << ", intercept=" << Form("%.1f", inter) << ")"
                           << (dbg.fail[0] ? Form(" [%s]", dbg.fail) : "")
                           << " gated=" << Form("%lld", dbg.n_gated)
                           << " short_win=" << Form("%lld", dbg.n_short_window)
                           << " long_band=" << Form("%lld", dbg.n_long_band)
                           << " slices=" << dbg.n_slices
                           << " min_per_slice=" << dbg.min_per_slice
                           << std::endl;
}
 
// Pass 1: per strip the LONG anchor = the modal peak, the SHORT anchor
// C_long/k and the offset from the ridge line; the fallback fills gaps.
static void MeasureRidges(std::vector<ChannelCal> &chans,
                          const StripPairSamples pairs[18],
                          const Int_t idx_l[18], const Int_t idx_r[18],
                          RidgeFit *ridge_dbg, GainMatchState &m) {
  for (Int_t s = 1; s <= 16; s++) {
    if (idx_l[s] < 0 || idx_r[s] < 0)
      continue;
    Bool_t l_is_long = (LongSide(s) == 'L');
    ChannelCal &c_long = chans[l_is_long ? idx_l[s] : idx_r[s]];
    ChannelCal &c_short = chans[l_is_long ? idx_r[s] : idx_l[s]];
    if (!IsCalibrated(c_long)) {
      Constants::DetailErr()
          << "  strip " << s
          << ": long side uncalibrated; skipping L/R gain match" << std::endl;
      continue;
    }
    Double_t peak_short = 0.0, peak_long = 0.0, offset = 0.0;
    FitStripRidge(s, pairs[s], c_long, ridge_dbg[s], peak_short, peak_long,
                  offset);
    if (peak_short <= 0 || peak_long <= 0)
      continue;
    // C_short/C_long = -1/slope is a preamp-gain ratio, so order unity.
    // Catches ridges too flat for the intercept test: ratio absurd.
    Double_t ratio = peak_short / peak_long;
    if (ratio < kRidgeRatioLo || ratio > kRidgeRatioHi) {
      Constants::DetailErr()
          << "  strip " << s << ": ridge ratio " << Form("%.2f", ratio)
          << " outside [" << kRidgeRatioLo << ", " << kRidgeRatioHi
          << "]; rejecting anchor " << Form("%.1f", peak_short) << " ADC"
          << std::endl;
      continue;
    }
    m.long_anchor_adc[s] = peak_long;
    m.short_anchor_adc[s] = peak_short;
    m.short_offset_adc[s] = offset;
    m.ridge_ok[s] = kTRUE;
    m.ratios_found[s % 2].push_back(ratio);
    m.offsets_found[s % 2].push_back(offset / peak_long);
    c_short.ridge_ratio = ratio;
    c_short.ridge_offset = offset / peak_long;
  }
}
 
// The fallback medians, per parity (the L/R preamps differ, odd ~1.2 and
// even ~0.9) and globally over the ratios.
static void RidgeFallbackMedians(GainMatchState &m) {
  // Fall back on the median RATIO (preamp property, transfers across
  // strips; ADC anchors don't).
  std::vector<Double_t> all_ratios;
  for (Int_t par = 0; par < 2; par++) {
    std::vector<Double_t> &v = m.ratios_found[par];
    all_ratios.insert(all_ratios.end(), v.begin(), v.end());
    if (v.size() < 2)
      continue;
    m.median_ratio[par] = Median(v);
  }
  if (!all_ratios.empty())
    m.global_ratio = Median(all_ratios);
  Constants::Detail() << "  ridge ratio medians: odd="
                      << Form("%.3f", m.median_ratio[1])
                      << " even=" << Form("%.3f", m.median_ratio[0])
                      << " all=" << Form("%.3f", m.global_ratio) << std::endl;
  // The offset fallback, likewise per parity (it is a chain property), in
  // units of C_long; 0 where the parity measured none.
  for (Int_t par = 0; par < 2; par++) {
    std::vector<Double_t> &v = m.offsets_found[par];
    if (v.size() < 2)
      continue;
    m.median_offset[par] = Median(v);
  }
  Constants::Detail() << "  ridge offset medians (x C_long): odd="
                      << Form("%.3f", m.median_offset[1])
                      << " even=" << Form("%.3f", m.median_offset[0])
                      << std::endl;
}
 
// Apply the anchors: the parity-median fallback where the ridge was not
// measurable, then the gains and the short offset, with the detail lines.
static void ApplyGainMatch(std::vector<ChannelCal> &chans,
                           const Int_t idx_l[18], const Int_t idx_r[18],
                           GainMatchState &m) {
  for (Int_t s = 1; s <= 16; s++) {
    if (idx_l[s] < 0 || idx_r[s] < 0)
      continue;
    Bool_t l_is_long = (LongSide(s) == 'L');
    ChannelCal &c_long = chans[l_is_long ? idx_l[s] : idx_r[s]];
    ChannelCal &c_short = chans[l_is_long ? idx_r[s] : idx_l[s]];
    if (!IsCalibrated(c_long))
      continue;
    if (m.short_anchor_adc[s] <= 0) {
      // No ridge: the modal peak is the best long anchor available.
      m.long_anchor_adc[s] = c_long.fit_adc;
      Double_t r =
          m.median_ratio[s % 2] > 0 ? m.median_ratio[s % 2] : m.global_ratio;
      if (r <= 0) {
        Constants::DetailErr() << "  strip " << s
                               << ": no ridge and no ratio fallback; keeping "
                                  "independent gains"
                               << std::endl;
        continue;
      }
      m.short_anchor_adc[s] = r * m.long_anchor_adc[s];
      m.short_offset_adc[s] = m.median_offset[s % 2] * m.long_anchor_adc[s];
      Constants::Detail()
          << "  strip " << s << " short_anchor=ratio fallback "
          << Form("%.3f", r)
          << " x C_long = " << Form("%.1f", m.short_anchor_adc[s]) << " ADC"
          << (m.median_ratio[s % 2] > 0 ? "" : " (global, parity had none)")
          << "  short offset=parity median "
          << Form("%.3f", m.median_offset[s % 2])
          << " x C_long = " << Form("%.1f", m.short_offset_adc[s]) << " ADC"
          << std::endl;
    }
    c_long.gain = 1.0 / m.long_anchor_adc[s];
    c_short.gain = 1.0 / m.short_anchor_adc[s];
    c_short.offset_adc = m.short_offset_adc[s];
    m.matched[s] = kTRUE;
    Constants::Detail() << "  strip " << s << " L/R match: long anchor="
                        << Form("%.1f", m.long_anchor_adc[s])
                        << " ADC  short anchor="
                        << Form("%.1f", m.short_anchor_adc[s])
                        << " ADC  short offset="
                        << Form("%.1f", m.short_offset_adc[s]) << " ADC"
                        << (m.ridge_ok[s] ? "" : " (fallback)") << std::endl;
  }
}
 
// Pass 2: check, do not correct. On the pairs on the ridge, eSum with the
// short offset removed peaks at 1.0 a.u.; rescaling would hide a failure.
static void CheckSummedBeamPeaks(std::vector<ChannelCal> &chans,
                                 const StripPairSamples pairs[18],
                                 const Int_t idx_l[18], const Int_t idx_r[18],
                                 const GainMatchState &m) {
  Double_t esum_peak[18] = {0};
  for (Int_t s = 1; s <= 16; s++) {
    if (!m.matched[s])
      continue;
    const StripPairSamples &p = pairs[s];
    Bool_t l_is_long = (LongSide(s) == 'L');
    Double_t g_long = Gain(chans[l_is_long ? idx_l[s] : idx_r[s]]);
    Double_t g_short = Gain(chans[l_is_long ? idx_r[s] : idx_l[s]]);
    std::vector<Float_t> esum;
    esum.reserve(p.gated_long.size());
    const Double_t sh_lo = kRidgeShortMinFrac * m.long_anchor_adc[s];
    const Double_t sh_hi = kRidgeShortMaxFrac * m.long_anchor_adc[s];
    for (Int_t j = 0; j < Int_t(p.gated_long.size()); j++) {
      if (p.gated_short[j] < sh_lo || p.gated_short[j] >= sh_hi)
        continue;
      esum.push_back(Float_t(
          g_long * Double_t(p.gated_long[j]) +
          g_short * (Double_t(p.gated_short[j]) - m.short_offset_adc[s])));
    }
    esum_peak[s] = GmFindPeak(esum, kGmEsumLo, kGmEsumHi, 0.0);
  }
  for (Int_t s = 1; s <= 16; s++) {
    if (!m.matched[s])
      continue;
    if (esum_peak[s] <= 0) {
      Constants::DetailErr()
          << "  strip " << s << ": no summed beam peak in (" << kGmEsumLo
          << ", " << kGmEsumHi << ") a.u." << std::endl;
      continue;
    }
    Double_t dev = esum_peak[s] - 1.0;
    Constants::Detail() << "  strip " << s
                        << " summed beam peak=" << Form("%.4f", esum_peak[s])
                        << " a.u. (" << Form("%+.1f%%", 100.0 * dev) << ")"
                        << (TMath::Abs(dev) > 0.05 ? "  <-- check ridge fit"
                                                   : "")
                        << std::endl;
  }
}
 
void ComputeLRGainMatch(std::vector<ChannelCal> &chans,
                        const StripPairSamples pairs[18],
                        const TString &plot_subdir) {
  RidgeFit ridge_dbg[18];
  Int_t idx_l[18], idx_r[18];
  IndexStripEnds(chans, idx_l, idx_r);
 
  GainMatchState m;
  MeasureRidges(chans, pairs, idx_l, idx_r, ridge_dbg, m);
 
  // Between measure and fallback: the detail lines keep their order.
  SaveRidgeFitPlots(pairs, ridge_dbg, plot_subdir);
 
  RidgeFallbackMedians(m);
  ApplyGainMatch(chans, idx_l, idx_r, m);
  CheckSummedBeamPeaks(chans, pairs, idx_l, idx_r, m);
}

void ReduceToAnchors(std::vector<ChannelCal> &chans,
                     std::vector<std::vector<Float_t>> &samples,
                     const StripPairSamples pairs[18],
                     const TString &plot_subdir) {
  Int_t n_chans = Int_t(chans.size());
 
  for (Int_t i = 0; i < n_chans; i++) {
    ChannelCal &c = chans[i];
    std::vector<Float_t> &v = samples[i];
    c.n_samples = Long64_t(v.size());
 
    if (Long64_t(v.size()) < kMinSamples) {
      c.fit_adc = 0;
      continue;
    }
    if (c.side == 'C') {
      // Cathode: the median (asymmetric tail, no clean peak).
      c.fit_adc = Median(v);
    } else {
      // The modal bin, never the tail-biased mean.
      Double_t mode = 0.0, rsigma = 0.0;
      RobustPeakSeed(v, mode, rsigma);
      c.fit_adc = mode;
    }
    Constants::Detail() << "  " << c.name << " anchor[ADC]=" << c.fit_adc
                        << " (n=" << c.n_samples << ")" << std::endl;
  }
 
  // L/R match (strips 1-16): both anchors from the charge-sharing ridge, the
  // modal peak above being the seed and the fallback.
  if (pairs)
    ComputeLRGainMatch(chans, pairs, plot_subdir);
}

void WriteCalibrationTree(TFile *dst, const std::vector<ChannelCal> &chans,
                          const StripAlignmentResult *align) {
  dst->cd();
  if (TObject *old = dst->Get("calibration"))
    old->Delete();
  TTree *cal = new TTree("calibration", "Per-channel normMUSIC calibration");
  Float_t gain[kMaxChannels] = {0};
  Float_t fit_adc[kMaxChannels] = {0};
  Long64_t fit_n[kMaxChannels] = {0};
  Bool_t ok[kMaxChannels] = {0};
  // GainLeft[k]/GainRight[k] scale the two ends of strip k+1;
  // GainStrip0/GainStrip17 the unsegmented strips.
  Float_t gain_left[16] = {0}, gain_right[16] = {0};
  Float_t gain_strip0 = 0.0f, gain_strip17 = 0.0f;
  Float_t gain_cathode = 0.0f;
  Float_t gain_grid = 0.0f;
  // OffsetLeft[k]/OffsetRight[k] are subtracted before the gain, only when
  // that end fired (the short-end offset); zero on the long ends.
  Float_t offset_left[16] = {0}, offset_right[16] = {0};
  // Per-strip ridge ratio and offset (x C_long) measured in this subfile,
  // 0 where not measurable; AggregateRidgeRatiosForRun medians them.
  Float_t ridge_ratio[18] = {0};
  Float_t ridge_offset[18] = {0};
  Float_t long_anchor[18] = {0};
  Int_t n_actual = TMath::Min(Int_t(chans.size()), kMaxChannels);
  for (Int_t k = 0; k < n_actual; k++) {
    const ChannelCal &c = chans[k];
    ok[k] = IsCalibrated(c) || c.gain > 0;
    gain[k] = ok[k] ? Float_t(Gain(c)) : 0.0f;
    fit_adc[k] = Float_t(c.fit_adc);
    fit_n[k] = c.n_samples;
    if (c.strip >= 1 && c.strip <= 16 && (c.side == 'L' || c.side == 'R')) {
      if (c.ridge_ratio > 0) {
        ridge_ratio[c.strip] = Float_t(c.ridge_ratio);
        ridge_offset[c.strip] = Float_t(c.ridge_offset);
      }
      // The long anchor in use (the modal peak): AggregateRidgeRatiosForRun
      // rebuilds the short gain and offset from it.
      if (c.side == LongSide(c.strip) && ok[k] && gain[k] > 0)
        long_anchor[c.strip] = Float_t(1.0 / gain[k]);
    }
    if (c.side == 'S' && c.strip == 0)
      gain_strip0 = gain[k];
    else if (c.side == 'S' && c.strip == 17)
      gain_strip17 = gain[k];
    else if (c.side == 'L' && c.strip >= 1 && c.strip <= 16) {
      gain_left[c.strip - 1] = gain[k];
      offset_left[c.strip - 1] = Float_t(c.offset_adc);
    } else if (c.side == 'R' && c.strip >= 1 && c.strip <= 16) {
      gain_right[c.strip - 1] = gain[k];
      offset_right[c.strip - 1] = Float_t(c.offset_adc);
    } else if (c.side == 'C')
      gain_cathode = gain[k];
    else if (c.side == 'G')
      gain_grid = gain[k];
  }
  cal->Branch("Gain", gain, Form("Gain[%d]/F", kMaxChannels));
  cal->Branch("Ok", ok, Form("Ok[%d]/O", kMaxChannels));
  cal->Branch("FitADC", fit_adc, Form("FitADC[%d]/F", kMaxChannels));
  cal->Branch("FitN", fit_n, Form("FitN[%d]/L", kMaxChannels));
  cal->Branch("GainLeft", gain_left, "GainLeft[16]/F");
  cal->Branch("GainRight", gain_right, "GainRight[16]/F");
  cal->Branch("GainStrip0", &gain_strip0, "GainStrip0/F");
  cal->Branch("GainStrip17", &gain_strip17, "GainStrip17/F");
  cal->Branch("GainCathode", &gain_cathode, "GainCathode/F");
  cal->Branch("GainGrid", &gain_grid, "GainGrid/F");
  cal->Branch("OffsetLeft", offset_left, "OffsetLeft[16]/F");
  cal->Branch("OffsetRight", offset_right, "OffsetRight[16]/F");
  cal->Branch("RidgeRatio", ridge_ratio, "RidgeRatio[18]/F");
  cal->Branch("RidgeOffset", ridge_offset, "RidgeOffset[18]/F");
  cal->Branch("LongAnchor", long_anchor, "LongAnchor[18]/F");
 
  // Per-strip two-point alignment: StripFactor on every end, StripOffset on
  // the long end, so beam = 1 and pile-up = 2 on the strip total.
  Float_t strip_factor[18], strip_offset[18];
  for (Int_t s = 0; s < 18; s++) {
    strip_factor[s] = 1.0f;
    strip_offset[s] = 0.0f;
  }
  if (align && align->ok)
    for (Int_t s = 0; s < 18; s++) {
      strip_factor[s] = Float_t(align->factors[s]);
      strip_offset[s] = Float_t(align->offsets[s]);
    }
  cal->Branch("StripFactor", strip_factor, "StripFactor[18]/F");
  cal->Branch("StripOffset", strip_offset, "StripOffset[18]/F");
  cal->Fill();
  cal->Write("calibration", TObject::kOverwrite);
}

void WriteCalibrationToEvents(const FileSpec &spec,
                              const std::vector<ChannelCal> &chans,
                              const StripAlignmentResult *align) {
  TString events_subpath = FileSet::EventsName(spec) + ".root";
  TFile *f = IO::OpenForWriting(events_subpath, "UPDATE");
  if (!f || f->IsZombie()) {
    std::cerr << "Cannot open " << events_subpath << " to write calibration"
              << std::endl;
    if (f)
      delete f;
    return;
  }
  WriteCalibrationTree(f, chans, align);
  Constants::Detail() << "  wrote calibration into " << events_subpath
                      << std::endl;
  f->Close();
  delete f;
}
 
// Per-channel calibrated overlay for one subfile, via AttachCalSidecar (the
// same path downstream macros use). The sidecar must already be on disk.
void SaveDynamicRangeOverlay(const FileSpec &spec, const TString &plot_subdir,
                             const TString &file_label) {
  const Int_t kNStrips = 18;
  const Double_t emin = Constants::cfg.STRIP_DE_MIN_NORMED;
  const Double_t emax = Constants::cfg.STRIP_DE_MAX_NORMED;
  // 0.025 a.u. bins: a few per beam-peak width, so eighteen overlaid curves
  // read as curves rather than as noise.
  const Int_t nbins = TMath::Max(20, Int_t((emax - emin) / 0.025 + 0.5));
  TH1D *h[kNStrips];
  for (Int_t s = 0; s < kNStrips; s++) {
    h[s] = new TH1D(Form("h_dynrange_%s_S%d", file_label.Data(), s),
                    ";#DeltaE [a.u.];Counts", nbins, emin, emax);
    h[s]->SetDirectory(nullptr);
  }
  TString sub = FileSet::EventsName(spec) + ".root";
  TFile *sf = IO::OpenForReading(sub);
  if (!sf || sf->IsZombie()) {
    if (sf)
      sf->Close();
    for (Int_t s = 0; s < kNStrips; s++)
      delete h[s];
    return;
  }
  TTree *tree = static_cast<TTree *>(sf->Get("events"));
  if (!tree) {
    sf->Close();
    delete sf;
    for (Int_t s = 0; s < kNStrips; s++)
      delete h[s];
    return;
  }
  EnergyView ev;
  ev.Attach(tree);
  if (!ev.is_normed) {
    sf->Close();
    delete sf;
    for (Int_t s = 0; s < kNStrips; s++)
      delete h[s];
    return;
  }
  // The grid alone on the same axis, so its beam classes read against the
  // strips' range.
  TH1D *hg = new TH1D(Form("h_dynrange_%s_grid", file_label.Data()),
                      ";Grid #DeltaE [a.u.];Counts", nbins, emin, emax);
  hg->SetDirectory(nullptr);
  Long64_t n = tree->GetEntries();
  for (Long64_t j = 0; j < n; j++) {
    tree->GetEntry(j);
    ev.Decode();
    for (Int_t s = 0; s < kNStrips; s++) {
      Double_t v = ev.Total(s);
      if (v > 0)
        h[s]->Fill(v);
    }
    if (ev.grid > 0.0)
      hg->Fill(ev.grid);
  }
  sf->Close();
  delete sf;
  std::vector<Int_t> colors = PlottingUtils::GetDefaultColors();
  Double_t y_top = 0;
  for (Int_t s = 0; s < kNStrips; s++) {
    Double_t m = h[s]->GetMaximum();
    if (m > y_top)
      y_top = m;
  }
  TCanvas *cv = PlottingUtils::GetConfiguredCanvas(kFALSE);
  cv->SetRightMargin(0.20);
  Bool_t first = kTRUE;
  for (Int_t s = 0; s < kNStrips; s++) {
    Int_t color = colors[s % Int_t(colors.size())];
    h[s]->SetLineColor(color);
    h[s]->SetLineWidth(2);
    h[s]->SetMaximum(1.15 * y_top);
    h[s]->Draw(first ? "HIST" : "HIST SAME");
    first = kFALSE;
  }
  TLegend *leg = PlottingUtils::AddLegend(0.81, 0.99, 0.10, 0.95);
  for (Int_t s = 0; s < kNStrips; s++)
    leg->AddEntry(h[s], Form("S%d", s), "l");
  leg->Draw();
  if (Constants::SavePlots())
    PlottingUtils::SaveFigure(cv, "dynamic_range_check", plot_subdir,
                              PlotSaveOptions::kLOG);
  delete cv;
  delete leg;
  TCanvas *cg = PlottingUtils::GetConfiguredCanvas(kTRUE);
  PlottingUtils::ConfigureAndDrawHistogram(hg, kBlue + 1);
  if (Constants::SavePlots())
    PlottingUtils::SaveFigure(cg, "grid_spectrum", plot_subdir,
                              PlotSaveOptions::kLOG);
  delete cg;
  delete hg;
  for (Int_t s = 0; s < kNStrips; s++)
    delete h[s];
}

Double_t SmoothedPeakIn(TH1D *proj, Double_t lo, Double_t hi) {
  Int_t b_lo = proj->FindBin(lo);
  Int_t b_hi = proj->FindBin(hi);
  Int_t b_max = -1;
  Double_t val_max = 0;
  for (Int_t b = b_lo; b <= b_hi; b++) {
    Double_t v = proj->GetBinContent(b);
    if (v > val_max) {
      val_max = v;
      b_max = b;
    }
  }
  if (b_max < 0 || b_max == b_lo || b_max == b_hi)
    return 0.0;
  Double_t peak = proj->GetBinCenter(b_max);
  if (b_max > 1 && b_max < proj->GetNbinsX()) {
    const Double_t ym = proj->GetBinContent(b_max - 1);
    const Double_t y0 = proj->GetBinContent(b_max);
    const Double_t yp = proj->GetBinContent(b_max + 1);
    const Double_t denom = ym - 2.0 * y0 + yp;
    if (denom < 0.0) {
      Double_t shift = 0.5 * (ym - yp) / denom;
      shift = TMath::Max(-0.5, TMath::Min(0.5, shift));
      peak += shift * proj->GetBinWidth(b_max);
    }
  }
  return peak;
}

StripAlignmentResult FindStripCentroidAlignment(const FileSpec &spec,
                                                const TString &plot_subdir,
                                                const TString &file_label) {
  const Int_t kNStrips = 18;
  const Int_t kNHistBins = 2000; // 0.005 a.u. bins
  const Double_t kHistMin = 0.0;
  const Double_t kHistMax = 10.0;
  const Int_t kSmoothTimes = 5;
  const Double_t kBeamLo = 0.6, kBeamHi = 1.5;   // beam peak window
  const Double_t kPileLoX = 1.7, kPileHiX = 2.3; // 2-particle window, x beam
  StripAlignmentResult result;
  for (Int_t s = 0; s < kNStrips; s++) {
    result.factors[s] = 1.0;
    result.offsets[s] = 0.0;
    result.centroids[s] = 0.0;
    result.pileups[s] = 0.0;
  }
 
  TString sub = FileSet::EventsName(spec) + ".root";
  TFile *sf = IO::OpenForReading(sub);
  if (!sf || sf->IsZombie()) {
    if (sf)
      sf->Close();
    delete sf;
    return result;
  }
  TTree *tree = static_cast<TTree *>(sf->Get("events"));
  if (!tree) {
    sf->Close();
    delete sf;
    return result;
  }
  EnergyView ev;
  ev.Attach(tree);
  if (!ev.is_normed) {
    std::cerr << "  " << file_label
              << ": calibration tree not found -- skipping alignment"
              << std::endl;
    sf->Close();
    delete sf;
    return result;
  }
 
  TH2D *h2 = new TH2D(Form("h2_strip_align_%s", file_label.Data()),
                      ";Strip number;#DeltaE [a.u.]", kNStrips, -0.5,
                      kNStrips - 0.5, kNHistBins, kHistMin, kHistMax);
  h2->SetDirectory(nullptr);
 
  Long64_t n = tree->GetEntries();
  Long64_t n_used = 0;
  for (Long64_t j = 0; j < n; j++) {
    tree->GetEntry(j);
    ev.Decode();
    Bool_t any = kFALSE;
    for (Int_t s = 0; s < kNStrips; s++) {
      Double_t v = ev.Total(s);
      if (v <= 0)
        continue;
      h2->Fill(Double_t(s), v);
      any = kTRUE;
    }
    if (any)
      n_used++;
  }
  sf->Close();
  delete sf;
 
  Constants::Detail() << "  strip alignment: " << n_used << " events decoded"
                      << std::endl;
 
  Bool_t beam_ok[kNStrips] = {kFALSE};
  Bool_t pile_ok[kNStrips] = {kFALSE};
  Int_t valid_strips = 0;
 
  for (Int_t s = 0; s < kNStrips; s++) {
    Int_t bin_ix = s + 1;
    TH1D *proj = h2->ProjectionY(
        Form("hproj_align_%s_s%d", file_label.Data(), s), bin_ix, bin_ix);
    proj->SetDirectory(nullptr);
    Long64_t n_entries = Long64_t(proj->GetEntries());
    if (n_entries < kMinSamples) {
      Constants::DetailErr()
          << "  strip " << s << ": too few entries for alignment (" << n_entries
          << ")" << std::endl;
      delete proj;
      continue;
    }
    proj->Smooth(kSmoothTimes);
 
    Double_t beam = SmoothedPeakIn(proj, kBeamLo, kBeamHi);
    if (beam <= 0) {
      Constants::DetailErr()
          << "  strip " << s << ": no beam peak in [" << kBeamLo << ", "
          << kBeamHi << "] a.u.; left unaligned" << std::endl;
      delete proj;
      continue;
    }
    // The 2-particle peak sits near twice the beam; a peak on the window
    // edge is rejected by SmoothedPeakIn.
    Double_t pile = SmoothedPeakIn(proj, kPileLoX * beam, kPileHiX * beam);
    delete proj;
 
    result.centroids[s] = beam;
    beam_ok[s] = kTRUE;
    if (pile > 0) {
      result.pileups[s] = pile;
      pile_ok[s] = kTRUE;
      result.factors[s] = 1.0 / (pile - beam);
      result.offsets[s] = 1.0 - result.factors[s] * beam;
    } else {
      result.factors[s] = 1.0 / beam;
      result.offsets[s] = 0.0;
    }
    valid_strips++;
    Constants::Detail() << "  strip " << s << " beam=" << Form("%.4f", beam)
                        << " pileup="
                        << (pile > 0 ? Form("%.4f", pile)
                                     : "none (factor only)")
                        << " a.u.  factor=" << Form("%.4f", result.factors[s])
                        << " offset=" << Form("%+.4f", result.offsets[s])
                        << " (n=" << n_entries << ")" << std::endl;
  }
  result.ok = (valid_strips >= 4) ? kTRUE : kFALSE;
 
  // Diagnostic plot: the pre-alignment totals with the two peaks each strip
  // was aligned on (beam orange, pile-up violet).
  {
    TCanvas *cv = PlottingUtils::GetConfiguredCanvas(kFALSE);
    PlottingUtils::ConfigureAndDraw2DHistogram(h2, cv);
    TGraph *g_beam_plot = new TGraph(kNStrips);
    TGraph *g_pile_plot = new TGraph(kNStrips);
    Int_t nb = 0, npk = 0;
    for (Int_t s = 0; s < kNStrips; s++) {
      if (beam_ok[s])
        g_beam_plot->SetPoint(nb++, Double_t(s), result.centroids[s]);
      if (pile_ok[s])
        g_pile_plot->SetPoint(npk++, Double_t(s), result.pileups[s]);
    }
    g_beam_plot->Set(nb);
    g_pile_plot->Set(npk);
    if (nb > 0) {
      g_beam_plot->SetMarkerStyle(20);
      g_beam_plot->SetMarkerColor(kOrange);
      g_beam_plot->Draw("P SAME");
    }
    if (npk > 0) {
      g_pile_plot->SetMarkerStyle(21);
      g_pile_plot->SetMarkerColor(kViolet + 2);
      g_pile_plot->Draw("P SAME");
    }
    if (Constants::SavePlots())
      PlottingUtils::SaveFigure(cv, "strip_alignment_check", plot_subdir,
                                PlotSaveOptions::kLINEAR);
    delete cv;
    delete g_beam_plot;
    delete g_pile_plot;
  }
 
  delete h2;
 
  return result;
}
 
void CalibrateBeam::CalibrateBeamOneSubfile(
    const FileSpec &spec, const std::vector<ChannelCal> &chans_template) {
  TString file_label = FileSet::FileLabel(spec);
  TString plot_subdir = "beam_calibration/" + file_label;
  Constants::Detail() << "Beam calibration: " << file_label << std::endl;
 
  // One gate per strip: strip total vs its neighbour's. A strip whose gate
  // cannot be fitted is excluded from calibration, not the whole subfile.
  BeamFit2D gate[18];
  Int_t n_gates = 0;
  {
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    for (Int_t s = 1; s <= 17; s++) {
      gate[s] =
          FindBeamGateStrips(spec, GatePartner(s), s, file_label, plot_subdir);
      if (gate[s].ok)
        n_gates++;
      else
        Constants::DetailErr()
            << "  " << file_label << ": strip " << s
            << " beam gate failed; that strip is not calibrated here"
            << std::endl;
    }
    // Strip 0 shares strip 1's gate: there is no strip before it, so the
    // (0, 1) pair is the only one available to either.
    gate[0] = gate[1];
  }
  if (n_gates == 0) {
    std::cerr << "  " << file_label << ": every per-strip beam gate failed"
              << std::endl;
    return;
  }
 
  std::vector<ChannelCal> chans = chans_template;
  std::vector<std::vector<Float_t>> samples;
  StripPairSamples pairs[18];
  CollectAnchorSamplesOneSubfile(spec, chans, gate, samples, pairs);
  {
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    ReduceToAnchors(chans, samples, pairs, plot_subdir);
  }
 
  // Uncalibrated channels silently get gain 0 (drop out of strip totals), so
  // log which/why; grep "[uncalibrated long]" — long side should never starve.
  for (Int_t i = 0; i < Int_t(chans.size()); i++) {
    const ChannelCal &c = chans[i];
    if (IsCalibrated(c))
      continue;
    TString kind;
    if (c.side == 'C')
      kind = "cathode";
    else if (c.side == 'G')
      kind = "grid";
    else if (c.side == 'S')
      kind = "unsegmented";
    else
      kind = (c.side == LongSide(c.strip)) ? "long" : "short";
    TString why;
    if (c.n_samples < kMinSamples)
      why =
          Form("too few beam samples (%lld < %lld)", c.n_samples, kMinSamples);
    else
      why = Form("bad exp anchor (fit_adc=%.1f)", c.fit_adc);
    Constants::DetailErr() << "  [uncalibrated " << kind << "] " << c.name
                           << " -> gain 0; " << why
                           << " (beam n=" << c.n_samples << ")" << std::endl;
  }
 
  for (Int_t i = 0; i < Int_t(chans.size()); i++)
    if (IsCalibrated(chans[i]))
      Constants::Detail() << "  " << chans[i].name << " gain=" << Gain(chans[i])
                          << " a.u./ADC" << std::endl;
 
  StripAlignmentResult align;
 
  // Write initial calibration tree: per-channel gains. The alignment step
  // needs this on disk so EnergyView can decode events in a.u.
  WriteCalibrationToEvents(spec, chans, &align);
 
  // Per-strip two-point alignment on the strip total: beam to 1.0, pile-up
  // to 2.0 (factor on every end, offset on the long end).
  {
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    align = FindStripCentroidAlignment(spec, plot_subdir, file_label);
  }
  if (align.ok) {
    WriteCalibrationToEvents(spec, chans, &align);
  }
 
  {
    std::lock_guard<std::mutex> lock(g_plot_mutex);
    SaveDynamicRangeOverlay(spec, plot_subdir, file_label);
  }
 
  // One line per subfile whatever the sample: what was measured and what fell
  // back. The per-strip and per-channel lines above are detail.
  Int_t n_cal = 0, n_uncal = 0, n_ridge = 0, n_aligned = 0;
  for (Int_t i = 0; i < Int_t(chans.size()); i++) {
    if (IsCalibrated(chans[i]))
      n_cal++;
    else
      n_uncal++;
    if (chans[i].ridge_ratio > 0.0)
      n_ridge++;
  }
  for (Int_t s = 1; s <= 16; s++)
    if (align.ok && align.factors[s] != 1.0)
      n_aligned++;
  // One string, one write, so concurrent workers' lines do not interleave.
  const TString line =
      Form("[calibration] %s: gates %d/17, channels %d calibrated (%d not), "
           "ridge measured on %d/16 strips, alignment on %d/16",
           file_label.Data(), n_gates, n_cal, n_uncal, n_ridge, n_aligned);
  std::cout << line << std::endl;
}

void CalibrateBeam::AggregateRidgeRatiosForRun(
    Int_t run, const std::vector<FileSpec> &specs) {
  std::vector<Double_t> per_strip[18], per_strip_off[18];
  for (Int_t k = 0; k < Int_t(specs.size()); k++) {
    TString sub = FileSet::EventsName(specs[k]) + ".root";
    TFile *cf = IO::OpenForReading(sub);
    if (!cf || cf->IsZombie()) {
      if (cf)
        delete cf;
      continue;
    }
    TTree *t = static_cast<TTree *>(cf->Get("calibration"));
    if (!t || t->GetEntries() < 1 || !t->GetBranch("RidgeRatio") ||
        !t->GetBranch("RidgeOffset")) {
      cf->Close();
      delete cf;
      continue;
    }
    Float_t rr[18] = {0}, ro[18] = {0};
    t->SetBranchAddress("RidgeRatio", rr);
    t->SetBranchAddress("RidgeOffset", ro);
    t->GetEntry(0);
    for (Int_t s = 1; s <= 16; s++)
      if (rr[s] > 0) {
        per_strip[s].push_back(Double_t(rr[s]));
        per_strip_off[s].push_back(Double_t(ro[s]));
      }
    cf->Close();
    delete cf;
  }
 
  // The ratio and the offset (both in units of C_long) travel together: a
  // subfile that measured the ridge measured both.
  Double_t med[18] = {0}, med_off[18] = {0};
  std::cout << "Run " << run << " ridge-ratio aggregation:" << std::endl;
  for (Int_t s = 1; s <= 16; s++) {
    std::vector<Double_t> &v = per_strip[s];
    std::vector<Double_t> &vo = per_strip_off[s];
    if (v.size() < 3) {
      std::cerr << "  strip " << s << ": only " << v.size()
                << " subfiles measured the ridge; leaving per-subfile gains"
                << std::endl;
      continue;
    }
    std::sort(v.begin(), v.end());
    std::sort(vo.begin(), vo.end());
    Int_t m = Int_t(v.size());
    // The IQR reads the sorted vectors, so it comes before the medians,
    // which only partially order them.
    Double_t lo = v[m / 4], hi = v[(3 * m) / 4];
    Double_t lo_off = vo[m / 4], hi_off = vo[(3 * m) / 4];
    med[s] = Median(v);
    med_off[s] = Median(vo);
    std::cout << "  strip " << s << " ratio=" << Form("%.4f", med[s])
              << "  IQR " << Form("%.4f", lo) << "-" << Form("%.4f", hi)
              << "  offset=" << Form("%.4f", med_off[s]) << " x C_long  IQR "
              << Form("%.4f", lo_off) << "-" << Form("%.4f", hi_off)
              << "  from " << m << " subfiles" << std::endl;
  }
 
  Int_t n_rewritten = 0;
  for (Int_t k = 0; k < Int_t(specs.size()); k++) {
    TString sub = FileSet::EventsName(specs[k]) + ".root";
    TFile *cf = IO::OpenForWriting(sub, "UPDATE");
    if (!cf || cf->IsZombie()) {
      if (cf)
        delete cf;
      continue;
    }
    TTree *t = static_cast<TTree *>(cf->Get("calibration"));
    if (!t || t->GetEntries() < 1 || !t->GetBranch("LongAnchor") ||
        !t->GetBranch("OffsetLeft")) {
      cf->Close();
      delete cf;
      continue;
    }
    Float_t gl[16] = {0}, gr[16] = {0}, ol[16] = {0}, orr[16] = {0},
            la[18] = {0}, rr[18] = {0}, ro[18] = {0};
    t->SetBranchAddress("GainLeft", gl);
    t->SetBranchAddress("GainRight", gr);
    t->SetBranchAddress("OffsetLeft", ol);
    t->SetBranchAddress("OffsetRight", orr);
    t->SetBranchAddress("LongAnchor", la);
    t->SetBranchAddress("RidgeRatio", rr);
    t->SetBranchAddress("RidgeOffset", ro);
    t->GetEntry(0);
    Bool_t changed = kFALSE;
    for (Int_t s = 1; s <= 16; s++) {
      if (med[s] <= 0 || la[s] <= 0)
        continue;
      // Short anchor and offset from the run medians and this subfile's own
      // long anchor, so per-subfile gain drift is preserved.
      Double_t anchor = med[s] * Double_t(la[s]);
      if (anchor <= 0)
        continue;
      const Float_t g = Float_t(1.0 / anchor);
      const Float_t off = Float_t(med_off[s] * Double_t(la[s]));
      if (LongSide(s) == 'L') {
        gr[s - 1] = g;
        orr[s - 1] = off;
      } else {
        gl[s - 1] = g;
        ol[s - 1] = off;
      }
      rr[s] = Float_t(med[s]);
      ro[s] = Float_t(med_off[s]);
      changed = kTRUE;
    }
    if (changed) {
      TTree *nt = t->CloneTree(0);
      nt->Fill();
      cf->cd();
      nt->Write("calibration", TObject::kOverwrite);
      n_rewritten++;
    }
    cf->Close();
    delete cf;
  }
  std::cout << "  rewrote short gains and offsets in " << n_rewritten
            << " subfiles" << std::endl;
}
 
void CalibrateBeam::Run(const TString &file_label) {
  const TString project_root = Paths::DatasetDir();
  InitUtils::SetROOTPreferences(PlotSaveFormat::kPNG,
                                Paths::ResultsDir() + "/plots",
                                Paths::ResultsDir() + "/root_files");
  gROOT->SetBatch(kTRUE);
 
  std::vector<FileSpec> specs;
  if (file_label.IsNull()) {
    specs = FileSet::BuildProcessedFileSpecs();
    if (specs.empty()) {
      std::cerr << "No file specs from FileSet::BuildProcessedFileSpecs()"
                << std::endl;
      return;
    }
  } else {
    FileSpec s = FileSet::ResolveFileSpec(file_label);
    if (s.run < 0) {
      std::cerr << "Could not resolve file label '" << file_label << "'"
                << std::endl;
      return;
    }
    specs.push_back(s);
  }
 
  std::vector<ChannelCal> chans = CalibrateBeam::BuildChannels();
 
  Int_t n_specs = Int_t(specs.size());
 
  std::set<Int_t> runs;
  for (Int_t k = 0; k < n_specs; k++)
    runs.insert(specs[k].run);
 
  Int_t n_workers =
      TMath::Min(Int_t(std::thread::hardware_concurrency()), n_specs);
  n_workers = TMath::Min(n_workers, Constants::cfg.MAX_FUSED_WORKERS);
  if (n_workers < 1)
    n_workers = 1;
  std::cout << "calibrate-beam: " << n_specs << " subfiles on " << n_workers
            << " workers" << std::endl;
 
  std::queue<Int_t> work;
  for (Int_t k = 0; k < n_specs; k++)
    work.push(k);
  std::mutex work_mutex;
 
  std::vector<std::thread> workers;
  for (Int_t w = 0; w < n_workers; w++) {
    workers.emplace_back([&]() {
      while (true) {
        Int_t k;
        {
          std::lock_guard<std::mutex> lk(work_mutex);
          if (work.empty())
            return;
          k = work.front();
          work.pop();
        }
        // The first files of the epoch draw (a single named file always).
        Constants::SetPlotsThisFile(Constants::InPlotSample(k));
        CalibrateBeamOneSubfile(specs[k], chans);
      }
    });
  }
  for (Int_t w = 0; w < Int_t(workers.size()); w++)
    workers[w].join();
}