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

modified: Armory/PCZRecon.h added a2c0

Browse Source 
modified:   TrackRecon.C added a2c0, changed a1c1 cfrac parameters based on optimisation for the alpha and source data. Might need to be changed for the reaction though dfo rsome reason the apramteres don't quite makes sense
	modified:   eloss_calculations/Eloss.py added a scaling factor to particles with A>10 based to make the beams go further into the detector based on reaction data
	modified:   pc_energy_calibration.dat halved the anode dE cal to account for the misamtch between calib and rpedict for the source data
This commit is contained in:
Vignesh Sitaraman 2026-08-17 11:05:22 -04:00
parent c83557f30c
commit cf584a6d7a
4 changed files with 304 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
Armory/PCZRecon.h
View File

@ -2,8 +2,9 @@
#define PCZRecon_h #define PCZRecon_h
// PC Z-position reconstruction, one section per anode/cathode topology: // PC Z-position reconstruction, one section per anode/cathode topology:
// A1C0 (anode only), A1C1 (anode + single cathode, charge division), A1C2 // A1C0 (single anode wire only), A2C0 (two-wire anode cluster, no cathode --
// (anode + two cathodes, "step ladder" correction). Each topology gets one // same math as A1C0, see that section), A1C1 (anode + single cathode, charge
// division), A1C2 (anode + two cathodes, "step ladder" correction). Each topology gets one
// well-defined entry point instead of the math being split across files by // well-defined entry point instead of the math being split across files by
// historical accident (A1C0/A1C1 used to live in TrackRecon.C itself; A1C2's // historical accident (A1C0/A1C1 used to live in TrackRecon.C itself; A1C2's
// underlying model lives in the separately-shared PC_StepLadder_Correction.h // underlying model lives in the separately-shared PC_StepLadder_Correction.h
@ -83,6 +84,31 @@ inline TVector3 a1c0_hybrid_pcz(const std::pair<TVector3, TVector3> &apwire, dou
return pc; return pc;
} }
// ---------------------------------------------------------------------
// A2C0: two-wire anode-cluster (charge-shared), no-cathode position
// reconstruction
// ---------------------------------------------------------------------
//
// Same math as A1C0 above -- pseudowire + phi-minimization, no dither --
// just handed a genuine 2-wire energy-weighted pseudowire (GetPseudoWire
// over a 2-wire anode cluster) instead of a single real wire. a1c0_wirePos
// already treats its `apwire` argument as an opaque pseudowire pair
// regardless of how many physical wires went into it, so this is a
// documented, named entry point rather than new math: call sites can say
// what topology they mean instead of reusing a1c0_wirePos silently and
// trusting a comment to explain why.
//
// Deliberately no dithered twin (no a2c0_hybrid_pcz): A2C0 is meant to feed
// the reaction-analysis plots at its raw, undithered resolution, not stand
// in for a1c0_hybrid_pcz's benchmark-truth-comparison role. If a "genuine
// A2C0" BenchMark validation block is ever wanted (mirroring the existing
// aClusters.size()==1 && cClusters.size()==0 A1C0 block in TrackRecon.C),
// add one there rather than adding dithering here.
inline TVector3 a2c0_wirePos(const std::pair<TVector3, TVector3> &apwire, double phi, bool isQQQ)
{
return a1c0_wirePos(apwire, phi, isQQQ);
}
// --------------------------------------------------------------------- // ---------------------------------------------------------------------
// A1C1: single-anode + single-cathode charge-division position reconstruction // A1C1: single-anode + single-cathode charge-division position reconstruction
// --------------------------------------------------------------------- // ---------------------------------------------------------------------

291
TrackRecon.C
View File

@ -144,6 +144,16 @@ static const double a1c1_k_17F[7] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
static const double a1c1_cfmin_27Al[7] = {0.42, 0.42, 0.42, 0.40, 0.42, 0.43, 0.43}; static const double a1c1_cfmin_27Al[7] = {0.42, 0.42, 0.42, 0.40, 0.42, 0.43, 0.43};
static const double a1c1_k_27Al[7] = {0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06}; static const double a1c1_k_27Al[7] = {0.06, 0.06, 0.06, 0.06, 0.06, 0.06, 0.06};
//low band for 17F data
static const double a1c1_cfmin2_17F[7] = {0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10};
static const double a1c1_k2_17F[7] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
static const double a1c1_cfmin2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // no low band
static const double a1c1_k2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
double a1c1_cfmin2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
double a1c1_k2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
// active per-cell set, populated by dataset in Begin() // active per-cell set, populated by dataset in Begin()
double a1c1_cfmin_cell[7] = {0.20, 0.20, 0.20, 0.20, 0.20, 0.20, 0.20}; double a1c1_cfmin_cell[7] = {0.20, 0.20, 0.20, 0.20, 0.20, 0.20, 0.20};
double a1c1_k_cell[7] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}; double a1c1_k_cell[7] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
@ -242,14 +252,6 @@ inline double evalEloss(TSpline3 *fwd, TSpline3 *inv, double E, double pathlen)
return std::isfinite(e) ? e : 0.0; return std::isfinite(e) ? e : 0.0;
} }
static const double a1c1_cfmin2_17F[7] = {0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10};
static const double a1c1_k2_17F[7] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
static const double a1c1_cfmin2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // no low band
static const double a1c1_k2_27Al[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
double a1c1_cfmin2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
double a1c1_k2_cell[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
// a1c1_zcorr / A1C1CellSol / A1C1Sol / solve_cell / a1c1_solve / SideChoice / // a1c1_zcorr / A1C1CellSol / A1C1Sol / solve_cell / a1c1_solve / SideChoice /
// a1c1_pick_side / a1c1_solve_pick / a1c1_cfrac_pcz / a1c0_wirePos / // a1c1_pick_side / a1c1_solve_pick / a1c1_cfrac_pcz / a1c0_wirePos /
// a1c0_hybrid_pcz / a1c2_zfix now live in Armory/PCZRecon.h (included above), // a1c0_hybrid_pcz / a1c2_zfix now live in Armory/PCZRecon.h (included above),
@ -395,6 +397,7 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events, void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events,
const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters); const std::vector<std::vector<std::tuple<int, double, double>>> &aClusters, const std::vector<std::vector<std::tuple<int, double, double>>> &cClusters);
void a1c1CalibDiagnostic(HistPlotter *plotter, const std::vector<Event> &PC_Events); void a1c1CalibDiagnostic(HistPlotter *plotter, const std::vector<Event> &PC_Events);
void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events);
void TrackRecon::Begin(TTree * /*tree*/) void TrackRecon::Begin(TTree * /*tree*/)
{ {
@ -839,7 +842,7 @@ inline void pcEnergyCalibrationAccumulate(const std::vector<Event> &PC_Events,
{ {
if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) > phi_win) if (TMath::Abs(sievent.pos.DeltaPhi(pcevent.pos)) > phi_win)
return; return;
if (TMath::Abs(sievent.Time1 - pcevent.Time1) > 150) // time coincidence if (sievent.Time1 - pcevent.Time1 < 150) // time coincidence
return; return;
double theta = (sievent.pos - source_pos).Theta(); double theta = (sievent.pos - source_pos).Theta();
if (theta <= 0.0 || !std::isfinite(theta)) if (theta <= 0.0 || !std::isfinite(theta))
@ -1315,7 +1318,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (diagnostic_eplots) if (diagnostic_eplots)
{ {
// if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al // if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al
if (tRing - static_cast<double>(pc.t[k]) < -150) // proton tests, 27Al if (tRing - static_cast<double>(pc.t[k]) < 150) // proton tests, 27Al
{ {
PCAQQQTimeCut = true; PCAQQQTimeCut = true;
plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ"); plotter->Fill2D("CalibratedQQQEvsPCE_R", 1000, 0, 10, 2000, 0, 30000, eRingMeV, pc.e[k], "hPCQQQ");
@ -1554,20 +1557,21 @@ Bool_t TrackRecon::Process(Long64_t entry)
} }
} }
if ((pcEnergyCalibLoaded || doPCEnergyCalibration) && cClusters.empty()) if (cClusters.empty())
{ {
for (const auto &aCl : aClusters) for (const auto &aCl : aClusters)
{ {
if (aCl.size() != 1) // a1c0: exactly one anode wire -- same convention as if (aCl.size() < 1 || aCl.size() > 2) // A1C0 (1 wire) or A2C0 (2 wires) --
continue; // reaction_ax_core and miscHistograms_oneWire's a1c0 loops continue; // reaction_ax_core / miscHistograms_oneWire's
// a1c0 convention, one wire wider for A2C0.
if (clusterHasExcludedAnode(aCl)) 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);
double apSumE = std::get<1>(aPw); double apSumE = std::get<1>(aPw);
double apTSMaxE = std::get<3>(aPw); double apTSMaxE = std::get<3>(aPw);
int anodeIdx = std::get<0>(aCl[0]); int anodeIdx = std::get<0>(aCl[0]); // representative wire index (tag/sanity-check only,
if (anodeIdx < 0 || anodeIdx >= 24) if (anodeIdx < 0 || anodeIdx >= 24) // not assumed to be "the" wire for A2C0's 2-wire cluster)
continue; continue;
const Event *bestSi = nullptr; const Event *bestSi = nullptr;
@ -1594,27 +1598,34 @@ Bool_t TrackRecon::Process(Long64_t entry)
if (!bestSi) if (!bestSi)
continue; continue;
TVector3 pc = a1c0_wirePos(apwire, bestSi->pos.Phi(), bestIsQQQ); // same A1C0 z reference as the benchmark bool isA2C0 = (aCl.size() == 2);
TVector3 pc = isA2C0 ? a2c0_wirePos(apwire, bestSi->pos.Phi(), bestIsQQQ)
: a1c0_wirePos(apwire, bestSi->pos.Phi(), bestIsQQQ); // same z reference as the benchmark
Event PCEventRaw(pc, apSumE, -1.0, apTSMaxE, -1.0);
PCEventRaw.multi1 = static_cast<int>(aCl.size());
PCEventRaw.multi2 = 0;
PCEventRaw.Anodech = anodeIdx;
PCEventRaw.Cathodech = -1;
PC_Events.push_back(PCEventRaw);
if (pcEnergyCalibLoaded) if (pcEnergyCalibLoaded)
{ {
// aCl is guaranteed size 1 by the filter above, so anodeIdx unambiguously double anodeCalibSum = 0.0;
// identifies the single wire this event's calibration applies to. for (const auto &w : aCl)
double anodeCalibSum = (anodeIdx >= 0 && anodeIdx < 24) {
? pcEnergySlope[anodeIdx] * std::get<1>(aCl[0]) int wi = std::get<0>(w);
: 0.0; if (wi >= 0 && wi < 24)
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 = 1; ev.multi1 = static_cast<int>(aCl.size());
ev.multi2 = 0; // no cathode -> a1c0 topology in pcCalibratedHistograms ev.multi2 = 0; // no cathode -> a1c0/a2c0 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);
} }
// 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.
} }
} }
@ -1734,6 +1745,7 @@ Bool_t TrackRecon::Process(Long64_t entry)
pcCalibratedHistograms(plotter, QQQ_Events, SX3_Events, PC_Events_calibrated); pcCalibratedHistograms(plotter, QQQ_Events, SX3_Events, PC_Events_calibrated);
a1c1CalibDiagnostic(plotter, PC_Events); // <-- new, unconditional a1c1CalibDiagnostic(plotter, PC_Events); // <-- new, unconditional
pcVertexByWireGeometry(plotter, QQQ_Events, SX3_Events, PC_Events); // <-- new, unconditional
auto hasPCCoincidence = [&](const TVector3 &pos) auto hasPCCoincidence = [&](const TVector3 &pos)
{ {
@ -1988,6 +2000,71 @@ void a1c1CalibDiagnostic(HistPlotter *plotter, const std::vector<Event> &PC_Even
} }
} }
void pcVertexByWireGeometry(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events)
{
static TRandom3 rand(0); // seeded once, not per call -- dithers A1C0's Z below
auto fillFor = [&](const std::vector<Event> &sis, bool isQQQ)
{
double phi_win = isQQQ ? TMath::Pi() / 4.0 : TMath::Pi() / 3.0; // same per-detector
double perp_max = isQQQ ? 6.0 : 10.0; // tolerances used
const std::string det = isQQQ ? "_QQQ" : "_SX3"; // elsewhere in this file
for (const auto &pcevent : PC_Events)
{
// Only topologies with an established pcz method below -- A2C1/A2C2 etc.
// don't have one yet, so they're skipped here rather than silently
// falling back to a raw, un-dispatched pos.Z().
bool knownTopo = (pcevent.multi1 == 1 && pcevent.multi2 == 2) ||
(pcevent.multi1 == 1 && pcevent.multi2 == 1) ||
(pcevent.multi1 == 1 && pcevent.multi2 == 0) ||
(pcevent.multi1 == 2 && pcevent.multi2 == 0);
if (!knownTopo)
continue;
for (const auto &si : sis)
{
if (TMath::Abs(si.pos.DeltaPhi(pcevent.pos)) > phi_win)
continue;
if (si.Time1 - pcevent.Time1 < 150) // loose time coincidence, same convention as elsewhere
continue;
double pcz;
bool a1c1_inband = false;
if (pcevent.multi1 == 1 && pcevent.multi2 == 2) // A1C2
pcz = a1c2_zfix(pcevent.pos.Z());
else if (pcevent.multi1 == 1 && pcevent.multi2 == 1) // A1C1
pcz = a1c1_cfrac_pcz(pcevent, si.pos, a1c1_inband);
else if (pcevent.multi1 == 1 && pcevent.multi2 == 0) // A1C0
pcz = rand.Gaus(pcevent.pos.Z(), dither_sigma_c0 / 2.0);
else // A2C0 (multi1==2, multi2==0) -- undithered by design
pcz = pcevent.pos.Z();
TVector3 x2(pcevent.pos.X(), pcevent.pos.Y(), pcz);
TVector3 vtx = beamVertex(si.pos, x2 - si.pos);
if (beamPerp(vtx) > perp_max)
continue;
if (vtx.Z() < -173.6 || vtx.Z() > 100)
continue;
std::string topo = "_a" + std::to_string(pcevent.multi1) + "c" + std::to_string(pcevent.multi2);
plotter->Fill2D("WireGeometry_dE_vs_VertexZ" + topo, 800, -400, 400, 800, 0, 40000, vtx.Z(), pcevent.Energy1, "WireGeometry");
plotter->Fill2D("WireGeometry_dE_vs_VertexZ" + topo + det, 800, -400, 400, 800, 0, 40000, vtx.Z(), pcevent.Energy1, "WireGeometry");
if (pcevent.multi1 == 1 && pcevent.multi2 == 1 && a1c1_inband)
{
plotter->Fill2D("WireGeometry_dE_vs_VertexZ_a1c1_inband", 800, -400, 400, 800, 0, 40000, vtx.Z(), pcevent.Energy1, "WireGeometry");
plotter->Fill2D("WireGeometry_dE_vs_VertexZ_a1c1_inband" + det, 800, -400, 400, 800, 0, 40000, vtx.Z(), pcevent.Energy1, "WireGeometry");
}
}
}
};
fillFor(QQQ_Events, true);
fillFor(SX3_Events, false);
}
void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events_calibrated) void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Events, const std::vector<Event> &SX3_Events, const std::vector<Event> &PC_Events_calibrated)
{ {
static TRandom3 rand(0); // seeded once, not per call -- for Si-side pixel/strip dithering below static TRandom3 rand(0); // seeded once, not per call -- for Si-side pixel/strip dithering below
@ -2066,7 +2143,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
for (const auto &qqqevent : QQQ_Events) for (const auto &qqqevent : QQQ_Events)
{ {
bool phicut = TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0; bool phicut = TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
bool timecut = TMath::Abs(qqqevent.Time1 - pcevent.Time1) < 150; bool timecut = (qqqevent.Time1 - pcevent.Time1) < 150;
if (!(phicut && timecut)) if (!(phicut && timecut))
continue; continue;
@ -2081,16 +2158,16 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
double Egu_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, qqqevent.Energy1, pcc.guard_cm); double Egu_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, qqqevent.Energy1, pcc.guard_cm);
double Eca_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, qqqevent.Energy1, pcc.cathode_cm); double Eca_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, qqqevent.Energy1, pcc.cathode_cm);
plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asProton" + topo, 800, 0, 2, 400, 0, 0.6, pcevent.Energy1, Egu_p - Eca_p, "hCalibPC"); plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asProton" + topo, 400, 0, 0.6, 400, 0, 0.6, pcevent.Energy1, Egu_p - Eca_p, "hCalibPC");
double Egu_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, qqqevent.Energy1, pcc.guard_cm); double Egu_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, qqqevent.Energy1, pcc.guard_cm);
double Eca_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, qqqevent.Energy1, pcc.cathode_cm); double Eca_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, qqqevent.Energy1, pcc.cathode_cm);
plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asAlpha" + topo, 800, 0, 2, 400, 0, 0.6, pcevent.Energy1, Egu_a - Eca_a, "hCalibPC"); plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asAlpha" + topo, 400, 0, 0.6, 400, 0, 0.6, pcevent.Energy1, Egu_a - Eca_a, "hCalibPC");
} }
for (const auto &sx3event : SX3_Events) for (const auto &sx3event : SX3_Events)
{ {
bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0; bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 4.0;
bool timecut = TMath::Abs(sx3event.Time1 - pcevent.Time1) < 150; bool timecut = (sx3event.Time1 - pcevent.Time1) < 150;
if (!(phicut && timecut)) if (!(phicut && timecut))
continue; continue;
@ -2106,11 +2183,11 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_
double Egu_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.guard_cm); double Egu_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.guard_cm);
double Eca_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.cathode_cm); double Eca_p = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.cathode_cm);
plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asProton" + topo, 800, 0, 2, 400, 0, 0.6, pcevent.Energy1, Egu_p - Eca_p, "hCalibPC"); plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asProton" + topo, 400, 0, 0.6, 400, 0, 0.6, pcevent.Energy1, Egu_p - Eca_p, "hCalibPC");
double Egu_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, sx3event.Energy1, pcc.guard_cm); double Egu_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, sx3event.Energy1, pcc.guard_cm);
double Eca_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, sx3event.Energy1, pcc.cathode_cm); double Eca_a = evalEloss(MeV_to_cm_spl, cm_to_MeV_spl, sx3event.Energy1, pcc.cathode_cm);
plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asAlpha" + topo, 800, 0, 2, 400, 0, 0.6, pcevent.Energy1, Egu_a - Eca_a, "hCalibPC"); plotter->Fill2D("Calib_dEgasPred_vs_dEgasCalib_asAlpha" + topo, 400, 0, 0.6, 400, 0, 0.6, pcevent.Energy1, Egu_a - Eca_a, "hCalibPC");
} }
} }
} }
@ -2196,9 +2273,13 @@ void PCSX3ClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
plotter->Fill2D("dE_E_Anodesx3B_a1c0", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E"); plotter->Fill2D("dE_E_Anodesx3B_a1c0", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 1 && pcevent.multi2 == 0) if (pcevent.multi1 == 1 && pcevent.multi2 == 0)
plotter->Fill2D("dE_E_Cathodesx3B_a1c0", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E"); plotter->Fill2D("dE_E_Cathodesx3B_a1c0", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 0)
plotter->Fill2D("dE_E_Anodesx3B_a2c0", 400, 0, 30, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, "PID_dE_E");
if (pcevent.multi1 == 2 && pcevent.multi2 == 0)
plotter->Fill2D("dE_E_Cathodesx3B_a2c0", 400, 0, 30, 800, 0, 10000, sx3event.Energy1, pcevent.Energy2, "PID_dE_E");
plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < -150), 100, -200, 200, 100, -200, 200, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles"); plotter->Fill2D("sx3phi_vs_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < 150), 100, -200, 200, 100, -200, 200, sx3event.pos.Phi() * 180 / M_PI, pcevent.pos.Phi() * 180 / M_PI, "Kinematics_Angles");
plotter->Fill1D("sx3phi_minus_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < -150), 100, -180, 180, (sx3event.pos.DeltaPhi(pcevent.pos)) * 180 / M_PI, "Kinematics_Angles"); plotter->Fill1D("sx3phi_minus_pcphi" + std::to_string(sx3event.Time1 - pcevent.Time1 < 150), 100, -180, 180, (sx3event.pos.DeltaPhi(pcevent.pos)) * 180 / M_PI, "Kinematics_Angles");
if (PCSX3TimeCut) if (PCSX3TimeCut)
{ {
@ -2716,6 +2797,11 @@ void PCQQQClusterAnalysis(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
plotter->Fill2D("dE_E_Anodesx3B_a1c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E"); plotter->Fill2D("dE_E_Anodesx3B_a1c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
plotter->Fill2D("dE_E_Cathodesx3B_a1c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E"); plotter->Fill2D("dE_E_Cathodesx3B_a1c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
} }
if (pcevent.multi1 == 2 && pcevent.multi2 == 0)
{
plotter->Fill2D("dE_E_Anodesx3B_a2c0", 400, 0, 30, 800, 0, 40000, qqqevent.Energy1, pcevent.Energy1, "PID_dE_E");
plotter->Fill2D("dE_E_Cathodesx3B_a2c0", 400, 0, 30, 800, 0, 10000, qqqevent.Energy1, pcevent.Energy2, "PID_dE_E");
}
if (phicut) if (phicut)
{ {
plotter->Fill2D("dE2_E_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta, "PID_dE_E"); plotter->Fill2D("dE2_E_AnodeQQQR_TC1PC1_pidlow" + std::to_string(lowercut_cath), 400, 0, 30, 800, 0, 4000, qqqevent.Energy1, pcevent.Energy1 * sinTheta, "PID_dE_E");
@ -3389,7 +3475,13 @@ void protonMiscHistograms(HistPlotter *plotter, const std::vector<Event> &QQQ_Ev
// if(qqqevent.Energy1 > 5.0) continue; //coarse gating // if(qqqevent.Energy1 > 5.0) continue; //coarse gating
for (const auto &pcevent : PC_Events) for (const auto &pcevent : PC_Events)
{ {
if (!(pcevent.multi1 == 1 && pcevent.multi2 <= 2)) // A1C0/A1C1/A1C2 (multi1==1, multi2 in {0,1,2}) plus A2C0 (multi1==2,
// multi2==0) -- the only no-cathode topology besides A1C0. multi1==2
// otherwise means A2C1/A2C2 (two-wire anode cluster WITH a cathode),
// which is intentionally still excluded here, same as before.
bool topoOK = (pcevent.multi1 == 1 && pcevent.multi2 <= 2) ||
(pcevent.multi1 == 2 && pcevent.multi2 == 0);
if (!topoOK)
continue; continue;
// if(pcevent.Energy1 > 11000) continue; //coarse gating // if(pcevent.Energy1 > 11000) continue; //coarse gating
@ -3767,6 +3859,65 @@ void protonMiscHistograms_sx3(HistPlotter *plotter, const std::vector<Event> &QQ
} }
} }
} // end A1C1 comparison loop } // end A1C1 comparison loop
for (const auto &pcevent : PC_Events)
{
bool topoOK = (pcevent.multi1 == 1 && pcevent.multi2 == 0) || // A1C0
(pcevent.multi1 == 2 && pcevent.multi2 == 0); // A2C0
if (!topoOK)
continue;
bool phicut = TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) <= TMath::Pi() / 3.0;
if (!phicut)
continue;
TVector3 x1(sx3event.pos);
TVector3 r_rhoMin = beamVertex(x1, pcevent.pos - x1); // no z-fix needed -- A1C0/A2C0's
double vertex_z = r_rhoMin.Z(); // pos.Z() is already the true wire z
if (beamPerp(r_rhoMin) > 10.0)
continue;
if (vertex_z < -173.6 || vertex_z > 100)
continue; // same beam-region acceptance as the A1C2/A1C1 loops above
double theta_s = (sx3event.pos - r_rhoMin).Theta();
double sinTheta_customV = TMath::Sin(theta_s);
double path_length_s = pathLengthCm(sx3event.pos, r_rhoMin);
// No cathode signal to pick an ejectile hypothesis from (there's no
// Energy2 to test against the 1400 threshold) -- proton table only,
// the same default the A1C2/A1C1 loops fall back to for their
// "_cathode_protons" tag.
double sx3Efix = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, path_length_s);
std::string tag = "_a" + std::to_string(pcevent.multi1) + "c0";
std::string pmlabel = "proton+miscsx3" + tag;
plotter->Fill2D("pmiscs_dE_E_Anodesx3" + tag, 400, 0, 10, 800, 0, 40000, sx3event.Energy1, pcevent.Energy1, pmlabel);
plotter->Fill2D("pmiscs_dE3_E_Anodesx3" + tag, 400, 0, 10, 400, 0, 40000, sx3event.Energy1, pcevent.Energy1 * sinTheta_customV * 3., pmlabel);
plotter->Fill1D("pmiscs_pcz" + tag, 600, -300, 300, pcevent.pos.Z(), pmlabel);
plotter->Fill2D("pmiscs_dE3_Ef_Anodesx3" + tag, 400, 0, 10, 400, 0, 40000, sx3Efix, pcevent.Energy1 * sinTheta_customV * 3, pmlabel);
plotter->Fill2D("pmiscs_Ef_vs_theta_sx3" + tag, 100, 0, 180, 800, 0, 20, theta_s * 180 / M_PI, sx3Efix, pmlabel);
plotter->Fill1D("pmiscs_VertexReconZ" + tag, 800, -400, 400, vertex_z, pmlabel);
plotter->Fill2D("pmiscs_VertexReconXY" + tag, 200, -100, 100, 200, -100, 100, r_rhoMin.X(), r_rhoMin.Y(), pmlabel);
plotter->Fill2D("pmiscs_VertexReconZ_vs_Ef" + tag, 800, -400, 400, 800, 0, 20, vertex_z, sx3Efix, pmlabel);
// Gas segmentation validation, mirroring the A1C2/A1C1 loops' dEgas family.
PCCollect pcc = pcCollectionPath(r_rhoMin, sx3event.pos);
if (pcc.ok)
{
double E_gu = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.guard_cm);
double E_ca = evalEloss(MeV_to_cm_p_spl, cm_to_MeVp_spl, sx3event.Energy1, pcc.cathode_cm);
double dE_pred = E_gu - E_ca;
plotter->Fill2D("pmiscs_dEgas_vs_Ef" + tag, 400, 0, 20, 400, 0, 0.6, sx3Efix, dE_pred, pmlabel);
if (pcevent.Anodech >= 0 && pcevent.Anodech < 24)
{
double anodeE_MeV = pcEnergySlope[pcevent.Anodech] * pcevent.Energy1;
plotter->Fill2D("pmiscs_dEgasCalib_vs_Ef" + tag, 400, 0, 20, 800, 0, 0.6, sx3Efix, anodeE_MeV, pmlabel);
plotter->Fill2D("pmiscs_dEgasCalib_vs_VertexZ" + tag, 800, -400, 400, 800, 0, 0.6, vertex_z, anodeE_MeV, pmlabel);
plotter->Fill2D("pmiscs_dEgasPred_vs_dEgasCalib" + tag, 800, 0, 2, 400, 0, 0.6, anodeE_MeV, dE_pred, pmlabel);
}
}
} // end A1C0/A2C0 loop
} // end sx3Events loop } // end sx3Events loop
} }
@ -3827,12 +3978,7 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
const AAEjectileMasses &ej_m, const std::string &globaltag) const AAEjectileMasses &ej_m, const std::string &globaltag)
{ {
const std::string sfx = "_" + det + globaltag; const std::string sfx = "_" + det + globaltag;
static TRandom3 rand(0); // seeded once (random seed via TUUID), not per call -- static TRandom3 rand(0);
// used only to dither a1c0's Z below, matching dither_sigma_c0's
// established use elsewhere in this file (e.g. a1c0_hybrid_pcz).
// a1c1 is deliberately left undithered: its cfrac-based sub-wire-pitch
// fraction already gives continuous Z, unlike a1c0's single-wire position.
for (const auto &sievent : Si_Events) for (const auto &sievent : Si_Events)
{ {
if (sievent.Energy1 < si_ecut) if (sievent.Energy1 < si_ecut)
@ -4032,46 +4178,51 @@ static void reaction_ax_core(HistPlotter *plotter, const std::vector<Event> &Si_
for (const auto &aCl : aClusters) for (const auto &aCl : aClusters)
{ {
if (aCl.size() != 1) // a1c0: exactly one anode wire, no cathode -- same if (aCl.size() < 1 || aCl.size() > 2)
continue; // convention as pcevent.multi1==1 && multi2==0 elsewhere continue;
if (clusterHasExcludedAnode(aCl)) 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);
double apSumE = std::get<1>(aPw); double apSumE = std::get<1>(aPw);
TVector3 pc = a1c0_wirePos(apwire, sievent.pos.Phi(), isQQQ); bool isA2C0 = (aCl.size() == 2);
const std::string a0tag = isA2C0 ? "a2c0" : "a1c0";
TVector3 pc = isA2C0 ? a2c0_wirePos(apwire, sievent.pos.Phi(), isQQQ)
: a1c0_wirePos(apwire, sievent.pos.Phi(), isQQQ);
if (TMath::Abs(sievent.pos.DeltaPhi(pc)) > phi_win) if (TMath::Abs(sievent.pos.DeltaPhi(pc)) > phi_win)
continue; continue;
std::string pmlabel = globaltag + "_" + rx + "+misc_" + det + "_a1c0"; std::string pmlabel = globaltag + "_" + rx + "+misc_" + det + "_" + a0tag;
plotter->Fill2D(rx + "_dE_E_Anode_a1c0" + sfx, 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, apSumE, pmlabel); plotter->Fill2D(rx + "_dE_E_Anode_" + a0tag + sfx, 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, apSumE, pmlabel);
TVector3 r_rhoMin_a1c0 = beamVertex(sievent.pos, pc - sievent.pos); TVector3 r_rhoMin_a0 = beamVertex(sievent.pos, pc - sievent.pos);
double beam_path_length_a1c0 = TMath::Abs(r_rhoMin_a1c0.Z() - z_entrance) * 0.1; double beam_path_length_a0 = TMath::Abs(r_rhoMin_a0.Z() - z_entrance) * 0.1;
double beam_energy_at_vertex_a1c0 = evalElossForward(beam_MeV_to_cm, beam_cm_to_MeV, beamE0, beam_path_length_a1c0); double beam_energy_at_vertex_a0 = evalElossForward(beam_MeV_to_cm, beam_cm_to_MeV, beamE0, beam_path_length_a0);
plotter->Fill2D(rx + "_dE_E_Anode_a1c0" + sfx + "_10MeV" + std::to_string(beam_energy_at_vertex_a1c0 < 10), 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, apSumE, pmlabel); plotter->Fill2D(rx + "_dE_E_Anode_" + a0tag + sfx + "_10MeV" + std::to_string(beam_energy_at_vertex_a0 < 10), 400, 0, dEa_max, 800, 0, 40000, sievent.Energy1, apSumE, pmlabel);
plotter->Fill2D(rx + "_dPhi_a1c0" + sfx, 100, -200, 200, 100, -200, 200, pc.Phi() * 180 / M_PI, sievent.pos.Phi() * 180 / M_PI, pmlabel); plotter->Fill2D(rx + "_dPhi_" + a0tag + sfx, 100, -200, 200, 100, -200, 200, pc.Phi() * 180 / M_PI, sievent.pos.Phi() * 180 / M_PI, pmlabel);
plotter->Fill1D(rx + "_rawZ_a1c0" + sfx, 600, -300, 300, pc.Z(), pmlabel); plotter->Fill1D(rx + "_rawZ_" + a0tag + sfx, 600, -300, 300, pc.Z(), pmlabel);
// Calibrated anode energy for the a1c0 wire, using the same pcEnergySlope int anodeCh_a0 = std::get<0>(aCl[0]);
// calibration already applied to A1C0 events elsewhere in this file (see the double anodeE_MeV_a0 = 0.0;
// pcEnergyCalibLoaded block above). aCl is guaranteed size 1 by the filter bool anyValidWire = false;
// above, so aCl[0] is unambiguously "the" wire for this event. for (const auto &w : aCl)
int anodeCh_a1c0 = std::get<0>(aCl[0]); {
double anodeE_MeV_a1c0 = (anodeCh_a1c0 >= 0 && anodeCh_a1c0 < 24) int wi = std::get<0>(w);
? pcEnergySlope[anodeCh_a1c0] * std::get<1>(aCl[0]) if (wi >= 0 && wi < 24)
: -1.0; {
if (anodeCh_a1c0 < 0 || anodeCh_a1c0 >= 24) anodeE_MeV_a0 += pcEnergySlope[wi] * std::get<1>(w);
anodeCh_a1c0 = -1; anyValidWire = true;
}
}
if (!anyValidWire)
anodeE_MeV_a0 = -1.0;
if (anodeCh_a0 < 0 || anodeCh_a0 >= 24)
anodeCh_a0 = -1;
// a1c0 Z is a deterministic function of wire position (a1c1_zcorr is just a double pcz_a0 = isA2C0 ? pc.Z() : rand.Gaus(pc.Z(), dither_sigma_c0 / 2.0);
// scale+offset) with no sub-wire-pitch information, unlike a1c1's cfrac -- so reconstructAndFill(pcz_a0, pc, apSumE, -1.0, anodeE_MeV_a0, -1.0, a0tag, "", anodeCh_a0);
// dither only the Z fed into reconstruction, matching dither_sigma_c0's use
// elsewhere (e.g. a1c0_hybrid_pcz). pc itself stays raw/undithered: _rawZ_a1c0,
// the phi cut, and _dPhi_a1c0 above are all meant to reflect the true wire position.
double pcz_a1c0_dith = rand.Gaus(pc.Z(), dither_sigma_c0 / 2.0);
reconstructAndFill(pcz_a1c0_dith, pc, apSumE, -1.0, anodeE_MeV_a1c0, -1.0, "a1c0", "", anodeCh_a1c0);
} }
} }
} }

5
eloss_calculations/Eloss.py
View File

@ -6,7 +6,7 @@ P_TORR = 250
TEMP_K = 293.15 TEMP_K = 293.15
R = 8.3144 R = 8.3144
MEV2U = 1.0 / 931.494 MEV2U = 1.0 / 931.494
P_CO2 = 4 P_CO2 = 3
# Gas Density Calculations # Gas Density Calculations
p_pa = P_TORR * 133.322 p_pa = P_TORR * 133.322
@ -77,6 +77,9 @@ def generate_lookup(z, mass_u, e_start_mev, label):
projectile.T(e_u) projectile.T(e_u)
# dedx returns MeV / (g/cm2) # dedx returns MeV / (g/cm2)
if(mass_u >=10.0):
loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2 * 0.89
else:
loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2 loss_mev = catima.dedx(projectile, gas_mix) * step_g_cm2
current_e_total = max(0.0, current_e_total - loss_mev) current_e_total = max(0.0, current_e_total - loss_mev)

96
pc_energy_calibration.dat
View File

@ -1,48 +1,48 @@
0 2.270700e-05 0.000000e+00 1 0 7.569000E-06 0.000000E+00 1
1 3.849247e-05 0.000000e+00 1 1 1.283082E-05 0.000000E+00 1
2 3.849247e-05 0.000000e+00 1 2 1.283082E-05 0.000000E+00 1
3 3.849247e-05 0.000000e+00 1 3 1.283082E-05 0.000000E+00 1
4 3.849247e-05 0.000000e+00 1 4 1.283082E-05 0.000000E+00 1
5 4.489228e-05 0.000000e+00 2 5 1.496409E-05 0.000000E+00 2
6 4.489228e-05 0.000000e+00 1 6 2.036882E-05 0.000000E+00 1
7 3.374396e-05 0.000000e+00 1 7 1.124799E-05 0.000000E+00 1
8 3.370297e-05 0.000000e+00 1 8 1.123432E-05 0.000000E+00 1
9 1.000000e+00 0.000000e+00 3 9 3.333333E-01 0.000000E+00 3
10 3.876003e-05 0.000000e+00 1 10 1.292001E-05 0.000000E+00 1
11 2.636256e-05 0.000000e+00 1 11 8.787520E-06 0.000000E+00 1
12 1.000000e+00 0.000000e+00 3 12 3.333333E-01 0.000000E+00 3
13 2.360190e-05 0.000000e+00 1 13 7.867300E-06 0.000000E+00 1
14 2.012292e-05 0.000000e+00 1 14 6.707640E-06 0.000000E+00 1
15 2.554712e-05 0.000000e+00 1 15 8.515706E-06 0.000000E+00 1
16 3.124900e-05 0.000000e+00 1 16 1.041633E-05 0.000000E+00 1
17 3.124900e-05 0.000000e+00 1 17 1.041633E-05 0.000000E+00 1
18 3.849247e-05 0.000000e+00 1 18 1.283082E-05 0.000000E+00 1
19 2.712940e-05 0.000000e+00 2 19 9.043133E-06 0.000000E+00 2
20 3.140230e-05 0.000000e+00 1 20 1.046743E-05 0.000000E+00 1
21 3.140230e-05 0.000000e+00 2 21 1.046743E-05 0.000000E+00 2
22 3.142680e-05 0.000000e+00 2 22 1.047560E-05 0.000000E+00 2
23 3.849247e-05 0.000000e+00 2 23 1.283082E-05 0.000000E+00 2
24 3.024648e-05 0.000000e+00 1 24 1.008216E-05 0.000000E+00 1
25 3.732136e-05 0.000000e+00 1 25 1.244045E-05 0.000000E+00 1
26 3.486437e-05 0.000000e+00 1 26 1.162146E-05 0.000000E+00 1
27 5.511389e-05 0.000000e+00 1 27 1.837130E-05 0.000000E+00 1
28 4.472116e-05 0.000000e+00 1 28 1.490705E-05 0.000000E+00 1
29 3.655778e-05 0.000000e+00 1 29 1.218593E-05 0.000000E+00 1
30 3.579373e-05 0.000000e+00 1 30 1.193124E-05 0.000000E+00 1
31 3.507960e-05 0.000000e+00 1 31 1.169320E-05 0.000000E+00 1
32 3.613826e-05 0.000000e+00 1 32 1.204609E-05 0.000000E+00 1
33 2.886743e-05 0.000000e+00 1 33 9.622476E-06 0.000000E+00 1
34 3.031249e-05 0.000000e+00 1 34 1.010416E-05 0.000000E+00 1
35 2.983830e-05 0.000000e+00 1 35 9.946100E-06 0.000000E+00 1
36 3.735619e-05 0.000000e+00 1 36 1.245206E-05 0.000000E+00 1
37 6.012566e-05 0.000000e+00 2 37 2.004189E-05 0.000000E+00 2
38 3.303600e-05 0.000000e+00 1 38 1.101200E-05 0.000000E+00 1
39 6.039465e-05 0.000000e+00 1 39 2.013155E-05 0.000000E+00 1
40 4.511714e-05 0.000000e+00 1 40 1.503905E-05 0.000000E+00 1
41 5.086580e-05 0.000000e+00 1 41 1.695527E-05 0.000000E+00 1
42 3.765319e-05 0.000000e+00 1 42 1.255106E-05 0.000000E+00 1
43 5.656884e-05 0.000000e+00 1 43 1.885628E-05 0.000000E+00 1
44 1.196229e-04 0.000000e+00 2 44 3.987430E-05 0.000000E+00 2
45 2.972753e-05 0.000000e+00 1 45 9.909176E-06 0.000000E+00 1
46 2.829864e-05 0.000000e+00 1 46 9.432880E-06 0.000000E+00 1
47 3.309860e-05 0.000000e+00 1 47 1.103287E-05 0.000000E+00 1