From caeff0d7c5e5c26505c0bd294b98be8d981a8098 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Mon, 10 Aug 2026 14:27:52 -0400 Subject: [PATCH] modified: TrackRecon.C rehaul of the calibration dynamic so that any channels that dont have sufficient sx3 a1c2 stats get calibrated by lining up their centroids with those that do. Still a bit iffy on the bad anode formalism but keeping it in for now. modified: pc_energy_calibration.dat modified: pccal/fit_pc_energy_calibration.C --- TrackRecon.C | 139 ++++++++---- pc_energy_calibration.dat | 96 ++++----- pccal/fit_pc_energy_calibration.C | 348 +++++++++++++++++++++++++++--- 3 files changed, 457 insertions(+), 126 deletions(-) diff --git a/TrackRecon.C b/TrackRecon.C index b684beb..946b41e 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -33,6 +33,7 @@ Int_t colors[40] = { #include #include #include +#include #include #include // getpid(), for a unique per-process pc_calib_raw/ filename #include @@ -467,6 +468,28 @@ double pcIntercept[48]; double pcEnergySlope[48]; bool pcEnergyCalibLoaded = false; +// Wires currently suspected to have unreliable anode calibration -- factor >3x +// fit outliers piling up against the calibration ceiling (6, 19, 21, 22, 23), +// plus wire 12, which has too few calibration points to fit at all. Unlike +// a1c1_dead_anode above, these wires have plenty of raw statistics; they're +// just not trusted yet, so this is a testable toggle rather than a permanent +// mask. Set DISABLE_BAD_ANODE_WIRES=1 in the environment to exclude them from +// every anode cluster (A1C1/A1C2/A1C0) across the whole analysis, so the +// impact on downstream histograms can be compared against the default (off). +static const std::set badAnodeWires = {6, 12, 19, 21, 22, 23}; +bool excludeBadAnodeWires = false; // set in Begin() from DISABLE_BAD_ANODE_WIRES +inline bool isAnodeWireExcluded(int wire) +{ + return excludeBadAnodeWires && badAnodeWires.count(wire) > 0; +} +inline bool clusterHasExcludedAnode(const std::vector> &cl) +{ + for (const auto &w : cl) + if (isAnodeWireExcluded(std::get<0>(w))) + return true; + return false; +} + HistPlotter *plotter; bool HitNonZero; @@ -563,6 +586,16 @@ void TrackRecon::Begin(TTree * /*tree*/) if (getenv("CATHODE_GAIN")) cathode_gain = std::atof(getenv("CATHODE_GAIN")); + if (getenv("DISABLE_BAD_ANODE_WIRES")) + { + excludeBadAnodeWires = (std::atoi(getenv("DISABLE_BAD_ANODE_WIRES")) != 0); + std::cout << "DISABLE_BAD_ANODE_WIRES = " << excludeBadAnodeWires + << " -- excluding wires: "; + for (int w : badAnodeWires) + std::cout << w << " "; + std::cout << (excludeBadAnodeWires ? "(active)" : "(list defined but not active)") << std::endl; + } + const double *cfmin_src = a1c1_cfmin_17F; const double *k_src = a1c1_k_17F; const double *cfmin2_src = a1c1_cfmin2_17F; @@ -648,12 +681,24 @@ void TrackRecon::Begin(TTree * /*tree*/) } // ------------Load independent PC energy calibration (ADC -> dE_gas MeV)-------------- /// + // Gas gain depends on pressure (17F ran at 250 torr, 27Al at 350 torr for the alpha+gas + // campaigns), so a single pooled slope table isn't valid across datasets even when the + // eloss tables used to build each dataset's calibration points were themselves correct. + // Prefer a dataset-specific file (pc_energy_calibration_.dat); fall back to the + // old shared name (pc_energy_calibration.dat) if that doesn't exist, so this doesn't break + // for anyone who hasn't split their calibration by dataset yet. for (int i = 0; i < 48; i++) { pcEnergySlope[i] = 1.0; } { - std::ifstream pcEnergyFile("pc_energy_calibration.dat"); + std::string pcEnergyFilename = "pc_energy_calibration_" + dataset + ".dat"; + std::ifstream pcEnergyFile(pcEnergyFilename); + if (!pcEnergyFile.is_open()) + { + pcEnergyFilename = "pc_energy_calibration.dat"; + pcEnergyFile.open(pcEnergyFilename); + } if (pcEnergyFile.is_open()) { std::string line; @@ -670,7 +715,7 @@ void TrackRecon::Begin(TTree * /*tree*/) } pcEnergyFile.close(); pcEnergyCalibLoaded = true; - std::cout << "Loaded independent PC energy calibration from pc_energy_calibration.dat" + std::cout << "Loaded independent PC energy calibration from " << pcEnergyFilename << " -- populating PC_Events_calibrated" << std::endl; } } @@ -893,7 +938,7 @@ inline void pcEnergyCalibrationAccumulate(const std::vector &PC_Events, // second reference point to pick a cfrac branch), and A1C2's own crossover z, // though unambiguous, is no better than the Si hit's position once one exists. // There's no case where skipping the Si hit gives a more trustworthy point. - auto considerSi = [&](const Event &sievent, double phi_win) + auto considerSi = [&](const Event &sievent, double phi_win, bool isSX3) { if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) > phi_win) return; @@ -918,15 +963,20 @@ inline void pcEnergyCalibrationAccumulate(const std::vector &PC_Events, double Ex = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, alpha_source_mev, d_ex); if (!std::isfinite(Ee) || Ee <= 0.0 || !std::isfinite(Ex) || Ex < 0.0 || Ee <= Ex) return; - if (pcevent.Anodech >= 0 && pcevent.Anodech < 24) + // Anode: restricted to A1C2 topology, gated on SX3 coincidence only -- + // the trusted combination. QQQ-coincident and A1C1 points no longer + // contribute anode calibration data (cathode, below, is unaffected). + if (isSX3 && pcevent.multi1 == 1 && pcevent.multi2 == 2 && + pcevent.Anodech >= 0 && pcevent.Anodech < 24) pcCalibData[pcevent.Anodech].push_back({pcevent.Energy1, Ee - Ex}); + // Cathode: unchanged -- still A1C1, still both QQQ- and SX3-coincident. if (pcevent.multi2 == 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24) pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, Ee - Ex}); }; for (const auto &qqqevent : QQQ_Events) - considerSi(qqqevent, TMath::Pi() / 4.0); + considerSi(qqqevent, TMath::Pi() / 4.0, false); for (const auto &sx3event : SX3_Events) - considerSi(sx3event, TMath::Pi() / 3.0); + considerSi(sx3event, TMath::Pi() / 3.0, true); } } @@ -1538,6 +1588,8 @@ Bool_t TrackRecon::Process(Long64_t entry) for (const auto &aCluster : aClusters) { + if (clusterHasExcludedAnode(aCluster)) + continue; if (aCluster.size() == 2) { double ae0 = std::get<1>(aCluster[0]); @@ -1605,10 +1657,11 @@ Bool_t TrackRecon::Process(Long64_t entry) if ((pcEnergyCalibLoaded || doPCEnergyCalibration) && cClusters.empty()) { - const TVector3 source_pos_a1c0(beam_axis_x, beam_axis_y, source_vertex); for (const auto &aCl : aClusters) { - if (aCl.empty()) + if (aCl.size() != 1) // a1c0: exactly one anode wire -- same convention as + continue; // reaction_ax_core and miscHistograms_oneWire's a1c0 loops + if (clusterHasExcludedAnode(aCl)) continue; auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE"); auto apwire = std::get<0>(aPw); @@ -1647,43 +1700,23 @@ Bool_t TrackRecon::Process(Long64_t entry) if (pcEnergyCalibLoaded) { - double anodeCalibSum = 0.0; - for (const auto &w : aCl) - { - int wi = std::get<0>(w); - if (wi >= 0 && wi < 24) - anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w); - } + // aCl is guaranteed size 1 by the filter above, so anodeIdx unambiguously + // identifies the single wire this event's calibration applies to. + double anodeCalibSum = (anodeIdx >= 0 && anodeIdx < 24) + ? pcEnergySlope[anodeIdx] * std::get<1>(aCl[0]) + : 0.0; Event ev(pc, anodeCalibSum, -1.0, apTSMaxE, -1.0); - ev.multi1 = static_cast(aCl.size()); - ev.multi2 = 0; // no cathode -> a{n}c0 topology in pcCalibratedHistograms + ev.multi1 = 1; + ev.multi2 = 0; // no cathode -> a1c0 topology in pcCalibratedHistograms ev.Anodech = anodeIdx; ev.Cathodech = -1; PC_Events_calibrated.push_back(ev); } - // Anode-wire calibration point -- source runs only. The fixed alpha-source - // energy is only valid there; proton-run A1C0 has no elastic tag to predict - // its energy, so it contributes to the display but not the fit. - if (doPCEnergyCalibration && source_run) - { - TVector3 ray_dir = (pc - source_pos_a1c0).Unit(); - TVector3 virt_out = source_pos_a1c0 + ray_dir * 2000.0; - PCCollect pcc = pcCollectionPath(source_pos_a1c0, virt_out); - if (pcc.ok) - { - double tc = pathLengthCm(source_pos_a1c0, virt_out); - double de = tc - pcc.guard_cm; // source -> guard wires, cm - double dx = tc - pcc.cathode_cm; // source -> cathode, cm - if (std::isfinite(de) && de > 0.0 && std::isfinite(dx) && dx > de) - { - double Ee = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, alpha_source_mev, de); - double Ex = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, alpha_source_mev, dx); - if (std::isfinite(Ee) && std::isfinite(Ex) && Ee > Ex && Ex >= 0.0) - pcCalibData[anodeIdx].push_back({apSumE, Ee - Ex}); - } - } - } + // NOTE: A1C0 no longer contributes calibration points here -- training is + // restricted to A1C2 gated on SX3 (see pcEnergyCalibrationAccumulate). This + // used to push a model-predicted (Ee-Ex) point per A1C0 hit into the same + // pcCalibData[] the fit reads, which bypassed that restriction entirely. } } @@ -1884,8 +1917,12 @@ void TrackRecon::Terminate() { gSystem->mkdir("pc_calib_raw", kTRUE); std::string runTypeTag = source_run ? "src_" : (ta_foil_run ? "ap_" : "other_"); - std::string tag = runTypeTag + (getenv("RUN_NUMBER") ? std::string("run") + getenv("RUN_NUMBER") - : dataset + "_pid" + std::to_string(getpid())); + // dataset is included unconditionally -- previously it was dropped whenever + // RUN_NUMBER was set (i.e. always, when launched from run_tr.sh), so a 27Al + // run and a 17F run at the same run number produced indistinguishable + // filenames and fit_pc_energy_calibration.C's dataset_filter could never + // actually separate them. + std::string tag = runTypeTag + dataset + (getenv("RUN_NUMBER") ? std::string("_run") + getenv("RUN_NUMBER") : std::string("_pid") + std::to_string(getpid())); std::string outname = "pc_calib_raw/points_" + tag + ".dat"; std::ofstream outfile(outname); outfile << std::scientific << std::setprecision(6); @@ -2033,15 +2070,21 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector &QQQ_ for (const auto &qqqevent : QQQ_Events) { - plotter->Fill2D("Calib_dE_AnodeE_vs_QQQE" + t, 400, 0, 10, 800, 0, 3, qqqevent.Energy1, pcevent.Energy1, "hCalibPC"); + plotter->Fill2D("Calib_dE_AnodeE_vs_QQQE" + t, 400, 0, 10, 800, 0, 2, qqqevent.Energy1, pcevent.Energy1, "hCalibPC"); + if (pcevent.Anodech >= 0 && pcevent.Anodech < 24) + plotter->Fill2D("Calib_dE_AnodeE_vs_QQQE" + t + "_anode" + std::to_string(pcevent.Anodech), + 400, 0, 10, 800, 0, 3, qqqevent.Energy1, pcevent.Energy1, "EdE_wire"); if (hasCathode) - plotter->Fill2D("Calib_dE_CathodeE_vs_QQQE" + t, 400, 0, 10, 800, 0, 3, qqqevent.Energy1, pcevent.Energy2, "hCalibPC"); + plotter->Fill2D("Calib_dE_CathodeE_vs_QQQE" + t, 400, 0, 10, 800, 0, 2, qqqevent.Energy1, pcevent.Energy2, "hCalibPC"); } for (const auto &sx3event : SX3_Events) { - plotter->Fill2D("Calib_dE_AnodeE_vs_SX3E" + t, 400, 0, 10, 800, 0, 3, sx3event.Energy1, pcevent.Energy1, "hCalibPC"); + plotter->Fill2D("Calib_dE_AnodeE_vs_SX3E" + t, 400, 0, 10, 800, 0, 2, sx3event.Energy1, pcevent.Energy1, "hCalibPC"); + if (pcevent.Anodech >= 0 && pcevent.Anodech < 24) + plotter->Fill2D("Calib_dE_AnodeE_vs_SX3E" + t + "_anode" + std::to_string(pcevent.Anodech), + 400, 0, 10, 800, 0, 3, sx3event.Energy1, pcevent.Energy1, "EdE_wire"); if (hasCathode) - plotter->Fill2D("Calib_dE_CathodeE_vs_SX3E" + t, 400, 0, 10, 800, 0, 3, sx3event.Energy1, pcevent.Energy2, "hCalibPC"); + plotter->Fill2D("Calib_dE_CathodeE_vs_SX3E" + t, 400, 0, 10, 800, 0, 2, sx3event.Energy1, pcevent.Energy2, "hCalibPC"); } } } @@ -3435,6 +3478,8 @@ void miscHistograms_oneWire(HistPlotter *plotter, const std::vector &QQQ_ { if (acluster.size() != 1) // this function is scoped to single-wire anode continue; // clusters -- same convention as a1c0 elsewhere + if (clusterHasExcludedAnode(acluster)) + continue; auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(acluster, "ANODE"); // if(apSumE<6000) continue; int a_number = acluster.size(); @@ -3549,7 +3594,7 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector &QQQ_Ev { TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz); TVector3 rv = beamVertex(qqqevent.pos, x2 - qqqevent.pos); - if (beamPerp(rv) > 6.0 ) + if (beamPerp(rv) > 6.0) return; double th = (qqqevent.pos - rv).Theta(); double pl = pathLengthCm(qqqevent.pos, rv); @@ -4106,6 +4151,8 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector &Si_ { if (aCl.size() != 1) // a1c0: exactly one anode wire, no cathode -- same continue; // convention as pcevent.multi1==1 && multi2==0 elsewhere + if (clusterHasExcludedAnode(aCl)) + continue; auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE"); auto apwire = std::get<0>(aPw); double apSumE = std::get<1>(aPw); diff --git a/pc_energy_calibration.dat b/pc_energy_calibration.dat index ba26901..fd6dc2c 100644 --- a/pc_energy_calibration.dat +++ b/pc_energy_calibration.dat @@ -1,48 +1,48 @@ -0 4.290670e-05 0.000000e+00 -1 2.507557e-05 0.000000e+00 -2 3.625846e-05 0.000000e+00 -3 3.092832e-05 0.000000e+00 -4 2.692658e-05 0.000000e+00 -5 7.107454e-05 0.000000e+00 -6 7.774720e-05 0.000000e+00 -7 3.750049e-05 0.000000e+00 -8 4.165947e-05 0.000000e+00 -9 8.174412e-05 0.000000e+00 -10 4.354938e-05 0.000000e+00 -11 3.077641e-05 0.000000e+00 -12 1.000000e+00 0.000000e+00 -13 2.667676e-05 0.000000e+00 -14 2.587430e-05 0.000000e+00 -15 2.952973e-05 0.000000e+00 -16 3.347797e-05 0.000000e+00 -17 2.047205e-05 0.000000e+00 -18 1.628894e-05 0.000000e+00 -19 9.460420e-05 0.000000e+00 -20 3.718270e-05 0.000000e+00 -21 1.936310e-04 0.000000e+00 -22 1.061558e-04 0.000000e+00 -23 1.386274e-04 0.000000e+00 -24 3.032978e-05 0.000000e+00 -25 3.747258e-05 0.000000e+00 -26 3.486437e-05 0.000000e+00 -27 5.510136e-05 0.000000e+00 -28 4.472009e-05 0.000000e+00 -29 3.655720e-05 0.000000e+00 -30 3.584356e-05 0.000000e+00 -31 3.508027e-05 0.000000e+00 -32 3.613794e-05 0.000000e+00 -33 2.886789e-05 0.000000e+00 -34 3.048314e-05 0.000000e+00 -35 2.988481e-05 0.000000e+00 -36 3.735497e-05 0.000000e+00 -37 1.406118e-04 0.000000e+00 -38 3.302271e-05 0.000000e+00 -39 6.039465e-05 0.000000e+00 -40 4.551638e-05 0.000000e+00 -41 5.183895e-05 0.000000e+00 -42 3.797358e-05 0.000000e+00 -43 5.768681e-05 0.000000e+00 -44 7.210543e-05 0.000000e+00 -45 2.972753e-05 0.000000e+00 -46 2.822351e-05 0.000000e+00 -47 3.376403e-05 0.000000e+00 +0 2.270700e-05 0.000000e+00 1 +1 3.849247e-05 0.000000e+00 1 +2 2.478595e-05 0.000000e+00 1 +3 2.418520e-05 0.000000e+00 1 +4 2.473039e-05 0.000000e+00 1 +5 6.110646e-05 0.000000e+00 2 +6 4.489228e-05 0.000000e+00 1 +7 3.374396e-05 0.000000e+00 1 +8 3.370297e-05 0.000000e+00 1 +9 1.000000e+00 0.000000e+00 3 +10 3.876003e-05 0.000000e+00 1 +11 2.636256e-05 0.000000e+00 1 +12 1.000000e+00 0.000000e+00 3 +13 2.360190e-05 0.000000e+00 1 +14 2.012292e-05 0.000000e+00 1 +15 2.554712e-05 0.000000e+00 1 +16 2.085631e-05 0.000000e+00 1 +17 3.124900e-05 0.000000e+00 1 +18 4.364870e-05 0.000000e+00 1 +19 2.712940e-05 0.000000e+00 2 +20 3.140230e-05 0.000000e+00 1 +21 6.818787e-05 0.000000e+00 2 +22 3.142680e-05 0.000000e+00 2 +23 7.073949e-05 0.000000e+00 2 +24 3.024648e-05 0.000000e+00 1 +25 3.732136e-05 0.000000e+00 1 +26 3.486437e-05 0.000000e+00 1 +27 5.511389e-05 0.000000e+00 1 +28 4.472116e-05 0.000000e+00 1 +29 3.655778e-05 0.000000e+00 1 +30 3.579373e-05 0.000000e+00 1 +31 3.507960e-05 0.000000e+00 1 +32 3.613826e-05 0.000000e+00 1 +33 2.886743e-05 0.000000e+00 1 +34 3.031249e-05 0.000000e+00 1 +35 2.983830e-05 0.000000e+00 1 +36 3.735619e-05 0.000000e+00 1 +37 6.012566e-05 0.000000e+00 2 +38 3.303600e-05 0.000000e+00 1 +39 6.039465e-05 0.000000e+00 1 +40 4.511714e-05 0.000000e+00 1 +41 5.086580e-05 0.000000e+00 1 +42 3.765319e-05 0.000000e+00 1 +43 5.656884e-05 0.000000e+00 1 +44 1.196229e-04 0.000000e+00 2 +45 2.972753e-05 0.000000e+00 1 +46 2.829864e-05 0.000000e+00 1 +47 3.309860e-05 0.000000e+00 1 diff --git a/pccal/fit_pc_energy_calibration.C b/pccal/fit_pc_energy_calibration.C index d544218..d7a8907 100644 --- a/pccal/fit_pc_energy_calibration.C +++ b/pccal/fit_pc_energy_calibration.C @@ -11,17 +11,31 @@ #include #include #include +#include #include #include #include #include #include +#include // Optional dataset_filter (e.g., "17F" or "27Al") prevents mixing different // gas pressure/temperature environments which causes gain smearing. // Leave blank ("") to pool all files. -void fit_pc_energy_calibration(const std::string& dataset_filter = "") +// +// known_dead_wires: wire indices (0-23 anode, 24-47 cathode) known to be +// non-functional regardless of what their calibration data looks like. Whatever +// slope their own points would produce -- least-squares or robust -- isn't a real +// calibration and is excluded here rather than trusted. Update this list as more +// dead channels are identified; currently anode 9 and 12. +void fit_pc_energy_calibration(const std::string& dataset_filter = "", + const std::vector& known_dead_wires = {9, 12}) { + std::vector isDead(48, false); + for (int w : known_dead_wires) + if (w >= 0 && w < 48) + isDead[w] = true; + std::vector> pts[48]; // [wire] -> (ADC, dE_gas MeV) TSystemDirectory dir("pc_calib_raw", "pc_calib_raw"); @@ -35,6 +49,7 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") int nFiles = 0; long long nOverflowCut = 0; + long long nOverflowCutPerWire[48] = {0}; TIter next(files); TSystemFile *f; @@ -63,6 +78,7 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") pts[wire].push_back({adc, dE_gas}); } else { nOverflowCut++; + nOverflowCutPerWire[wire]++; } } } @@ -73,6 +89,214 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") << " run file(s) from pc_calib_raw/ (Filter: '" << dataset_filter << "')" << std::endl; std::cout << "fit_pc_energy_calibration: cut " << nOverflowCut << " points due to ADC >= 64k overflow." << std::endl; + std::cout << "fit_pc_energy_calibration: per-wire overflow breakdown (wire: cut / kept):" << std::endl; + for (int wire = 0; wire < 48; ++wire) + { + long long kept = static_cast(pts[wire].size()); + if (nOverflowCutPerWire[wire] == 0 && kept == 0) + continue; + double fracCut = (nOverflowCutPerWire[wire] + kept > 0) + ? 100.0 * nOverflowCutPerWire[wire] / (nOverflowCutPerWire[wire] + kept) + : 0.0; + std::cout << " " << (wire < 24 ? "anode " : "cathode ") << (wire < 24 ? wire : wire - 24) + << ": " << nOverflowCutPerWire[wire] << " cut / " << kept << " kept (" + << fracCut << "% cut)" << std::endl; + } + + // --- Adaptive lower-bound cut (anode only) ----------------------------- + // Raw ADC per wire can be genuinely bimodal (pileup, accidental coincidences, + // wrong-topology leakage) with a lower population that would otherwise bias + // both the least-squares slope and the median-based peak-matching fallback. + // Rather than one fixed ADC cutoff -- which can't be right for every wire at + // once when gains differ (the valley between the two populations sits at a + // genuinely different absolute ADC per wire) -- this finds each wire's own + // dominant peak from a coarse histogram of its own points, preferring the + // highest-ADC prominent peak over a lower one, and cuts at the valley just + // below it. Points below that valley are dropped before any fit/median runs. + const double kPeakProminenceFrac = 0.25; // a local max must reach this fraction + // of the tallest bin to count as a peak + const int kFloorHistBins = 60; + auto findAdaptiveFloor = [&](const std::vector &adcVals) -> double + { + if (adcVals.size() < 10) + return 0.0; // too few points to say anything about shape -- don't cut + double maxADC = *std::max_element(adcVals.begin(), adcVals.end()); + if (maxADC <= 0.0) + return 0.0; + std::vector hist(kFloorHistBins, 0); + double binW = maxADC / kFloorHistBins; + for (double v : adcVals) + { + int b = std::min(kFloorHistBins - 1, static_cast(v / binW)); + if (b >= 0) hist[b]++; + } + int tallest = *std::max_element(hist.begin(), hist.end()); + if (tallest <= 0) + return 0.0; + + // Scan from the highest-ADC bin down; the first local maximum that's + // prominent enough is taken as "the" peak. + int peakBin = -1; + for (int b = kFloorHistBins - 1; b >= 1; --b) + { + if (hist[b] < kPeakProminenceFrac * tallest) + continue; + bool isLocalMax = hist[b] >= hist[b - 1] && (b == kFloorHistBins - 1 || hist[b] >= hist[b + 1]); + if (isLocalMax) + { + peakBin = b; + break; + } + } + if (peakBin <= 0) + return 0.0; // no clear peak below the top edge -- nothing to cut against + + // Walk down from the peak to the valley: the floor is the bin where the + // count stops falling and starts rising again (the start of a lower + // population), or the histogram runs out. + int valleyBin = peakBin; + for (int b = peakBin - 1; b >= 0; --b) + { + if (hist[b] > hist[valleyBin]) + break; + valleyBin = b; + } + return valleyBin * binW; + }; + + std::vector> ptsFit[48]; // math uses this; pts[] stays raw for display + for (int wire = 0; wire < 48; ++wire) + ptsFit[wire] = pts[wire]; + + std::vector floorADC(48, 0.0); + for (int wire = 0; wire < 24; ++wire) // anode only + { + std::vector adcVals; + adcVals.reserve(pts[wire].size()); + for (const auto &p : pts[wire]) + adcVals.push_back(p.first); + double floor = findAdaptiveFloor(adcVals); + if (floor <= 0.0) + continue; + floorADC[wire] = floor; + size_t before = pts[wire].size(); + std::vector> kept; + kept.reserve(pts[wire].size()); + for (const auto &p : pts[wire]) + if (p.first >= floor) + kept.push_back(p); + ptsFit[wire] = std::move(kept); + std::cout << "fit_pc_energy_calibration: anode wire " << wire << " adaptive floor = " + << floor << " ADC -- kept " << ptsFit[wire].size() << "/" << before << " point(s)" << std::endl; + } + + // --- Pass 1: per-wire stats and the standard least-squares fit --- + // Split out from plotting so we can look at ALL wires' fitted slopes before + // deciding whether any individual one needs the robust fallback below. + std::vector slope_lsq(48, 1.0), maxX_arr(48, 0.0), maxY_arr(48, 0.0), n_arr(48, 0.0); + std::vector ok_arr(48, false); + + for (int wire = 0; wire < 48; ++wire) + { + double n = static_cast(ptsFit[wire].size()); + n_arr[wire] = n; + bool ok = (n >= 2); + + double maxX = 0.0, maxY = 0.0; // display bounds: from the raw (unfiltered) data + for (const auto &p : pts[wire]) + { + if (p.first > maxX) maxX = p.first; + if (p.second > maxY) maxY = p.second; + } + maxX_arr[wire] = maxX; + maxY_arr[wire] = maxY; + + double sxx = 0, sxy = 0; // fit sums: from the filtered data + for (const auto &p : ptsFit[wire]) + { + sxx += p.first * p.first; + sxy += p.first * p.second; + } + + if (ok) + { + if (std::isfinite(sxx) && std::abs(sxx) > 1e-12) + slope_lsq[wire] = sxy / sxx; + else + ok = false; + } + ok_arr[wire] = ok; + } + + auto medianOf = [](std::vector v) -> double + { + if (v.empty()) + return 0.0; + std::sort(v.begin(), v.end()); + size_t n = v.size(); + return (n % 2 == 1) ? v[n / 2] : 0.5 * (v[n / 2 - 1] + v[n / 2]); + }; + + // Median least-squares slope on each side (anode vs cathode gain structures differ, + // so they're compared separately, not pooled) -- the reference point for flagging + // any single wire's fit as an outlier. + std::vector anodeGoodSlopes, cathodeGoodSlopes; + for (int wire = 0; wire < 48; ++wire) + { + if (!ok_arr[wire] || isDead[wire]) + continue; + (wire < 24 ? anodeGoodSlopes : cathodeGoodSlopes).push_back(slope_lsq[wire]); + } + double medianAnodeSlope = medianOf(anodeGoodSlopes); + double medianCathodeSlope = medianOf(cathodeGoodSlopes); + + // A least-squares slope more than this factor away from its side's median (or a + // wire with too few points to fit at all) gets a peak-matched slope instead: that + // wire's own ADC peak lined up against the A1C2 consensus dE_gas peak pooled from + // every trusted wire on its side (computed below). Flagged wires are marked + // method=2 in the output file and "PEAK-MATCHED (a1c2 consensus)" in their + // diagnostic plot -- check those plots rather than trusting this blindly, since a + // genuinely multi-modal raw distribution needs a different fix, not a different average. + const double kOutlierRatioHi = 1.75; + const double kOutlierRatioLo = 1.0 / kOutlierRatioHi; + + std::cout << "fit_pc_energy_calibration: median anode slope = " << medianAnodeSlope + << ", median cathode slope = " << medianCathodeSlope << std::endl; + + // A1C2 consensus dE_gas peak: the target quantity (Ee - Ex from source-to-SX3 + // geometry) doesn't depend on which wire recorded the ADC -- it's the same + // physical prediction for every event. Pooling it across every wire that DID get + // a trusted direct fit ("calibrated using the a1c2 method") gives a far more + // precise reference than any single low-statistics wire's own handful of points + // could, which is what the wires below actually get matched against. + std::vector isOutlier(48, false); + for (int wire = 0; wire < 48; ++wire) + { + if (!ok_arr[wire] || isDead[wire]) + continue; + double medianSide = (wire < 24) ? medianAnodeSlope : medianCathodeSlope; + double ratio = (medianSide > 0.0) ? (slope_lsq[wire] / medianSide) : 1.0; + isOutlier[wire] = (medianSide > 0.0) && (ratio > kOutlierRatioHi || ratio < kOutlierRatioLo); + } + std::vector anodeConsensusPts, cathodeConsensusPts; + int nTrustedAnode = 0, nTrustedCathode = 0; + for (int wire = 0; wire < 48; ++wire) + { + if (!ok_arr[wire] || isDead[wire] || isOutlier[wire]) + continue; + auto &bucket = (wire < 24) ? anodeConsensusPts : cathodeConsensusPts; + for (const auto &p : ptsFit[wire]) + bucket.push_back(p.second); + (wire < 24 ? nTrustedAnode : nTrustedCathode)++; + } + double consensusAnodeDEgas = medianOf(anodeConsensusPts); + double consensusCathodeDEgas = medianOf(cathodeConsensusPts); + std::cout << "fit_pc_energy_calibration: A1C2 consensus anode dE_gas peak = " + << consensusAnodeDEgas << " MeV from " << anodeConsensusPts.size() + << " point(s) across " << nTrustedAnode << " trusted anode wire(s)" << std::endl; + std::cout << "fit_pc_energy_calibration: consensus cathode dE_gas peak = " + << consensusCathodeDEgas << " MeV from " << cathodeConsensusPts.size() + << " point(s) across " << nTrustedCathode << " trusted cathode wire(s)" << std::endl; // --- Setup ROOT Diagnostic Graphics --- gROOT->SetBatch(kTRUE); // Run silently without popping up windows @@ -90,49 +314,77 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") TCanvas *cCathodes = new TCanvas("cCathodes", "Cathode Calibrations", 1800, 1200); cCathodes->Divide(6, 4, 0.01, 0.01); - std::ofstream outfile("pc_energy_calibration.dat"); + std::string outFilename = dataset_filter.empty() + ? "pc_energy_calibration.dat" + : "pc_energy_calibration_" + dataset_filter + ".dat"; + std::ofstream outfile(outFilename); outfile << std::scientific << std::setprecision(6); - // --- Fit and Plot each wire --- + // --- Pass 2: finalize each wire's slope (lsq / robust fallback / identity / known-dead), plot, write --- for (int wire = 0; wire < 48; ++wire) { + double n = n_arr[wire]; + bool ok = ok_arr[wire]; double slope = 1.0, intercept = 0.0; - double n = static_cast(pts[wire].size()); - bool ok = (n >= 2); - - // 1. Calculate Standard Linear Least Squares Fit (y = mx) where intercept is forced to 0 - double maxX = 0.0, maxY = 0.0; - double sxx = 0, sxy = 0; // Only need sum(x^2) and sum(x*y) for fixed-0 intercept + int method = 0; // 0 = identity/no data, 1 = least-squares, 2 = peak-matched (a1c2 consensus), 3 = known dead - for (const auto &p : pts[wire]) + if (isDead[wire]) { - if (p.first > maxX) maxX = p.first; - if (p.second > maxY) maxY = p.second; - sxx += p.first * p.first; - sxy += p.first * p.second; + method = 3; + std::cerr << "fit_pc_energy_calibration: wire " << wire + << " is in known_dead_wires -- writing identity regardless of what its " + << n << " raw point(s) would otherwise fit to." << std::endl; } - - if (ok) + else if (ok && !isOutlier[wire]) { - if (std::isfinite(sxx) && std::abs(sxx) > 1e-12) + slope = slope_lsq[wire]; + method = 1; + } + else if (n >= 1) + { + // Either the least-squares slope is a >1.75x outlier vs its side's median, or + // there weren't even enough points (n<2) to attempt a direct fit at all. Either + // way: take this wire's own ADC peak (median of whatever points it has) and line + // it up with the A1C2-consensus dE_gas peak from the trusted wires on its side, + // rather than trusting a noisy few-point regression or falling back to identity. + std::vector xs; + xs.reserve(ptsFit[wire].size()); + for (const auto &p : ptsFit[wire]) + xs.push_back(p.first); + double medX = medianOf(xs); + double targetY = (wire < 24) ? consensusAnodeDEgas : consensusCathodeDEgas; + + if (medX > 1e-9 && targetY > 0.0) { - slope = sxy / sxx; - intercept = 0.0; // Forced mathematically + slope = targetY / medX; + method = 2; + std::cerr << "fit_pc_energy_calibration: wire " << wire << " " + << (ok ? "least-squares slope (" + std::to_string(slope_lsq[wire]) + ") is an outlier vs its side's median" + : "has too few points (" + std::to_string(static_cast(n)) + ") for a direct fit") + << " -- using peak-matched slope (own ADC median=" << medX + << " vs consensus dE_gas=" << targetY << ") = " << slope + << " instead. Check pc_calib_plots/wire_" << Form("%02d", wire) + << ".png to confirm this makes sense." << std::endl; } - else + else if (ok) { - ok = false; + // No usable consensus target (e.g. this side has no trusted wires at all) -- + // fall back to whatever least-squares gave, even though it was flagged. + slope = slope_lsq[wire]; + method = 1; } + // else: n>=1 but medX<=0 and no lsq fit either -- falls through to identity below. } - if (!ok) { - std::cerr << "fit_pc_energy_calibration: wire " << wire << " has too few points (" << n - << ") to fit -- writing identity (slope=1, intercept=0)" << std::endl; + if (method == 0 && !isDead[wire]) { + std::cerr << "fit_pc_energy_calibration: wire " << wire << " has " << n + << " point(s) and no usable consensus target -- writing identity (slope=1, intercept=0)" << std::endl; } - outfile << wire << " " << slope << " " << intercept << "\n"; + outfile << wire << " " << slope << " " << intercept << " " << method << "\n"; // 2. Generate and Fill 2D Density Histogram + double maxX = maxX_arr[wire], maxY = maxY_arr[wire]; if (maxX <= 0) maxX = 64000.0; if (maxY <= 0) maxY = 10.0; @@ -147,21 +399,41 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") } // 3. Setup Fit Line and Stats Box + // Known-dead wires never get a fit line drawn, even if they have raw points + // (wire 9 does) -- whatever structure is in that data isn't trusted as real + // calibration, but the raw scatter is still worth keeping for reference/crosstalk + // diagnosis. TF1 *fitLine = nullptr; TPaveText *pt = nullptr; - - if (n > 0) { + TLine *floorLine = nullptr; + bool drawFit = (n > 0) && !isDead[wire]; + + if (drawFit) { fitLine = new TF1(Form("fit_%d", wire), "[0]*x", 0, maxX * 1.05); // Formula is strictly y = m*x fitLine->SetParameter(0, slope); - fitLine->SetLineColor(kRed); + fitLine->SetLineColor(method == 2 ? kMagenta : kRed); fitLine->SetLineWidth(2); - pt = new TPaveText(0.15, 0.75, 0.55, 0.88, "NDC"); + pt = new TPaveText(0.15, 0.72, 0.55, 0.88, "NDC"); pt->SetFillColor(kWhite); pt->SetBorderSize(1); pt->AddText(Form("N = %.0f", n)); pt->AddText(Form("m = %.2e", slope)); pt->AddText("b = 0 (Fixed)"); + if (method == 2) + pt->AddText("PEAK-MATCHED (a1c2 consensus)"); + if (floorADC[wire] > 0.0) + pt->AddText(Form("floor cut @ %.0f ADC", floorADC[wire])); + } + // Adaptive floor marker: drawn whenever a cut was found, even for wires + // that ended up on the identity/dead path, so it's visible on every plot + // where it was computed, not just the ones that used it in a fit. + if (wire < 24 && floorADC[wire] > 0.0) + { + floorLine = new TLine(floorADC[wire], 0, floorADC[wire], maxY * 1.05); + floorLine->SetLineColor(kGreen + 2); + floorLine->SetLineStyle(2); + floorLine->SetLineWidth(2); } // ------------------------------------------------------------ @@ -170,13 +442,19 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") TCanvas cTemp("cTemp", "cTemp", 800, 600); cTemp.SetGridx(); cTemp.SetGridy(); - if (n > 0) { + if (drawFit) { h2->Draw("COLZ"); fitLine->Draw("SAME"); pt->Draw(); + if (floorLine) floorLine->Draw("SAME"); + } else if (isDead[wire]) { + h2->SetTitle(wName + (n > 0 ? " (KNOWN DEAD -- excluded from fit)" : " (KNOWN DEAD, NO DATA)")); + h2->Draw(n > 0 ? "COLZ" : ""); + if (floorLine) floorLine->Draw("SAME"); } else { h2->SetTitle(wName + " (DEAD/NO DATA)"); h2->Draw(); + if (floorLine) floorLine->Draw("SAME"); } cTemp.SaveAs(Form("pc_calib_plots/wire_%02d.png", wire)); @@ -187,11 +465,17 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") pad->SetGridx(); pad->SetGridy(); - if (n > 0) + if (drawFit) { h2->DrawClone("COLZ"); fitLine->DrawClone("SAME"); pt->DrawClone(); + if (floorLine) floorLine->DrawClone("SAME"); + } + else if (n > 0) + { + h2->DrawClone("COLZ"); // known-dead wire with data: show the raw scatter, no fit + if (floorLine) floorLine->DrawClone("SAME"); } else { @@ -207,6 +491,6 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "") fOut->Close(); - std::cout << "fit_pc_energy_calibration: wrote pc_energy_calibration.dat" << std::endl; + std::cout << "fit_pc_energy_calibration: wrote " << outFilename << std::endl; std::cout << "fit_pc_energy_calibration: individual high-res PNGs saved to pc_calib_plots/" << std::endl; } \ No newline at end of file