From 3e9c86d8dc7b9caaf4e7c033387727528e871970 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Tue, 28 Jul 2026 11:32:30 -0400 Subject: [PATCH] modified: .gitignore modified: TrackRecon.C new file: pc_energy_calibration.dat modified: pccal/fit_pc_energy_calibration.C modified: run_27Al.sh modified: run_tr.sh --- .gitignore | 1 + TrackRecon.C | 283 ++++++++++++++++++++---------- pc_energy_calibration.dat | 48 +++++ pccal/fit_pc_energy_calibration.C | 194 ++++++++++++++++---- run_27Al.sh | 19 +- run_tr.sh | 10 +- 6 files changed, 413 insertions(+), 142 deletions(-) create mode 100644 pc_energy_calibration.dat diff --git a/.gitignore b/.gitignore index 66b5d0f..fdcad5f 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,4 @@ Armory/EventBuilder EventBuilder anasen_fem/anode_to_cathode_1d_17.43.csv CLAUDE.md +pc_calib_raw/ \ No newline at end of file diff --git a/TrackRecon.C b/TrackRecon.C index d3f1788..0b875da 100644 --- a/TrackRecon.C +++ b/TrackRecon.C @@ -44,15 +44,16 @@ bool process_alpha_proton_scattering = false, doMiscHistograms = true, doPCSX3ClusterAnalysis = true, doPCQQQClusterAnalysis = true, - doOldAnalysis = false, + doOldAnalysis = true, do27AlapAnalysis = false, BenchMark = true, onwire_analysis = true, diagnostic_eplots = true, - diagnostic_tplots = false, + diagnostic_tplots = true, reactiondata = false, doPCEnergyCalibration = false, - ta_foil_run = false; + ta_foil_run = false, + source_run = false; // --- Geometry, Calibration, & Model Variables --- double source_vertex = 53.0, @@ -153,7 +154,17 @@ static std::vector a1c1_dead_cathode_17F = {}; // 0,13,15 can be recover static std::vector a1c1_dead_anode_27Al = {0, 12, 19}; static std::vector a1c1_dead_cathode_27Al = {13}; -std::vector> pcCalibData[48]; +// Calibration points are streamed straight to disk as they're generated +// (opened in Begin(), written by pcCalibWritePoint(), closed in Terminate()) +// instead of buffered in memory for the whole run -- a full run can produce +// many millions of points, and holding them all in a std::vector until +// Terminate() grows unbounded and can exhaust memory on long runs. +std::ofstream pcCalibOutFile; +inline void pcCalibWritePoint(int wire, double adc, double dE_gas) +{ + if (pcCalibOutFile.is_open()) + pcCalibOutFile << wire << " " << adc << " " << dE_gas << "\n"; +} std::vector *a1c1_dead_anode = &a1c1_dead_anode_17F; // active set, chosen in Begin() std::vector *a1c1_dead_cathode = &a1c1_dead_cathode_17F; @@ -507,20 +518,27 @@ void TrackRecon::Begin(TTree * /*tree*/) std::cout << "Analyzing dataset as reactiondata" << std::endl; } - if (getenv("RUN_NUMBER")) + // Run classification is by OUT_DIR, not RUN_NUMBER: run numbers collide + // across datasets/blocks (e.g. 17F's alpha+gas block and 27Al's proton + // block both use runs 18-21), so a RUN_NUMBER-keyed lookup misclassifies + // whichever dataset didn't originally populate kTaFoilRuns. + std::string outdir = getenv("OUT_DIR") ? getenv("OUT_DIR") : ""; + ta_foil_run = (outdir == "Output_p"); + source_run = (outdir == "Output_a"); + if (ta_foil_run && getenv("RUN_NUMBER")) { int run_number = std::atoi(getenv("RUN_NUMBER")); for (const auto &r : kTaFoilRuns) { if (r.run == run_number) { - ta_foil_run = true; ta_foil_z_mm = r.z_mm; break; } } } - std::cout << "Ta foil: " << (ta_foil_run ? ("present, run in proton-scattering campaign, z=" + std::to_string(ta_foil_z_mm) + " mm") : std::string("not applicable (RUN_NUMBER unset or not a proton-scattering run)")) << std::endl; + std::cout << "OUT_DIR=" << outdir << " -> ta_foil_run=" << ta_foil_run + << " (z=" << ta_foil_z_mm << " mm), source_run=" << source_run << std::endl; // if (getenv("PC_ENERGY_CALIBRATION")) // doPCEnergyCalibration = std::atoi(getenv("PC_ENERGY_CALIBRATION")) != 0; @@ -530,6 +548,25 @@ void TrackRecon::Begin(TTree * /*tree*/) if (getenv("source_vertex")) source_vertex = (double)std::atof(std::string(getenv("source_vertex")).c_str()); + if (doPCEnergyCalibration) + { + gSystem->mkdir("pc_calib_raw", kTRUE); // kTRUE = create parents, no-op if it exists + // Tag source-run vs proton-run files distinctly (src_/ap_ prefix) so the + // aggregator can treat them differently: source runs sit at a handful of + // known, fixed positions -- each run is a single clean, additive + // calibration point per wire, not thousands of samples of the same thing -- + // while proton runs sample a continuously varying vertex/angle within a + // single run, where the per-event spread IS the useful signal and must + // stay pooled at the individual-point level. + 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())); + std::string outname = "pc_calib_raw/points_" + tag + ".dat"; + pcCalibOutFile.open(outname); + pcCalibOutFile << std::scientific << std::setprecision(6); + std::cout << "PC energy calibration: streaming raw points to " << outname << std::endl; + } + if (getenv("CO2percent")) co2pc = std::atoi(getenv("CO2percent")); std::cout << "CO2 percent set to " << co2pc << std::endl; @@ -787,7 +824,7 @@ inline double evalElossForward(TSpline3 *fwd, TSpline3 *inv, double E, double pa return e; } -inline void pcEnergyCalibrationAccumulate(const std::vector &PC_Events) +inline void pcEnergyCalibrationAccumulate(const std::vector &PC_Events, const std::vector &SX3_Events, const std::vector &QQQ_Events) { const TVector3 source_pos(beam_axis_x, beam_axis_y, source_vertex); for (const auto &pcevent : PC_Events) @@ -795,26 +832,66 @@ inline void pcEnergyCalibrationAccumulate(const std::vector &PC_Events) if (!(pcevent.multi1 >= 1 && pcevent.multi2 >= 1)) continue; - TVector3 interaction = pcevent.pos; - if (pcevent.multi1 == 1 && pcevent.multi2 == 1) + // Use an Si-hit-derived z (pczguess) instead of the PC's own cfrac + // sub-cell reconstruction: after depositing dE_gas in the PC, the alpha + // continues on and strikes SX3 or QQQ, giving a genuinely INDEPENDENT + // trajectory-angle measurement from the known, fixed source position -- + // same pczguess pattern used elsewhere in this file for A1C0 benchmarks + // (each detector keeps its own established formula/geometry constants + // for consistency with the rest of the file). This sidesteps the + // PC-only reconstruction entirely (no more circular dependence on the + // PC's own charge-division z, and no more dependence on anode/cathode + // multiplicity for z precision), so the push gates below only need to + // protect ADC purity, not z precision. Try SX3 first, then QQQ. + bool foundSi = false; + double pcz = 0.0; + for (const auto &sx3event : SX3_Events) { - bool inband = false; - double pcz = a1c1_cfrac_pcz(pcevent, source_pos, inband); - if (!inband) + if (!(std::isfinite(sx3event.Time1) && std::isfinite(pcevent.Time1)) || TMath::Abs(sx3event.Time1 - pcevent.Time1) > 150.0) continue; - interaction.SetZ(pcz); + if (TMath::Abs(sx3event.pos.DeltaPhi(pcevent.pos)) > TMath::Pi() / 3.0) + continue; + double sx3theta = TMath::ATan2(88.0, sx3event.pos.Z() - source_vertex); + pcz = 37.0 / TMath::Tan(sx3theta) + source_vertex; + foundSi = true; + break; } + if (!foundSi) + { + for (const auto &qqqevent : QQQ_Events) + { + if (!(std::isfinite(qqqevent.Time1) && std::isfinite(pcevent.Time1)) || TMath::Abs(qqqevent.Time1 - pcevent.Time1) > 150.0) + continue; + if (TMath::Abs(qqqevent.pos.DeltaPhi(pcevent.pos)) > TMath::Pi() / 4.0) + continue; + double qqqTheta = (qqqevent.pos - TVector3(0, 0, source_vertex)).Theta(); + pcz = 37.0 / TMath::Tan(qqqTheta) + source_vertex; + foundSi = true; + break; + } + } + if (!foundSi) + continue; + TVector3 interaction(pcevent.pos.X(), pcevent.pos.Y(), pcz); double path_length = pathLengthCm(source_pos, interaction); double e_remaining = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, pc_calib_alpha_source_mev, path_length); double dE_gas = pc_calib_alpha_source_mev - e_remaining; if (!std::isfinite(dE_gas) || dE_gas <= 0.0) continue; - if (pcevent.Anodech >= 0 && pcevent.Anodech < 24) - pcCalibData[pcevent.Anodech].push_back({pcevent.Energy1, dE_gas}); - if (pcevent.Cathodech >= 0 && pcevent.Cathodech < 24) - pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, dE_gas}); + // Energy1/Energy2 are per-wire ADC. When multi1/multi2 > 1 the charge is + // shared across several wires, so a single wire's Energy here is only + // its (event-dependent) share of the total -- not comparable to the + // full predicted dE_gas. Restrict each side to clean, single-wire + // clusters so ADC and dE_gas both refer to the whole deposited charge. + // z now comes from SX3 independently of multi1/multi2, so anode and + // cathode purity are gated independently (mirrors the proton + // accumulator's tryEvent). + if (pcevent.multi1 == 1 && pcevent.Anodech >= 0 && pcevent.Anodech < 24) + pcCalibWritePoint(pcevent.Anodech, pcevent.Energy1, dE_gas); + if (pcevent.multi2 >= 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24) + pcCalibWritePoint(24 + pcevent.Cathodech, pcevent.Energy2, dE_gas); } } @@ -848,7 +925,7 @@ inline double predictElasticEnergy(Kinematics &kin, double angle3_deg, double t3 // reading the PC's own (uncalibrated) energy, predicts the alpha's energy // from elastic kinematics at the reconstructed angle/vertex and walks that // forward through the gas to a predicted dE_gas -- same target quantity as -// the source calibration, merged into the same pcCalibData[wire] clouds. +// the source calibration, streamed into the same per-wire output file. inline void pcEnergyCalibrationAccumulateProton(const std::vector &PC_Events, const std::vector &QQQ_Events, const std::vector &SX3_Events) { if (!ta_foil_run) @@ -902,9 +979,9 @@ inline void pcEnergyCalibrationAccumulateProton(const std::vector &PC_Eve // Single-wire only, so each wire's cloud maps its own charge (see the // source-run accumulator for the same reasoning). if (pcevent.multi1 == 1 && pcevent.Anodech >= 0 && pcevent.Anodech < 24) - pcCalibData[pcevent.Anodech].push_back({pcevent.Energy1, dE_gas}); - if (pcevent.Cathodech >= 0 && pcevent.Cathodech < 24) - pcCalibData[24 + pcevent.Cathodech].push_back({pcevent.Energy2, dE_gas}); + pcCalibWritePoint(pcevent.Anodech, pcevent.Energy1, dE_gas); + if (pcevent.multi2 == 1 && pcevent.Cathodech >= 0 && pcevent.Cathodech < 24) + pcCalibWritePoint(24 + pcevent.Cathodech, pcevent.Energy2, dE_gas); }; for (const auto &pcevent : PC_Events) @@ -1396,6 +1473,38 @@ Bool_t TrackRecon::Process(Long64_t entry) for (const auto &aCluster : aClusters) { + // A2 anode charge-sharing ratio: for exactly-two-fired-wire (multi-anode) + // clusters, the ratio of the smaller to the larger wire's energy. This + // characterises how the deposited charge splits between the two straddled + // anode wires (~0 = one wire dominates, ~1 = even split) and is filled once + // per anode cluster (before the cathode loop) so it isn't inflated by the + // cathode multiplicity. + if (aCluster.size() == 2) + { + double ae0 = std::get<1>(aCluster[0]); + double ae1 = std::get<1>(aCluster[1]); + double alo = std::min(ae0, ae1); + double ahi = std::max(ae0, ae1); + if (ahi > 0.0) + { + double aratio = alo / ahi; + plotter->Fill1D("A2_anode_ratio", 120, 0, 1.2, aratio, "hGMPC"); + plotter->Fill2D("A2_anode_ratio_vs_sum", 800, 0, 40000, 120, 0, 1.2, ae0 + ae1, aratio, "hGMPC"); + plotter->Fill2D("A2_anode_ratio_vs_lowerIndex", 24, 0, 24, 120, 0, 1.2, + std::min(std::get<0>(aCluster[0]), std::get<0>(aCluster[1])), aratio, "hGMPC"); + } + // Raw-vs-raw sanity check for the "A2 total = 2x A1 total" question, BEFORE + // any energy-calibration math touches it: if the raw (gain-matched-only) + // 2-wire sum already peaks at ~2x the raw single-wire distribution, the + // effect is upstream of the calibration code (neighbor-wire crosstalk on + // the raw signal, or a genuine double-hit) rather than introduced by how + // the calibration sums wires. + plotter->Fill1D("Raw_A2_AnodeSum", 800, 0, 40000, ae0 + ae1, "hGMPC"); + } + else if (aCluster.size() == 1) + { + plotter->Fill1D("Raw_A1_AnodeSum", 800, 0, 40000, std::get<1>(aCluster[0]), "hGMPC"); + } for (const auto &cCluster : cClusters) { if (aCluster.size() == 0) @@ -1419,69 +1528,51 @@ Bool_t TrackRecon::Process(Long64_t entry) if (pcEnergyCalibLoaded) { Event PCEventCalibrated = PCEvent; - + // Calibrate EACH anode wire's charge with its own slope and sum the + // results, rather than applying the first wire's factor to the raw + // cluster sum: the alpha's total gas dE is the charge summed over the + // wires it straddles, and whether it lands on one wire or two is + // phi-correlated, so a single-factor-on-sum made the calibrated anode + // energy depend on phi (unphysical: dE depends only on theta,z). + // + // The intercept, though, is a per-EVENT baseline (pedestal/threshold + // offset from the single-wire fit, where one wire's ADC stood for the + // whole dE_gas) -- not a per-wire quantity. Adding intercept[wi] once + // per fired wire double-counts that baseline for multi-wire clusters, + // reintroducing a multiplicity- (hence indirectly phi-) dependent bias + // through the back door. Apply it exactly once, from the primary + // (max-energy, i.e. aCluster[0]) wire. double anodeCalibSum = 0.0; - double primaryIntercept = 0.0; - - // Vector to track individual calibrated fractional energies for A2 charge sharing - std::vector calibWireEnergies; + std::vector calibWireEnergies; // per-wire slope*ADC (no intercept), for the A2 ratio below calibWireEnergies.reserve(aCluster.size()); - - if (!aCluster.empty()) - { - // Use the intercept of the primary (maximum energy) wire exactly ONCE. - // The first element of aCluster from Make_Clusters is the maximum energy wire. - int primaryWire = std::get<0>(aCluster[0]); - if (primaryWire >= 0 && primaryWire < 24) - { - primaryIntercept = pcEnergyIntercept[primaryWire]; - } - } - - // Apply per-wire gain slopes to the fractional charge, but DO NOT add the intercept yet for (const auto &w : aCluster) { int wi = std::get<0>(w); - if (wi >= 0 && wi < 24) - { - // Note: std::get<1>(w) already has the relative matching 'pcSlope' applied - double wCalibE = pcEnergySlope[wi] * std::get<1>(w); - anodeCalibSum += wCalibE; - calibWireEnergies.push_back(wCalibE); - } + double wCalibE = (wi >= 0 && wi < 24) ? pcEnergySlope[wi] * std::get<1>(w) : 0.0; + anodeCalibSum += wCalibE; + calibWireEnergies.push_back(wCalibE); } - - // Apply the baseline offset (intercept) exactly once for the whole event + int primaryAnodeWire = std::get<0>(aCluster[0]); + double primaryIntercept = (primaryAnodeWire >= 0 && primaryAnodeWire < 24) ? pcEnergyIntercept[primaryAnodeWire] : 0.0; PCEventCalibrated.Energy1 = anodeCalibSum + primaryIntercept; - - // Cathode uses the single max wire (cpMaxE) -- indexed by z, so it's already phi-consistent + // Cathode uses the single max wire (cpMaxE) -- indexed by z, so it's + // already phi-consistent; leave it as-is. PCEventCalibrated.Energy2 = pcEnergySlope[24 + PCEvent.Cathodech] * cpMaxE + pcEnergyIntercept[24 + PCEvent.Cathodech]; - PC_Events_calibrated.push_back(PCEventCalibrated); - // --------------------------------------------------------- - // A2 Anode Charge-Sharing Diagnostics - // --------------------------------------------------------- - if (aCluster.size() == 2 && calibWireEnergies.size() == 2) + // Calibrated-energy A2 charge-sharing ratio: same diagnostic as the raw + // version above (A2_anode_ratio), but on the per-wire CALIBRATED shares, + // so miscalibration between the two wires shows up as a ratio pulled + // away from what the raw-ADC ratio would give. Checks phi/energy + // dependence directly on the quantity that actually feeds Energy1. + if (calibWireEnergies.size() == 2) { - double e1 = calibWireEnergies[0]; - double e2 = calibWireEnergies[1]; - - double eSmaller = std::min(e1, e2); - double eLarger = std::max(e1, e2); - - // Ratio goes from 0.0 (almost all charge on one wire) to 1.0 (perfectly shared) + double eSmaller = std::min(calibWireEnergies[0], calibWireEnergies[1]); + double eLarger = std::max(calibWireEnergies[0], calibWireEnergies[1]); double ratio = (eLarger > 0.0) ? (eSmaller / eLarger) : 0.0; - - // Plotting the 1D ratio distribution plotter->Fill1D("Calib_A2_AnodeRatio", 200, 0.0, 1.0, ratio, "hCalibPC"); - - // Checking phi-dependence: Should peak at ratio=1.0 strictly between wires - // and approach 0.0 when passing directly over a wire plotter->Fill2D("Calib_A2_AnodeRatio_vs_Phi", 360, -180, 180, 200, 0.0, 1.0, PCEvent.pos.Phi() * 180.0 / M_PI, ratio, "hCalibPC"); - - // Energy dependence: Verifies if charge induction geometry changes via delta-rays/track-length plotter->Fill2D("Calib_A2_AnodeRatio_vs_TotalE", 400, 0, 10, 200, 0.0, 1.0, PCEventCalibrated.Energy1, ratio, "hCalibPC"); } @@ -1541,14 +1632,18 @@ Bool_t TrackRecon::Process(Long64_t entry) if (pcEnergyCalibLoaded) { - // Per-wire-then-sum, phi-independent (same reasoning as the crossover branch). + // Per-wire-then-sum, intercept applied once from the primary wire -- + // same reasoning as the crossover branch above. 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) + pcEnergyIntercept[wi]; + anodeCalibSum += pcEnergySlope[wi] * std::get<1>(w); } + int primaryAnodeWireA1C0 = std::get<0>(aCl[0]); + double primaryInterceptA1C0 = (primaryAnodeWireA1C0 >= 0 && primaryAnodeWireA1C0 < 24) ? pcEnergyIntercept[primaryAnodeWireA1C0] : 0.0; + anodeCalibSum += primaryInterceptA1C0; 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 @@ -1560,20 +1655,28 @@ Bool_t TrackRecon::Process(Long64_t entry) // 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 && !ta_foil_run) + if (doPCEnergyCalibration && source_run) { double path = pathLengthCm(source_pos_a1c0, pc); double e_rem = evalElossForward(MeV_to_cm_spl, cm_to_MeV_spl, pc_calib_alpha_source_mev, path); double dE_gas = pc_calib_alpha_source_mev - e_rem; if (std::isfinite(dE_gas) && dE_gas > 0.0) - pcCalibData[anodeIdx].push_back({apSumE, dE_gas}); + pcCalibWritePoint(anodeIdx, apSumE, dE_gas); } } } if (doPCEnergyCalibration) { - pcEnergyCalibrationAccumulate(PC_Events); + // pcEnergyCalibrationAccumulate assumes source_pos = the FIXED alpha-source + // position -- only true for the source runs. On a proton-scattering run, + // source_vertex is the beam-entrance placeholder, not a real source, so + // running this there silently wrote bogus points (huge path length, usually + // clamped to the full 5.486 MeV) into the same pool as the real + // kinematics-derived proton points. Gate it to non-proton-campaign runs; + // pcEnergyCalibrationAccumulateProton already self-gates the other way. + if (source_run) + pcEnergyCalibrationAccumulate(PC_Events, SX3_Events, QQQ_Events); pcEnergyCalibrationAccumulateProton(PC_Events, QQQ_Events, SX3_Events); } @@ -1745,26 +1848,14 @@ void TrackRecon::Terminate() { plotter->FlushToDisk(10); - if (doPCEnergyCalibration) + if (doPCEnergyCalibration && pcCalibOutFile.is_open()) { - gSystem->mkdir("pc_calib_raw", kTRUE); // kTRUE = create parents, no-op if it exists - std::string tag = getenv("RUN_NUMBER") ? std::string("run") + getenv("RUN_NUMBER") - : dataset + "_pid" + std::to_string(getpid()); - std::string outname = "pc_calib_raw/points_" + tag + ".dat"; - std::ofstream outfile(outname); - outfile << std::scientific << std::setprecision(6); - long long nPoints = 0; - for (int wire = 0; wire < 48; ++wire) - { - for (const auto &p : pcCalibData[wire]) - { - outfile << wire << " " << p.first << " " << p.second << "\n"; - ++nPoints; - } - } - outfile.close(); - std::cout << "PC energy calibration: appended " << nPoints << " raw points to " << outname - << " -- run pccal/fit_pc_energy_calibration.C once all calibration runs are done" + // Points were streamed directly to disk as they were generated (see + // pcCalibWritePoint / Begin()) instead of buffered in memory for the + // whole run, so there's nothing left to write here -- just close it. + pcCalibOutFile.close(); + std::cout << "PC energy calibration: closed raw points file -- run " + << "pccal/fit_pc_energy_calibration.C once all calibration runs are done" << " to (re)produce pc_energy_calibration.dat" << std::endl; } } @@ -1890,7 +1981,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector &QQQ_ const bool hasCathode = (pcevent.Cathodech >= 0); const double totalE = hasCathode ? (pcevent.Energy1 + pcevent.Energy2) : pcevent.Energy1; if (hasCathode) - plotter->Fill2D("Calib_AnodeE_vs_CathodeE_a1c1andup", 800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC"); + plotter->Fill2D("Calib_AnodeE_vs_CathodeE_a1c1andup",800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC"); for (const std::string &t : {std::string(""), topo}) { plotter->Fill2D("Calib_AnodeE_vs_AnodeIndex" + t, 24, 0, 24, 800, 0, 3, pcevent.Anodech, pcevent.Energy1, "hCalibPC"); @@ -1901,7 +1992,7 @@ void pcCalibratedHistograms(HistPlotter *plotter, const std::vector &QQQ_ { plotter->Fill2D("Calib_CathodeE_vs_CathodeIndex" + t, 24, 0, 24, 800, 0, 3, pcevent.Cathodech, pcevent.Energy2, "hCalibPC"); plotter->Fill1D("Calib_CathodeE" + t, 800, 0, 3, pcevent.Energy2, "hCalibPC"); - plotter->Fill2D("Calib_AnodeE_vs_CathodeE" + t, 800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC"); + plotter->Fill2D("Calib_AnodeE_vs_CathodeE" + t,800, 0, 3, 800, 0, 3, pcevent.Energy1, pcevent.Energy2, "hCalibPC"); } for (const auto &qqqevent : QQQ_Events) diff --git a/pc_energy_calibration.dat b/pc_energy_calibration.dat new file mode 100644 index 0000000..9e82e0a --- /dev/null +++ b/pc_energy_calibration.dat @@ -0,0 +1,48 @@ +0 8.024604e-05 0.000000e+00 +1 7.374085e-05 0.000000e+00 +2 9.583027e-05 0.000000e+00 +3 9.008263e-05 0.000000e+00 +4 8.452456e-05 0.000000e+00 +5 1.835247e-04 0.000000e+00 +6 5.125204e-05 0.000000e+00 +7 1.081530e-04 0.000000e+00 +8 1.269708e-04 0.000000e+00 +9 2.018195e-04 0.000000e+00 +10 1.089829e-04 0.000000e+00 +11 8.956875e-05 0.000000e+00 +12 1.000000e+00 0.000000e+00 +13 7.299374e-05 0.000000e+00 +14 7.270788e-05 0.000000e+00 +15 8.406974e-05 0.000000e+00 +16 8.997978e-05 0.000000e+00 +17 4.542671e-05 0.000000e+00 +18 4.755692e-05 0.000000e+00 +19 2.138104e-04 0.000000e+00 +20 9.264217e-05 0.000000e+00 +21 1.701897e-04 0.000000e+00 +22 2.162763e-04 0.000000e+00 +23 3.219591e-04 0.000000e+00 +24 7.364520e-05 0.000000e+00 +25 9.343009e-05 0.000000e+00 +26 1.013807e-04 0.000000e+00 +27 1.495574e-04 0.000000e+00 +28 1.166817e-04 0.000000e+00 +29 9.196472e-05 0.000000e+00 +30 1.030198e-04 0.000000e+00 +31 1.007519e-04 0.000000e+00 +32 1.026175e-04 0.000000e+00 +33 8.109499e-05 0.000000e+00 +34 8.137172e-05 0.000000e+00 +35 8.502646e-05 0.000000e+00 +36 1.042813e-04 0.000000e+00 +37 3.213148e-04 0.000000e+00 +38 9.069578e-05 0.000000e+00 +39 1.600775e-04 0.000000e+00 +40 1.357355e-04 0.000000e+00 +41 1.812490e-04 0.000000e+00 +42 1.186314e-04 0.000000e+00 +43 1.196526e-04 0.000000e+00 +44 2.486180e-04 0.000000e+00 +45 8.347962e-05 0.000000e+00 +46 7.408129e-05 0.000000e+00 +47 8.388388e-05 0.000000e+00 diff --git a/pccal/fit_pc_energy_calibration.C b/pccal/fit_pc_energy_calibration.C index 1de6b87..1f9a050 100644 --- a/pccal/fit_pc_energy_calibration.C +++ b/pccal/fit_pc_energy_calibration.C @@ -1,29 +1,26 @@ -// Aggregates every run's raw PC energy-calibration points (written by -// TrackRecon.C's Terminate() into pc_calib_raw/points_*.dat, one file per run -// so parallel jobs never collide) into a single per-wire linear fit, pooling -// BOTH the 17F and 27Al datasets together -- the PC's wire-level gain isn't -// expected to differ meaningfully between the two campaigns, so calibration -// statistics from either dataset's alpha-source or proton-scattering runs are -// combined into one fit rather than kept separate. -// -// Run non-interactively from the top-level ANASEN-Analysis directory once all -// desired calibration runs have completed: -// root -q -l -b -e '.L pccal/fit_pc_energy_calibration.C' -e 'fit_pc_energy_calibration()' -// -// Writes pc_energy_calibration.dat (wire slope intercept, 48 rows), which -// TrackRecon.C's Begin() loads automatically on the next run to populate -// PC_Events_calibrated. Re-run this any time pc_calib_raw/ has new points -- -// it always re-reads and re-fits everything currently on disk, so it's safe -// to call repeatedly as more calibration runs accumulate. #include #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include +#include +#include -void fit_pc_energy_calibration() +// 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 = "") { std::vector> pts[48]; // [wire] -> (ADC, dE_gas MeV) @@ -37,55 +34,180 @@ void fit_pc_energy_calibration() } int nFiles = 0; + long long nOverflowCut = 0; TIter next(files); TSystemFile *f; + + // Read all points from all files into the single pooled array while ((f = (TSystemFile *)next())) { TString name = f->GetName(); if (f->IsDirectory() || !name.BeginsWith("points_") || !name.EndsWith(".dat")) continue; + + // Apply dataset safeguard if requested + if (!dataset_filter.empty() && !name.Contains(dataset_filter.c_str())) + continue; + std::ifstream infile(std::string("pc_calib_raw/") + name.Data()); if (!infile.is_open()) continue; + int wire; double adc, dE_gas; while (infile >> wire >> adc >> dE_gas) + { if (wire >= 0 && wire < 48) - pts[wire].push_back({adc, dE_gas}); + { + if (adc < 64000.0) { + pts[wire].push_back({adc, dE_gas}); + } else { + nOverflowCut++; + } + } + } ++nFiles; } - std::cout << "fit_pc_energy_calibration: read " << nFiles << " run file(s) from pc_calib_raw/" << std::endl; + + std::cout << "fit_pc_energy_calibration: read " << nFiles + << " 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; + + // --- Setup ROOT Diagnostic Graphics --- + gROOT->SetBatch(kTRUE); // Run silently without popping up windows + gStyle->SetOptStat(0); + gStyle->SetPalette(kBird); + + // Create output directory for individual PNGs + gSystem->mkdir("pc_calib_plots", kTRUE); + + TFile *fOut = new TFile("pc_calib_diagnostics.root", "RECREATE"); + + TCanvas *cAnodes = new TCanvas("cAnodes", "Anode Calibrations", 1800, 1200); + cAnodes->Divide(6, 4, 0.01, 0.01); + + TCanvas *cCathodes = new TCanvas("cCathodes", "Cathode Calibrations", 1800, 1200); + cCathodes->Divide(6, 4, 0.01, 0.01); std::ofstream outfile("pc_energy_calibration.dat"); outfile << std::scientific << std::setprecision(6); + + // --- Fit and Plot each wire --- for (int wire = 0; wire < 48; ++wire) { double slope = 1.0, intercept = 0.0; - bool ok = pts[wire].size() >= 2; + 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 + + for (const auto &p : pts[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; + } + if (ok) { - double sx = 0, sy = 0, sxx = 0, sxy = 0; - for (const auto &p : pts[wire]) + if (std::isfinite(sxx) && std::abs(sxx) > 1e-12) { - sx += p.first; - sy += p.second; - sxx += p.first * p.first; - sxy += p.first * p.second; + slope = sxy / sxx; + intercept = 0.0; // Forced mathematically } - double n = static_cast(pts[wire].size()); - double denom = n * sxx - sx * sx; - ok = std::isfinite(denom) && std::abs(denom) > 1e-12; - if (ok) + else { - slope = (n * sxy - sx * sy) / denom; - intercept = (sy - slope * sx) / n; + ok = false; } } - if (!ok) - std::cerr << "fit_pc_energy_calibration: wire " << wire << " has too few points (" << pts[wire].size() + + 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; + } + outfile << wire << " " << slope << " " << intercept << "\n"; + + // 2. Generate and Fill 2D Density Histogram + if (maxX <= 0) maxX = 64000.0; + if (maxY <= 0) maxY = 10.0; + + TString wName = Form("%s %02d", wire < 24 ? "Anode" : "Cathode", wire < 24 ? wire : wire - 24); + TH2D *h2 = new TH2D(Form("h2_wire_%d", wire), + wName + "; ADC; dE_{gas} (MeV)", + 150, 0, maxX * 1.05, + 150, 0, maxY * 1.05); + + for (const auto &p : pts[wire]) { + h2->Fill(p.first, p.second); + } + + // 3. Setup Fit Line and Stats Box + TF1 *fitLine = nullptr; + TPaveText *pt = nullptr; + + if (n > 0) { + 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->SetLineWidth(2); + + pt = new TPaveText(0.15, 0.75, 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)"); + } + + // ------------------------------------------------------------ + // 4. Save High-Res Individual PNG + // ------------------------------------------------------------ + TCanvas cTemp("cTemp", "cTemp", 800, 600); + cTemp.SetGridx(); + cTemp.SetGridy(); + if (n > 0) { + h2->Draw("COLZ"); + fitLine->Draw("SAME"); + pt->Draw(); + } else { + h2->SetTitle(wName + " (DEAD/NO DATA)"); + h2->Draw(); + } + cTemp.SaveAs(Form("pc_calib_plots/wire_%02d.png", wire)); + + // ------------------------------------------------------------ + // 5. Draw onto Global PDF Canvas + // ------------------------------------------------------------ + TVirtualPad* pad = (wire < 24) ? cAnodes->cd(wire + 1) : cCathodes->cd(wire - 24 + 1); + pad->SetGridx(); + pad->SetGridy(); + + if (n > 0) + { + h2->DrawClone("COLZ"); + fitLine->DrawClone("SAME"); + pt->DrawClone(); + } + else + { + h2->DrawClone(); + } + + // Write histogram to ROOT file + fOut->cd(); + h2->Write(); } + outfile.close(); + + fOut->Close(); + std::cout << "fit_pc_energy_calibration: wrote pc_energy_calibration.dat" << 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: global contact sheets saved to pc_calib_anodes.pdf and pc_calib_cathodes.pdf" << std::endl; +} \ No newline at end of file diff --git a/run_27Al.sh b/run_27Al.sh index 2edb596..d3819d3 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -6,12 +6,13 @@ export OUT_DIR="Output_27Al" export reactiondata=1 export CO2percent=3 export pressure_in_torr=250 -export CATHODE_GAIN=1 +export CATHODE_GAIN=3.0 +export source_vertex=-200.0 export BEAM_AXIS_X=-15 export BEAM_AXIS_Y=-5 # Clean up previous runs -rm -f ${OUT_DIR}/results_run*.root output_27Al.root +rm -f ${OUT_DIR}/*.root echo "Pre-compiling TrackRecon.C safely on a single core..." root -q -l -b -e '.L TrackRecon.C++O' @@ -23,7 +24,6 @@ process_run() { local infile="../ANASEN_analysis/data/${DATASET}_Data/${prefix}${wrun}_mapped.root" local out="${outdir}/results_run${wrun}.root" - # Ensure the directory exists so ROOT doesn't fail silently mkdir -p "$outdir" root -q -l -b -x "$infile" \ @@ -39,19 +39,28 @@ process_run() { export -f process_run echo "Starting parallel processing..." -parallel --bar -j 4 process_run ::: {50..52} +time parallel --bar -j 8 process_run ::: {50..59} +time parallel --bar -j 8 process_run ::: 62 63 66 67 73 74 +# time parallel --bar -j 1 run_once {1} ::: 68 +# time parallel --bar -j 6 run_once {1} ::: {78..89} echo "Merging files..." hadd -k -j 4 ${OUT_DIR}/output_27Al.root ${OUT_DIR}/results_run*.root +rootbrowse ${OUT_DIR}/output_27Al.root + unset DATASET +unset PREFIX +unset OUT_DIR unset reactiondata unset CO2percent unset pressure_in_torr unset CATHODE_GAIN +unset source_vertex unset A1C1_LOWBAND_RFACTOR unset A1C1_Z_SCALE_QQQ unset A1C1_Z_OFF_QQQ unset A1C1_Z_OFF_SX3 unset BEAM_AXIS_X -unset BEAM_AXIS_Y \ No newline at end of file +unset BEAM_AXIS_Y +echo "Script execution finished." \ No newline at end of file diff --git a/run_tr.sh b/run_tr.sh index 422c8ce..906a130 100644 --- a/run_tr.sh +++ b/run_tr.sh @@ -63,7 +63,7 @@ if [[ 1 -eq 1 ]]; then export pressure_in_torr=350 rm -f ${OUT_DIR}/all.root echo "Processing 27Al alpha+gas runs..." - export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope" + # export source_vertex=-5.36; export timecut_low=12.0; export timecut_high=119.0; process_run 9 "$slope" unset timecut_high export source_vertex=53.44; export timecut_low=400.0; process_run 12 "$slope" unset Gain @@ -89,7 +89,7 @@ if [[ 1 -eq 1 ]]; then fi # --- Block 5: 27Al Protons+Gas Runs (15, 17-22) --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then # export CO2percent=4 export DATASET="27Al" @@ -101,15 +101,15 @@ if [[ 1 -eq 0 ]]; then echo "Starting parallel processing for 27Al proton runs..." # process_run 18 - parallel --bar -j 8 process_run ::: 15 {17..22} - # parallel --bar -j 8 process_run ::: {17..22} + # parallel --bar -j 8 process_run ::: 15 {17..22} + parallel --bar -j 8 process_run ::: {17..22} hadd -j 4 -k ${OUT_DIR}/Al_protons.root ${OUT_DIR}/results_run0{15..22}.root unset CATHODE_GAIN # exit fi # --- Block 6: 17F Proton Data --- -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export DATASET="17F" export PREFIX="ProtonRun_" export OUT_DIR="Output_p"