diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Action.yaml b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Action.yaml new file mode 100644 index 00000000..f4c2bc69 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Action.yaml @@ -0,0 +1,34 @@ +algorithm: + name: 'clas12::ProtonEnergyLossCorrection' + +actions: + + - name: Transform + type: transformer + rank: scalar + inputs: + - name: pid + type: int + - name: status + type: int + - name: run + type: int + - name: px + type: double + cast: float + - name: py + type: double + cast: float + - name: pz + type: double + cast: float + outputs: + - name: px + type: double + cast: float + - name: py + type: double + cast: float + - name: pz + type: double + cast: float \ No newline at end of file diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc new file mode 100644 index 00000000..82d21734 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -0,0 +1,353 @@ +#include "Algorithm.h" + +#include +#include + +#include + +namespace iguana::clas12::rga { + + REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); + + namespace { + // Keep an angle in degrees within [0, 360). + double WrapDeg360(double x) { + while(x >= 360.0) { x -= 360.0; } + while(x < 0.0) { x += 360.0; } + return x; + } + } + + // ----------------------------------------------------------------------------- + // Detector-region helpers (based on REC::Particle status) + // ----------------------------------------------------------------------------- + // + // CLAS12 status definitions: + // - FD tracks: abs(status) in [2000, 4000) + // - CD tracks: abs(status) in [4000, 5000) + // + bool ProtonEnergyLossCorrection::IsFD(int status) { + int s = std::abs(status); + return (s >= 2000 && s < 4000); + } + + bool ProtonEnergyLossCorrection::IsCD(int status) { + int s = std::abs(status); + return (s >= 4000 && s < 5000); + } + + // ----------------------------------------------------------------------------- + // math helpers (p, theta, phi) + // ----------------------------------------------------------------------------- + double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) { + return std::sqrt(px * px + py * py + pz * pz); + } + + // theta from momentum vector, in degrees. + // cos(theta) = pz / |p| + double ProtonEnergyLossCorrection::ThetaDeg(double px, double py, double pz) { + double r = Pmag(px, py, pz); + if(r <= 0.0) { + return 0.0; + } + + double c = pz / r; + if(c > 1.0) { + c = 1.0; + } + if(c < -1.0) { + c = -1.0; + } + + return (180.0 / M_PI) * std::acos(c); + } + + + // Returns phi in [0,360). + double ProtonEnergyLossCorrection::PhiDeg(double px, double py) { + double phi = (180.0 / M_PI) * std::atan2(px, py); + phi = phi - 90.0; + if(phi < 0.0) { + phi = 360.0 + phi; + } + phi = 360.0 - phi; + + return WrapDeg360(phi); + } + + // Convert from (p, theta_deg, phi_deg) back to Cartesian (px,py,pz). + void ProtonEnergyLossCorrection::SphericalToCartesian(double p, double theta_deg, + double phi_deg, double& px, double& py, double& pz) { + double theta = theta_deg * (M_PI / 180.0); + double phi = phi_deg * (M_PI / 180.0); + + px = p * std::sin(theta) * std::cos(phi); + py = p * std::sin(theta) * std::sin(phi); + pz = p * std::cos(theta); + } + + // Evaluate polynomial c0 + c1*x + c2*x^2 + ... + double ProtonEnergyLossCorrection::EvalPoly(Poly const& p, double x) { + double sum = 0.0; + double xn = 1.0; + for(double c : p.c) { + sum += c * xn; + xn *= x; + } + return sum; + } + + // ----------------------------------------------------------------------------- + // Period lookup + // ----------------------------------------------------------------------------- + // + // Given a run number, find the first period def whose run_ranges contain it. + // If no match is found, return nullptr and the algorithm leaves the row unchanged. + ProtonEnergyLossCorrection::PeriodDef + const* ProtonEnergyLossCorrection::FindPeriod(int run) const { + for(auto const& kv : m_periods) { + auto const& def = kv.second; + for(auto const& rr : def.run_ranges) { + if(run >= rr.first && run <= rr.second) { + return &def; + } + } + } + return nullptr; + } + + // ----------------------------------------------------------------------------- + // ConfigHook: read Config.yaml and populate m_periods + // ----------------------------------------------------------------------------- + // + // YAML schema expected: + // + // clas12: + // ProtonEnergyLossCorrection: + // periods: + // : + // run_ranges: + // - [min, max] + // - [min, max] + // FD: + // A_p: [ ... ] + // B_p: [ ... ] + // ... + // CD: + // A_p: [ ... ] + // ... + // + // For each period, load FD and CD RegionCoeffs. Each RegionCoeffs holds + // polynomial parameterizations in theta for A/B/C of p/theta/phi corrections. + // many A/B/C may be zero!! + void ProtonEnergyLossCorrection::ConfigHook() { + m_periods.clear(); + + // algorithm's installed Config.yaml + std::string const cfg_path = + GetConfig()->FindFile("algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml"); + + YAML::Node root = YAML::LoadFile(cfg_path); + YAML::Node node = root["clas12"]["ProtonEnergyLossCorrection"]; + YAML::Node periods_node = node["periods"]; + + // Helper: load a polynomial array by name (e.g. "A_p") from a region node. + auto load_poly = [](YAML::Node const& region_node, char const* name) -> Poly { + YAML::Node a = region_node[name]; + Poly p; + p.c.reserve((size_t)a.size()); + for(auto v : a) { + p.c.push_back(v.as()); + } + return p; + }; + + // Helper: load all A/B/C for p/theta/phi for a given detector region node. + auto load_region = [&](YAML::Node const& rnode) -> RegionCoeffs { + RegionCoeffs rc; + + rc.A_p = load_poly(rnode, "A_p"); + rc.B_p = load_poly(rnode, "B_p"); + rc.C_p = load_poly(rnode, "C_p"); + + rc.A_theta = load_poly(rnode, "A_theta"); + rc.B_theta = load_poly(rnode, "B_theta"); + rc.C_theta = load_poly(rnode, "C_theta"); + + rc.A_phi = load_poly(rnode, "A_phi"); + rc.B_phi = load_poly(rnode, "B_phi"); + rc.C_phi = load_poly(rnode, "C_phi"); + + return rc; + }; + + // Iterate periods dictionary: key -> data + for(auto it = periods_node.begin(); it != periods_node.end(); ++it) { + std::string key = it->first.as(); + YAML::Node pnode = it->second; + + PeriodDef def; + + for(auto rr : pnode["run_ranges"]) { + def.run_ranges.push_back({rr[0].as(), rr[1].as()}); + } + + def.fd = load_region(pnode["FD"]); + def.cd = load_region(pnode["CD"]); + + m_periods[key] = def; + } + } + + // ----------------------------------------------------------------------------- + // StartHook: cache bank indices + // ----------------------------------------------------------------------------- + void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) { + m_b_rec_particle = GetBankIndex(banks, "REC::Particle"); + m_b_run_config = GetBankIndex(banks, "RUN::config"); + } + + // ----------------------------------------------------------------------------- + // RunHook: apply correction to each matching REC::Particle row + // ----------------------------------------------------------------------------- + bool ProtonEnergyLossCorrection::RunHook(hipo::banklist& banks) const { + auto& rec = GetBank(banks, m_b_rec_particle, "REC::Particle"); + auto& run = GetBank(banks, m_b_run_config, "RUN::config"); + + int runnum = run.getInt("run", 0); + + // Loop all rows and correct in place. + for(auto const& row : rec.getRowList()) { + + // Select protons only. + int pid = rec.getInt("pid", row); + if(pid != 2212) { + continue; + } + + // Select FD or CD only (based on status). + int status = rec.getInt("status", row); + if(!IsFD(status) && !IsCD(status)) { + continue; + } + + // Read momentum components. + double px = rec.getFloat("px", row); + double py = rec.getFloat("py", row); + double pz = rec.getFloat("pz", row); + + // Transform returns corrected values (or unchanged if not applicable). + auto [px_new, py_new, pz_new] = Transform(pid, status, runnum, px, py, pz); + + // Write back to bank. + rec.putFloat("px", row, static_cast(px_new)); + rec.putFloat("py", row, static_cast(py_new)); + rec.putFloat("pz", row, static_cast(pz_new)); + } + + return true; + } + + void ProtonEnergyLossCorrection::StopHook() { + } + + // ----------------------------------------------------------------------------- + // Transform: core correction logic (FD vs CD functional forms) + // ----------------------------------------------------------------------------- + std::tuple + ProtonEnergyLossCorrection::Transform(int const pid, int const status, int const run, + vector_element_t const px_in, vector_element_t const py_in, + vector_element_t const pz_in) const { + // only protons are corrected. + if(pid != 2212) { + return {px_in, py_in, pz_in}; + } + + // Detector region classification from status. + bool is_fd = IsFD(status); + bool is_cd = IsCD(status); + if(!is_fd && !is_cd) { + return {px_in, py_in, pz_in}; + } + + // Period lookup from run number. If run is not recognized, do nothing. + PeriodDef const* period = FindPeriod(run); + if(!period) { + return {px_in, py_in, pz_in}; + } + + // Compute kinematics from input momentum. + double px = (double)px_in; + double py = (double)py_in; + double pz = (double)pz_in; + + double p = Pmag(px, py, pz); + if(p <= 0.0) { + return {px_in, py_in, pz_in}; + } + double theta = ThetaDeg(px, py, pz); // degrees + double phi = PhiDeg(px, py); // degrees in [0,360) + + // Choose FD vs CD coefficients. + RegionCoeffs const& coeffs = is_fd ? period->fd : period->cd; + + // Evaluate A/B/C polynomials at this theta. + double A_p = EvalPoly(coeffs.A_p, theta); + double B_p = EvalPoly(coeffs.B_p, theta); + double C_p = EvalPoly(coeffs.C_p, theta); + + double A_theta = EvalPoly(coeffs.A_theta, theta); + double B_theta = EvalPoly(coeffs.B_theta, theta); + double C_theta = EvalPoly(coeffs.C_theta, theta); + + double A_phi = EvalPoly(coeffs.A_phi, theta); + double B_phi = EvalPoly(coeffs.B_phi, theta); + double C_phi = EvalPoly(coeffs.C_phi, theta); + + // Apply corrections. + // + // Forward Detector (FD): + // p_new = p + A_p + B_p/p + C_p/p^2 + // theta_new = theta + A_theta + B_theta/theta + C_theta/theta^2 + // phi_new = phi + A_phi + B_phi/phi + C_phi/phi^2 + // + // Central Detector (CD): + // p_new = p + A_p + B_p*p + C_p*p^2 + // not all A_, B_ and C_ need be nonzero + double p_new = p; + double theta_new = theta; + double phi_new = phi; + + if(is_fd) { + p_new += A_p + (B_p / p) + (C_p / (p * p)); + + if(theta != 0.0) { + theta_new += A_theta + (B_theta / theta) + (C_theta / (theta * theta)); + } + if(phi != 0.0) { + phi_new += A_phi + (B_phi / phi) + (C_phi / (phi * phi)); + } + } else { + // is_cd + p_new += A_p + (B_p * p) + (C_p * p * p); + + if(theta != 0.0) { + theta_new += A_theta + (B_theta / theta) + (C_theta / (theta * theta)); + } + if(phi != 0.0) { + phi_new += A_phi + (B_phi / phi) + (C_phi / (phi * phi)); + } + } + + phi_new = WrapDeg360(phi_new); + + // Convert back to Cartesian using the corrected spherical variables. + double px_out = 0.0; + double py_out = 0.0; + double pz_out = 0.0; + SphericalToCartesian(p_new, theta_new, phi_new, px_out, py_out, pz_out); + + return {px_out, py_out, pz_out}; + } + +} \ No newline at end of file diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h new file mode 100644 index 00000000..a9a333f1 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -0,0 +1,117 @@ +#pragma once + +// ProtonEnergyLossCorrection +// +// Apply momentum and angular corrections to reconstructed protons in CLAS12 RGA +// data; runs event-by-event and modifies REC::Particle in place. +// For each row in REC::Particle: +// - If pid == 2212 (proton) AND status indicates FD or CD (should always), +// compute p, theta, phi from (px,py,pz), apply period-dependent corrections, +// then write back corrected (px,py,pz). +// +// Period dependence +// ----------------- +// The coefficients used for the corrections depend on the run number. +// The mapping from run -> (coefficients) is defined in Config.yaml. +// Config.yaml provides one or more "periods" blocks, each with: +// - run_ranges: list of [min,max] run-number ranges +// - FD coefficients +// - CD coefficients +// + +#include "iguana/algorithms/Algorithm.h" +#include "iguana/algorithms/TypeDefs.h" + +#include +#include +#include +#include +#include + +namespace iguana::clas12::rga { + + class ProtonEnergyLossCorrection : public Algorithm { + DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection) + + public: + // @action_function{scalar transformer} + // + // Inputs: + // pid : PDG ID (correct only protons: 2212) + // status : REC::Particle status code (used to identify detector) + // run : run number (selects the period and coefficients) + // px,py,pz : momentum components (GeV) + // + // Output: + // corrected (px,py,pz). + std::tuple Transform( + int const pid, + int const status, + int const run, + vector_element_t const px, + vector_element_t const py, + vector_element_t const pz) const; + + private: + void ConfigHook() override; + void StartHook(hipo::banklist& banks) override; + bool RunHook(hipo::banklist& banks) const override; + void StopHook() override; + + // Cached bank indices (resolved once in StartHook) + hipo::banklist::size_type m_b_rec_particle{}; + hipo::banklist::size_type m_b_run_config{}; + + // Simple polynomial container: c0 + c1*x + c2*x^2 + ... + struct Poly { + std::vector c; + }; + + // Coefficients for a single detector region (FD or CD). + // For each of p, theta, phi store (A,B,C) polynomials in theta: A(theta), B(theta), C(theta) + // Then the correction formula uses those A,B,C values. + struct RegionCoeffs { + Poly A_p; + Poly B_p; + Poly C_p; + + Poly A_theta; + Poly B_theta; + Poly C_theta; + + Poly A_phi; + Poly B_phi; + Poly C_phi; + }; + + // Full period definition: + // - run_ranges define which runs map to this period + // - fd / cd hold coefficients for Forward Detector and Central Detector + struct PeriodDef { + std::vector> run_ranges; + RegionCoeffs fd; + RegionCoeffs cd; + }; + + // Map period key -> period definition (loaded from YAML) + std::map m_periods{}; + + // Helper functions + static bool IsFD(int status); + static bool IsCD(int status); + + static double EvalPoly(Poly const& p, double x); + + static double Pmag(double px, double py, double pz); + static double ThetaDeg(double px, double py, double pz); + static double PhiDeg(double px, double py); + + // Convert (p, theta_deg, phi_deg) back to (px,py,pz). + static void SphericalToCartesian(double p, double theta_deg, double phi_deg, + double& px, double& py, double& pz); + + // Given a run number, return the matching run period (e.g. Fa18 Inb) or nullptr. + PeriodDef const* FindPeriod(int run) const; + }; + +} \ No newline at end of file diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml new file mode 100644 index 00000000..4e4a7432 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml @@ -0,0 +1,127 @@ +clas12: + ProtonEnergyLossCorrection: + + periods: + + rga_sp18_inb: + run_ranges: + - [3031, 3087] + - [3306, 3817] + - [4003, 4325] + FD: + A_p: [0.0146275, -0.00124929, 3.64154e-05] + B_p: [-0.00743169, 0.000458648, -6.45703e-06] + C_p: [0.0175282, -0.00128554, 3.5249e-05] + A_theta: [0.0] + B_theta: [0.0] + C_theta: [0.0] + A_phi: [0.0] + B_phi: [0.0] + C_phi: [0.0] + CD: + A_p: [-0.229055, 0.00924571, -9.09927e-05] + B_p: [0.371002, -0.0146818, 0.000146548] + C_p: [-0.174565, 0.00680452, -6.9e-05] + A_theta: [0.0] + B_theta: [0.0] + C_theta: [0.0] + A_phi: [0.0] + B_phi: [0.0] + C_phi: [0.0] + + rga_sp18_out: + run_ranges: + - [3103, 3293] + - [3820, 3987] + FD: + A_p: [0.00523188, -9.43614e-05] + B_p: [-0.00887291, 0.000759277] + C_p: [0.0] + A_theta: [0.0] + B_theta: [0.0] + C_theta: [0.0] + A_phi: [0.0] + B_phi: [0.0] + C_phi: [0.0] + CD: + A_p: [-0.204359, 0.00857339, -8.79867e-05] + B_p: [0.402543, -0.0168624, 0.000178539] + C_p: [-0.217865, 0.00908787, -9.77617e-05] + A_theta: [0.0] + B_theta: [0.0] + C_theta: [0.0] + A_phi: [0.0] + B_phi: [0.0] + C_phi: [0.0] + + rga_fa18_inb: + run_ranges: + - [4763, 5419] + FD: + A_p: [0.0099626, -0.0002414, -2.0e-06] + B_p: [-0.01428267, 0.00042833, 1.081e-05] + C_p: [0.01197102, -0.00055673, 7.85e-06] + A_theta: [0.0683831, -0.0083821, 0.0001670] + B_theta: [-0.15834256, 0.02630760, -0.00064126] + C_theta: [0.11587509, -0.01679559, 0.00038915] + A_phi: [0.0416510, -0.0064212, 0.0000622] + B_phi: [0.28414191, -0.00047647, 0.00010357] + C_phi: [-0.25690893, 0.00886707, -0.00016081] + CD: + A_p: [-0.2383991, 0.0124992, -0.0001646] + B_p: [0.60123885, -0.03128464, 0.00041314] + C_p: [-0.44080146, 0.02209857, -0.00028224] + A_theta: [0.1000890, -0.0039222, 0.0000359] + B_theta: [-0.0130680, 0.0004545, -0.0000026] + C_theta: [0.0] + A_phi: [0.0776934, -0.0059632, 0.0000749] + B_phi: [-0.31582008, 0.01649220, -0.00018505] + C_phi: [0.10909746, -0.00530642, 0.00005627] + + rga_fa18_out: + run_ranges: + - [5423, 5666] + FD: + A_p: [0.0135790, -0.0005303] + B_p: [-0.02165929, 0.00121123] + C_p: [0.0] + A_theta: [-0.3715486, 0.0272810, -0.0006278, 0.0000040] + B_theta: [2.00009939, -0.20781779, 0.00721092, -0.00008343] + C_theta: [0.0] + A_phi: [-0.9701486, 0.1213124, -0.0049215, 0.0000640] + B_phi: [2.85034691, -0.34405076, 0.01347377, -0.00016663] + C_phi: [0.0] + CD: + A_p: [-0.1927861, 0.0099546, -0.0001299] + B_p: [0.44307822, -0.02309469, 0.00030784] + C_p: [-0.32938000, 0.01648659, -0.00021181] + A_theta: [0.0581473, -0.0021818, 0.0000181] + B_theta: [0.00915748, -0.00040748, 0.00000562] + C_theta: [0.0] + A_phi: [-0.0733814, 0.0010335, -0.0000044] + B_phi: [-0.06127800, 0.00492239, -0.00005683] + C_phi: [0.02586507, -0.00160176, 0.00001642] + + rga_sp19_inb: + run_ranges: + - [6616, 6783] + FD: + A_p: [0.0095205, -0.0001914, -0.0000031] + B_p: [-0.01365658, 0.00036322, 0.00001217] + C_p: [0.01175256, -0.00053407, 0.00000742] + A_theta: [0.0723069, -0.0085078, 0.0001702] + B_theta: [-0.16048057, 0.02561073, -0.00062158] + C_theta: [0.10954630, -0.01566605, 0.00036132] + A_phi: [0.0486986, -0.0067579, 0.0000638] + B_phi: [0.26803189, 0.00016245, 0.00010433] + C_phi: [-0.24522460, 0.00826646, -0.00015640] + CD: + A_p: [-0.2716918, 0.0142491, -0.0001862] + B_p: [0.65945101, -0.03431360, 0.00045036] + C_p: [-0.46602726, 0.02335623, -0.00029720] + A_theta: [0.2550377, -0.0107983, 0.0001116] + B_theta: [-0.14022533, 0.00596067, -0.00006172] + C_theta: [0.0] + A_phi: [-0.5459156, 0.0219868, -0.0002349] + B_phi: [0.74223687, -0.03037065, 0.00032761] + C_phi: [-0.29798258, 0.01246744, -0.00013525] \ No newline at end of file diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc new file mode 100644 index 00000000..a0422e48 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc @@ -0,0 +1,178 @@ +#include "Validator.h" + +#include +#include +#include +#include + +namespace iguana::clas12::rga { + + REGISTER_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator); + + // Compute theta (deg) from momentum components. + double ProtonEnergyLossCorrectionValidator::ThetaDegFromPxPyPz(double px, double py, double pz) { + double pt = std::sqrt(px * px + py * py); + double th = std::atan2(pt, pz); + return th * (180.0 / M_PI); + } + + // Convert theta (deg) to a bin index. + // + // Returns: + // -1 if theta is outside [kThetaMinDeg, kThetaMaxDeg] + // otherwise an index in [0, kNBins-1]. + int ProtonEnergyLossCorrectionValidator::ThetaBinIndex(double theta_deg) { + if(theta_deg < kThetaMinDeg) { + return -1; + } + if(theta_deg > kThetaMaxDeg) { + return -1; + } + if(theta_deg >= kThetaMaxDeg) { + return kNBins - 1; + } + + int idx = (int)((theta_deg - kThetaMinDeg) / kThetaStepDeg); + if(idx < 0 || idx >= kNBins) { + return -1; + } + return idx; + } + + void ProtonEnergyLossCorrectionValidator::StartHook(hipo::banklist& banks) { + // Cache bank indices + m_b_particle = GetBankIndex(banks, "REC::Particle"); + m_b_config = GetBankIndex(banks, "RUN::config"); + + // Build and start an AlgorithmSequence + m_algo_seq = std::make_unique(); + m_algo_seq->Add("clas12::rga::ProtonEnergyLossCorrection"); + m_algo_seq->Start(banks); + + // Reset counters/accumulators. + for(auto& b : m_bins) { + b.n = 0; + b.sum_p_before = 0.0; + b.sum_p_after = 0.0; + } + m_total_protons_in_range = 0; + m_total_protons_all = 0; + } + + bool ProtonEnergyLossCorrectionValidator::RunHook(hipo::banklist& banks) const { + auto& particle = GetBank(banks, m_b_particle, "REC::Particle"); + auto& config = GetBank(banks, m_b_config, "RUN::config"); + (void)config; // run number is not needed for summary of \Delta p. + + // record before values, run the algorithm, compute after values. + // store the particle row index and p_before to we can match the same + // rows after the algorithm modifies the bank in place. + struct Entry { + int row = -1; + int bin = -1; + double p_before = 0.0; + }; + + std::vector entries; + entries.reserve((size_t)particle.getRows()); + + int const nrows = particle.getRows(); + + // snapshot before for all protons. + for(int i = 0; i < nrows; ++i) { + int pid = particle.getInt("pid", i); + if(pid != 2212) { + continue; + } + + double px = (double)particle.getFloat("px", i); + double py = (double)particle.getFloat("py", i); + double pz = (double)particle.getFloat("pz", i); + + double p = std::sqrt(px * px + py * py + pz * pz); + double theta_deg = ThetaDegFromPxPyPz(px, py, pz); + int bin = ThetaBinIndex(theta_deg); + + Entry e; + e.row = i; + e.bin = bin; + e.p_before = p; + entries.push_back(e); + } + + // Run algorithm under test + m_algo_seq->Run(banks); + + // Accumulate after values. + std::scoped_lock lock(m_mutex); + + for(auto const& e : entries) { + m_total_protons_all++; + + if(e.bin < 0) { + continue; + } + m_total_protons_in_range++; + + int i = e.row; + + double px = (double)particle.getFloat("px", i); + double py = (double)particle.getFloat("py", i); + double pz = (double)particle.getFloat("pz", i); + + double p_after = std::sqrt(px * px + py * py + pz * pz); + + auto& B = m_bins[(size_t)e.bin]; + B.n++; + B.sum_p_before += e.p_before; + B.sum_p_after += p_after; + } + + return true; + } + + void ProtonEnergyLossCorrectionValidator::StopHook() { + // human-readable summary table + std::cout << "\n"; + std::cout << "ProtonEnergyLossCorrectionValidator summary\n"; + std::cout << " theta bins: " << kThetaMinDeg << " to " << kThetaMaxDeg + << " (deg) in steps of " << kThetaStepDeg << " (deg)\n"; + std::cout << " total protons (all theta): " << m_total_protons_all << "\n"; + std::cout << " total protons in theta range: " << m_total_protons_in_range << "\n\n"; + + std::cout << std::fixed << std::setprecision(6); + std::cout << " BinRangeTheta(deg) N (GeV) (GeV) (GeV)\n"; + std::cout << " -------------------------------------------------------------------------------------\n"; + + for(int ib = 0; ib < kNBins; ++ib) { + double lo = kThetaMinDeg + (double)ib * kThetaStepDeg; + double hi = lo + kThetaStepDeg; + + auto const& B = m_bins[(size_t)ib]; + + if(B.n <= 0) { + std::cout << " [" << std::setw(4) << lo << "," << std::setw(4) << hi << ")" + << " " << std::setw(8) << 0 + << " " << std::setw(14) << 0.0 + << " " << std::setw(13) << 0.0 + << " " << std::setw(12) << 0.0 + << "\n"; + continue; + } + + double mean_before = B.sum_p_before / (double)B.n; + double mean_after = B.sum_p_after / (double)B.n; + double delta = mean_after - mean_before; + + std::cout << " [" << std::setw(4) << lo << "," << std::setw(4) << hi << ")" + << " " << std::setw(8) << B.n + << " " << std::setw(14) << mean_before + << " " << std::setw(13) << mean_after + << " " << std::setw(12) << delta + << "\n"; + } + + std::cout << "\n"; + } + +} \ No newline at end of file diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h new file mode 100644 index 00000000..a7cd1ed6 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h @@ -0,0 +1,72 @@ +#pragma once + +// ProtonEnergyLossCorrectionValidator +// +// For each processed event: +// 1) Find all protons (pid==2212) in REC::Particle. +// 2) Compute their theta before correction, place them into theta bins. +// 3) Run the ProtonEnergyLossCorrection algorithm in place. +// 4) Recompute momentum magnitude p after correction. +// 5) Accumulate mean p_before and mean p_after per theta bin. +// +// Output +// ------ +// In StopHook(), print a table: +// theta bin range, N, , , +// + +#include "iguana/algorithms/Validator.h" +#include "iguana/algorithms/AlgorithmSequence.h" + +#include +#include +#include + +namespace iguana::clas12::rga { + + class ProtonEnergyLossCorrectionValidator : public Validator { + DEFINE_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator, clas12::rga::ProtonEnergyLossCorrectionValidator) + + private: + void StartHook(hipo::banklist& banks) override; + bool RunHook(hipo::banklist& banks) const override; + void StopHook() override; + + private: + // Cached bank indices. + hipo::banklist::size_type m_b_particle{}; + hipo::banklist::size_type m_b_config{}; + + std::unique_ptr m_algo_seq; + + // Theta binning configuration (degrees). + // Bins: [5,10), [10,15), ... + static constexpr double kThetaMinDeg = 5.0; + static constexpr double kThetaMaxDeg = 70.0; + static constexpr double kThetaStepDeg = 5.0; + static constexpr int kNBins = (int)((kThetaMaxDeg - kThetaMinDeg) / kThetaStepDeg); + + // Simple per-bin accumulators. + struct BinAccum { + long long n = 0; + double sum_p_before = 0.0; + double sum_p_after = 0.0; + }; + + // accumulate over events -> these must be mutable. + mutable std::array m_bins{}; + mutable long long m_total_protons_in_range = 0; + mutable long long m_total_protons_all = 0; + + // Protect shared accumulation in multi-threaded contexts. + mutable std::mutex m_mutex{}; + + private: + // Map theta (deg) to bin index. + static int ThetaBinIndex(double theta_deg); + + // Compute theta (deg) from (px,py,pz). + static double ThetaDegFromPxPyPz(double px, double py, double pz); + }; + +} \ No newline at end of file diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 60a53f0a..84efb862 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -82,6 +82,14 @@ algo_dict = [ 'has_validator': false, 'test_args': {'banks': [ 'RECFT::Particle' ]}, }, + { + 'name': 'clas12::rga::ProtonEnergyLossCorrection', + 'has_validator': true, + 'validator_needs_ROOT': false, + 'test_args': { + 'banks': [ 'REC::Particle', 'RUN::config' ], + }, + }, { 'name': 'clas12::rga::MomentumCorrection', 'has_config': false,