rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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.

Browse Source 
modified:   pc_energy_calibration.dat
	modified:   pccal/fit_pc_energy_calibration.C
This commit is contained in:
Vignesh Sitaraman 2026-08-10 14:27:52 -04:00
parent 6ff00bedc2
commit caeff0d7c5
3 changed files with 457 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

139
TrackRecon.C
View File

@ -33,6 +33,7 @@ Int_t colors[40] = {
#include <iostream> #include <iostream>
#include <sstream> #include <sstream>
#include <vector> #include <vector>
#include <set>
#include <array> #include <array>
#include <unistd.h> // getpid(), for a unique per-process pc_calib_raw/ filename #include <unistd.h> // getpid(), for a unique per-process pc_calib_raw/ filename
#include <map> #include <map>
@ -467,6 +468,28 @@ double pcIntercept[48];
double pcEnergySlope[48]; double pcEnergySlope[48];
bool pcEnergyCalibLoaded = false; 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<int> 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<std::tuple<int, double, double>> &cl)
{
for (const auto &w : cl)
if (isAnodeWireExcluded(std::get<0>(w)))
return true;
return false;
}
HistPlotter *plotter; HistPlotter *plotter;
bool HitNonZero; bool HitNonZero;
@ -563,6 +586,16 @@ void TrackRecon::Begin(TTree * /*tree*/)
if (getenv("CATHODE_GAIN")) if (getenv("CATHODE_GAIN"))
cathode_gain = std::atof(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 *cfmin_src = a1c1_cfmin_17F;
const double *k_src = a1c1_k_17F; const double *k_src = a1c1_k_17F;
const double *cfmin2_src = a1c1_cfmin2_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)-------------- /// // ------------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_<dataset>.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++) for (int i = 0; i < 48; i++)
{ {
pcEnergySlope[i] = 1.0; 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()) if (pcEnergyFile.is_open())
{ {
std::string line; std::string line;
@ -670,7 +715,7 @@ void TrackRecon::Begin(TTree * /*tree*/)
} }
pcEnergyFile.close(); pcEnergyFile.close();
pcEnergyCalibLoaded = true; 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; << " -- populating PC_Events_calibrated" << std::endl;
} }
} }
@ -893,7 +938,7 @@ inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events,
// second reference point to pick a cfrac branch), and A1C2's own crossover z, // 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. // 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. // 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) if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) > phi_win)
return; return;
@ -918,15 +963,20 @@ inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events,
double Ex = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, alpha_source_mev, d_ex); 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) if (!std::isfinite(Ee) || Ee <= 0.0 || !std::isfinite(Ex) || Ex < 0.0 || Ee <= Ex)
return; 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}); 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) if (pcevent.multi2 == 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24)
pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, Ee - Ex}); pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, Ee - Ex});
}; };
for (const auto &qqqevent : QQQ_Events) 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) 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) for (const auto &aCluster : aClusters)
{ {
if (clusterHasExcludedAnode(aCluster))
continue;
if (aCluster.size() == 2) if (aCluster.size() == 2)
{ {
double ae0 = std::get<1>(aCluster[0]); double ae0 = std::get<1>(aCluster[0]);
@ -1605,10 +1657,11 @@ Bool_t TrackRecon::Process(Long64_t entry)
if ((pcEnergyCalibLoaded || doPCEnergyCalibration) && cClusters.empty()) if ((pcEnergyCalibLoaded || doPCEnergyCalibration) && cClusters.empty())
{ {
const TVector3 source_pos_a1c0(beam_axis_x, beam_axis_y, source_vertex);
for (const auto &aCl : aClusters) 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; continue;
auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE"); auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE");
auto apwire = std::get<0>(aPw); auto apwire = std::get<0>(aPw);
@ -1647,43 +1700,23 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (pcEnergyCalibLoaded) if (pcEnergyCalibLoaded)
{ {
double anodeCalibSum = 0.0; // aCl is guaranteed size 1 by the filter above, so anodeIdx unambiguously
for (const auto &w : aCl) // identifies the single wire this event's calibration applies to.
{ double anodeCalibSum = (anodeIdx >= 0 && anodeIdx < 24)
int wi = std::get<0>(w); ? pcEnergySlope[anodeIdx] * std::get<1>(aCl[0])
if (wi >= 0 && wi < 24) : 0.0;
anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w);
}
Event ev(pc, anodeCalibSum, -1.0, apTSMaxE, -1.0); Event ev(pc, anodeCalibSum, -1.0, apTSMaxE, -1.0);
ev.multi1 = static_cast<int>(aCl.size()); ev.multi1 = 1;
ev.multi2 = 0; // no cathode -> a{n}c0 topology in pcCalibratedHistograms ev.multi2 = 0; // no cathode -> a1c0 topology in pcCalibratedHistograms
ev.Anodech = anodeIdx; ev.Anodech = anodeIdx;
ev.Cathodech = -1; ev.Cathodech = -1;
PC_Events_calibrated.push_back(ev); PC_Events_calibrated.push_back(ev);
} }
// Anode-wire calibration point -- source runs only. The fixed alpha-source // NOTE: A1C0 no longer contributes calibration points here -- training is
// energy is only valid there; proton-run A1C0 has no elastic tag to predict // restricted to A1C2 gated on SX3 (see pcEnergyCalibrationAccumulate). This
// its energy, so it contributes to the display but not the fit. // used to push a model-predicted (Ee-Ex) point per A1C0 hit into the same
if (doPCEnergyCalibration && source_run) // pcCalibData[] the fit reads, which bypassed that restriction entirely.
{
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});
}
}
}
} }
} }
@ -1884,8 +1917,12 @@ void TrackRecon::Terminate()
{ {
gSystem->mkdir("pc_calib_raw", kTRUE); gSystem->mkdir("pc_calib_raw", kTRUE);
std::string runTypeTag = source_run ? "src_" : (ta_foil_run ? "ap_" : "other_"); 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 is included unconditionally -- previously it was dropped whenever
: dataset + "_pid" + std::to_string(getpid())); // 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::string outname = "pc_calib_raw/points_" + tag + ".dat";
std::ofstream outfile(outname); std::ofstream outfile(outname);
outfile << std::scientific << std::setprecision(6); outfile << std::scientific << std::setprecision(6);
@ -2033,15 +2070,21 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
for (const auto &qqqevent : QQQ_Events) 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) 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) 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) 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<Event> &QQQ_
{ {
if (acluster.size() != 1) // this function is scoped to single-wire anode if (acluster.size() != 1) // this function is scoped to single-wire anode
continue; // clusters -- same convention as a1c0 elsewhere continue; // clusters -- same convention as a1c0 elsewhere
if (clusterHasExcludedAnode(acluster))
continue;
auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(acluster, "ANODE"); auto [apwire, apSumE, apMaxE, apTSMaxE] = pwinstance.GetPseudoWire(acluster, "ANODE");
// if(apSumE<6000) continue; // if(apSumE<6000) continue;
int a_number = acluster.size(); int a_number = acluster.size();
@ -3549,7 +3594,7 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
{ {
TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz); TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz);
TVector3 rv = beamVertex(qqqevent.pos, x2 - qqqevent.pos); TVector3 rv = beamVertex(qqqevent.pos, x2 - qqqevent.pos);
if (beamPerp(rv) > 6.0 ) if (beamPerp(rv) > 6.0)
return; return;
double th = (qqqevent.pos - rv).Theta(); double th = (qqqevent.pos - rv).Theta();
double pl = pathLengthCm(qqqevent.pos, rv); double pl = pathLengthCm(qqqevent.pos, rv);
@ -4106,6 +4151,8 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
{ {
if (aCl.size() != 1) // a1c0: exactly one anode wire, no cathode -- same if (aCl.size() != 1) // a1c0: exactly one anode wire, no cathode -- same
continue; // convention as pcevent.multi1==1 && multi2==0 elsewhere continue; // convention as pcevent.multi1==1 && multi2==0 elsewhere
if (clusterHasExcludedAnode(aCl))
continue;
auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE"); auto aPw = pwinstance.GetPseudoWire(aCl, "ANODE");
auto apwire = std::get<0>(aPw); auto apwire = std::get<0>(aPw);
double apSumE = std::get<1>(aPw); double apSumE = std::get<1>(aPw);

96
pc_energy_calibration.dat
View File

@ -1,48 +1,48 @@
0 4.290670e-05 0.000000e+00 0 2.270700e-05 0.000000e+00 1
1 2.507557e-05 0.000000e+00 1 3.849247e-05 0.000000e+00 1
2 3.625846e-05 0.000000e+00 2 2.478595e-05 0.000000e+00 1
3 3.092832e-05 0.000000e+00 3 2.418520e-05 0.000000e+00 1
4 2.692658e-05 0.000000e+00 4 2.473039e-05 0.000000e+00 1
5 7.107454e-05 0.000000e+00 5 6.110646e-05 0.000000e+00 2
6 7.774720e-05 0.000000e+00 6 4.489228e-05 0.000000e+00 1
7 3.750049e-05 0.000000e+00 7 3.374396e-05 0.000000e+00 1
8 4.165947e-05 0.000000e+00 8 3.370297e-05 0.000000e+00 1
9 8.174412e-05 0.000000e+00 9 1.000000e+00 0.000000e+00 3
10 4.354938e-05 0.000000e+00 10 3.876003e-05 0.000000e+00 1
11 3.077641e-05 0.000000e+00 11 2.636256e-05 0.000000e+00 1
12 1.000000e+00 0.000000e+00 12 1.000000e+00 0.000000e+00 3
13 2.667676e-05 0.000000e+00 13 2.360190e-05 0.000000e+00 1
14 2.587430e-05 0.000000e+00 14 2.012292e-05 0.000000e+00 1
15 2.952973e-05 0.000000e+00 15 2.554712e-05 0.000000e+00 1
16 3.347797e-05 0.000000e+00 16 2.085631e-05 0.000000e+00 1
17 2.047205e-05 0.000000e+00 17 3.124900e-05 0.000000e+00 1
18 1.628894e-05 0.000000e+00 18 4.364870e-05 0.000000e+00 1
19 9.460420e-05 0.000000e+00 19 2.712940e-05 0.000000e+00 2
20 3.718270e-05 0.000000e+00 20 3.140230e-05 0.000000e+00 1
21 1.936310e-04 0.000000e+00 21 6.818787e-05 0.000000e+00 2
22 1.061558e-04 0.000000e+00 22 3.142680e-05 0.000000e+00 2
23 1.386274e-04 0.000000e+00 23 7.073949e-05 0.000000e+00 2
24 3.032978e-05 0.000000e+00 24 3.024648e-05 0.000000e+00 1
25 3.747258e-05 0.000000e+00 25 3.732136e-05 0.000000e+00 1
26 3.486437e-05 0.000000e+00 26 3.486437e-05 0.000000e+00 1
27 5.510136e-05 0.000000e+00 27 5.511389e-05 0.000000e+00 1
28 4.472009e-05 0.000000e+00 28 4.472116e-05 0.000000e+00 1
29 3.655720e-05 0.000000e+00 29 3.655778e-05 0.000000e+00 1
30 3.584356e-05 0.000000e+00 30 3.579373e-05 0.000000e+00 1
31 3.508027e-05 0.000000e+00 31 3.507960e-05 0.000000e+00 1
32 3.613794e-05 0.000000e+00 32 3.613826e-05 0.000000e+00 1
33 2.886789e-05 0.000000e+00 33 2.886743e-05 0.000000e+00 1
34 3.048314e-05 0.000000e+00 34 3.031249e-05 0.000000e+00 1
35 2.988481e-05 0.000000e+00 35 2.983830e-05 0.000000e+00 1
36 3.735497e-05 0.000000e+00 36 3.735619e-05 0.000000e+00 1
37 1.406118e-04 0.000000e+00 37 6.012566e-05 0.000000e+00 2
38 3.302271e-05 0.000000e+00 38 3.303600e-05 0.000000e+00 1
39 6.039465e-05 0.000000e+00 39 6.039465e-05 0.000000e+00 1
40 4.551638e-05 0.000000e+00 40 4.511714e-05 0.000000e+00 1
41 5.183895e-05 0.000000e+00 41 5.086580e-05 0.000000e+00 1
42 3.797358e-05 0.000000e+00 42 3.765319e-05 0.000000e+00 1
43 5.768681e-05 0.000000e+00 43 5.656884e-05 0.000000e+00 1
44 7.210543e-05 0.000000e+00 44 1.196229e-04 0.000000e+00 2
45 2.972753e-05 0.000000e+00 45 2.972753e-05 0.000000e+00 1
46 2.822351e-05 0.000000e+00 46 2.829864e-05 0.000000e+00 1
47 3.376403e-05 0.000000e+00 47 3.309860e-05 0.000000e+00 1

346
pccal/fit_pc_energy_calibration.C
View File

@ -11,17 +11,31 @@
#include <TSystem.h> #include <TSystem.h>
#include <TROOT.h> #include <TROOT.h>
#include <TPaveText.h> #include <TPaveText.h>
#include <TLine.h>
#include <fstream> #include <fstream>
#include <iostream> #include <iostream>
#include <vector> #include <vector>
#include <iomanip> #include <iomanip>
#include <cmath> #include <cmath>
#include <algorithm>
// Optional dataset_filter (e.g., "17F" or "27Al") prevents mixing different // Optional dataset_filter (e.g., "17F" or "27Al") prevents mixing different
// gas pressure/temperature environments which causes gain smearing. // gas pressure/temperature environments which causes gain smearing.
// Leave blank ("") to pool all files. // 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<int>& known_dead_wires = {9, 12})
{ {
std::vector<bool> isDead(48, false);
for (int w : known_dead_wires)
if (w >= 0 && w < 48)
isDead[w] = true;
std::vector<std::pair<double, double>> pts[48]; // [wire] -> (ADC, dE_gas MeV) std::vector<std::pair<double, double>> pts[48]; // [wire] -> (ADC, dE_gas MeV)
TSystemDirectory dir("pc_calib_raw", "pc_calib_raw"); 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; int nFiles = 0;
long long nOverflowCut = 0; long long nOverflowCut = 0;
long long nOverflowCutPerWire[48] = {0};
TIter next(files); TIter next(files);
TSystemFile *f; TSystemFile *f;
@ -63,6 +78,7 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "")
pts[wire].push_back({adc, dE_gas}); pts[wire].push_back({adc, dE_gas});
} else { } else {
nOverflowCut++; 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; << " run file(s) from pc_calib_raw/ (Filter: '" << dataset_filter << "')" << std::endl;
std::cout << "fit_pc_energy_calibration: cut " << nOverflowCut std::cout << "fit_pc_energy_calibration: cut " << nOverflowCut
<< " points due to ADC >= 64k overflow." << std::endl; << " 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<long long>(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<double> &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<int> hist(kFloorHistBins, 0);
double binW = maxADC / kFloorHistBins;
for (double v : adcVals)
{
int b = std::min(kFloorHistBins - 1, static_cast<int>(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<std::pair<double, double>> ptsFit[48]; // math uses this; pts[] stays raw for display
for (int wire = 0; wire < 48; ++wire)
ptsFit[wire] = pts[wire];
std::vector<double> floorADC(48, 0.0);
for (int wire = 0; wire < 24; ++wire) // anode only
{
std::vector<double> 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<std::pair<double, double>> 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<double> slope_lsq(48, 1.0), maxX_arr(48, 0.0), maxY_arr(48, 0.0), n_arr(48, 0.0);
std::vector<bool> ok_arr(48, false);
for (int wire = 0; wire < 48; ++wire)
{
double n = static_cast<double>(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<double> 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<double> 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<bool> 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<double> 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 --- // --- Setup ROOT Diagnostic Graphics ---
gROOT->SetBatch(kTRUE); // Run silently without popping up windows 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); TCanvas *cCathodes = new TCanvas("cCathodes", "Cathode Calibrations", 1800, 1200);
cCathodes->Divide(6, 4, 0.01, 0.01); 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); 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) 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 slope = 1.0, intercept = 0.0;
double n = static_cast<double>(pts[wire].size()); int method = 0; // 0 = identity/no data, 1 = least-squares, 2 = peak-matched (a1c2 consensus), 3 = known dead
bool ok = (n >= 2);
// 1. Calculate Standard Linear Least Squares Fit (y = mx) where intercept is forced to 0 if (isDead[wire])
double maxX = 0.0, maxY = 0.0;
double sxx = 0, sxy = 0; // Only need sum(x^2) and sum(x*y) for fixed-0 intercept
for (const auto &p : pts[wire])
{ {
if (p.first > maxX) maxX = p.first; method = 3;
if (p.second > maxY) maxY = p.second; std::cerr << "fit_pc_energy_calibration: wire " << wire
sxx += p.first * p.first; << " is in known_dead_wires -- writing identity regardless of what its "
sxy += p.first * p.second; << n << " raw point(s) would otherwise fit to." << std::endl;
} }
else if (ok && !isOutlier[wire])
if (ok)
{ {
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<double> 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; slope = targetY / medX;
intercept = 0.0; // Forced mathematically 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<int>(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) { if (method == 0 && !isDead[wire]) {
std::cerr << "fit_pc_energy_calibration: wire " << wire << " has too few points (" << n std::cerr << "fit_pc_energy_calibration: wire " << wire << " has " << n
<< ") to fit -- writing identity (slope=1, intercept=0)" << std::endl; << " 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 // 2. Generate and Fill 2D Density Histogram
double maxX = maxX_arr[wire], maxY = maxY_arr[wire];
if (maxX <= 0) maxX = 64000.0; if (maxX <= 0) maxX = 64000.0;
if (maxY <= 0) maxY = 10.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 // 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; TF1 *fitLine = nullptr;
TPaveText *pt = nullptr; TPaveText *pt = nullptr;
TLine *floorLine = nullptr;
bool drawFit = (n > 0) && !isDead[wire];
if (n > 0) { if (drawFit) {
fitLine = new TF1(Form("fit_%d", wire), "[0]*x", 0, maxX * 1.05); // Formula is strictly y = m*x fitLine = new TF1(Form("fit_%d", wire), "[0]*x", 0, maxX * 1.05); // Formula is strictly y = m*x
fitLine->SetParameter(0, slope); fitLine->SetParameter(0, slope);
fitLine->SetLineColor(kRed); fitLine->SetLineColor(method == 2 ? kMagenta : kRed);
fitLine->SetLineWidth(2); 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->SetFillColor(kWhite);
pt->SetBorderSize(1); pt->SetBorderSize(1);
pt->AddText(Form("N = %.0f", n)); pt->AddText(Form("N = %.0f", n));
pt->AddText(Form("m = %.2e", slope)); pt->AddText(Form("m = %.2e", slope));
pt->AddText("b = 0 (Fixed)"); 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); TCanvas cTemp("cTemp", "cTemp", 800, 600);
cTemp.SetGridx(); cTemp.SetGridx();
cTemp.SetGridy(); cTemp.SetGridy();
if (n > 0) { if (drawFit) {
h2->Draw("COLZ"); h2->Draw("COLZ");
fitLine->Draw("SAME"); fitLine->Draw("SAME");
pt->Draw(); 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 { } else {
h2->SetTitle(wName + " (DEAD/NO DATA)"); h2->SetTitle(wName + " (DEAD/NO DATA)");
h2->Draw(); h2->Draw();
if (floorLine) floorLine->Draw("SAME");
} }
cTemp.SaveAs(Form("pc_calib_plots/wire_%02d.png", wire)); 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->SetGridx();
pad->SetGridy(); pad->SetGridy();
if (n > 0) if (drawFit)
{ {
h2->DrawClone("COLZ"); h2->DrawClone("COLZ");
fitLine->DrawClone("SAME"); fitLine->DrawClone("SAME");
pt->DrawClone(); 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 else
{ {
@ -207,6 +491,6 @@ void fit_pc_energy_calibration(const std::string& dataset_filter = "")
fOut->Close(); 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; std::cout << "fit_pc_energy_calibration: individual high-res PNGs saved to pc_calib_plots/" << std::endl;
} }