From f0e0a8c40833143d3a7eefc680ac15b4498f8f6e Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Mon, 9 Feb 2026 18:52:12 -0500 Subject: [PATCH 01/13] Add clas12::rga::ProtonEnergyLossCorrection algorithm skeleton --- .../ProtonEnergyLossCorrection/Action.yaml | 34 ++ .../ProtonEnergyLossCorrection/Algorithm.cc | 313 ++++++++++++++++++ .../ProtonEnergyLossCorrection/Algorithm.h | 74 +++++ .../ProtonEnergyLossCorrection/Config.yaml | 0 src/iguana/algorithms/meson.build | 7 + 5 files changed, 428 insertions(+) create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Action.yaml create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml 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 000000000..f4c2bc690 --- /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 000000000..93a7b81d6 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -0,0 +1,313 @@ +#include "Algorithm.h" + +#include +#include + +#include + +namespace iguana::clas12 { + + REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); + + // -------------------------- + // Small math helpers + // -------------------------- + 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); + } + + double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) { + return std::sqrt(px * px + py * py + pz * pz); + } + + 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); + } + + // This reproduces your Java phi_calculation(x,y) EXACTLY: + // phi = toDegrees(atan2(x, y)); + // phi = phi - 90; + // if (phi < 0) phi = 360 + phi; + // phi = 360 - phi; + // Returns phi in [0,360). + double ProtonEnergyLossCorrection::PhiDegLikeJava(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; + // keep in [0,360) + while (phi >= 360.0) { + phi -= 360.0; + } + while (phi < 0.0) { + phi += 360.0; + } + return phi; + } + + 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); + } + + 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; + } + + ProtonEnergyLossCorrection::PeriodDef const* ProtonEnergyLossCorrection::FindPeriod(int run) const { + for (auto const& kv : periods_) { + auto const& def = kv.second; + for (auto const& rr : def.run_ranges) { + if (run >= rr.first && run <= rr.second) { + return &def; + } + } + } + return nullptr; + } + + // -------------------------- + // Hooks + // -------------------------- + void ProtonEnergyLossCorrection::ConfigHook() { + periods_.clear(); + + // The Algorithm base class provides config path machinery; most iguana algos + // load YAML from the standard algorithm config directory. + // We use the standard YAML-CPP loader on the resolved config file. + // + // If your branch uses a helper like GetConfigYAML(), swap this accordingly. + // For now we use the conventional path from Algorithm::GetConfigPath(). + std::string cfg_path = GetConfigPath(); + + YAML::Node root = YAML::LoadFile(cfg_path); + + if (!root["clas12"] || !root["clas12"]["ProtonEnergyLossCorrection"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node clas12:ProtonEnergyLossCorrection"); + } + + YAML::Node node = root["clas12"]["ProtonEnergyLossCorrection"]; + if (!node["periods"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node periods"); + } + + auto periods_node = node["periods"]; + for (auto it = periods_node.begin(); it != periods_node.end(); ++it) { + std::string key = it->first.as(); + YAML::Node pnode = it->second; + + PeriodDef def; + + if (!pnode["run_ranges"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: missing run_ranges for period " + key); + } + for (auto rr : pnode["run_ranges"]) { + if (!rr.IsSequence() || rr.size() != 2) { + throw std::runtime_error("ProtonEnergyLossCorrection: run_ranges entry must be [min,max] for period " + key); + } + def.run_ranges.push_back({rr[0].as(), rr[1].as()}); + } + + auto load_poly = [](YAML::Node const& n, char const* name) -> Poly { + Poly p; + if (!n[name]) { + throw std::runtime_error(std::string("ProtonEnergyLossCorrection: missing coeff node ") + name); + } + YAML::Node a = n[name]; + if (!a.IsSequence()) { + throw std::runtime_error(std::string("ProtonEnergyLossCorrection: coeff node not sequence: ") + name); + } + for (auto v : a) { + p.c.push_back(v.as()); + } + return p; + }; + + 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; + }; + + if (!pnode["FD"] || !pnode["CD"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: missing FD or CD node for period " + key); + } + + def.fd = load_region(pnode["FD"]); + def.cd = load_region(pnode["CD"]); + + periods_[key] = def; + } + } + + void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) { + b_rec_particle = GetBankIndex(banks, "REC::Particle"); + b_run_config = GetBankIndex(banks, "RUN::config"); + } + + bool ProtonEnergyLossCorrection::RunHook(hipo::banklist& banks) const { + auto& rec = GetBank(banks, b_rec_particle, "REC::Particle"); + auto& run = GetBank(banks, b_run_config, "RUN::config"); + + int runnum = run.getInt("run", 0); + + for (auto const& row : rec.getRowList()) { + int pid = rec.getInt("pid", row); + if (pid != 2212) { + continue; + } + + int status = rec.getInt("status", row); + if (!IsFD(status) && !IsCD(status)) { + continue; + } + + double px = rec.getFloat("px", row); + double py = rec.getFloat("py", row); + double pz = rec.getFloat("pz", row); + + auto [px_new, py_new, pz_new] = Transform(pid, status, runnum, px, py, pz); + + 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 + // -------------------------- + 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 + if (pid != 2212) { + return {px_in, py_in, pz_in}; + } + + bool is_fd = IsFD(status); + bool is_cd = IsCD(status); + if (!is_fd && !is_cd) { + return {px_in, py_in, pz_in}; + } + + PeriodDef const* period = FindPeriod(run); + if (!period) { + return {px_in, py_in, pz_in}; + } + + double px = px_in; + double py = py_in; + double pz = 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); + double phi = PhiDegLikeJava(px, py); + + RegionCoeffs const& c = is_fd ? period->fd : period->cd; + + // Build A/B/C from theta polynomials + double A_p = EvalPoly(c.A_p, theta); + double B_p = EvalPoly(c.B_p, theta); + double C_p = EvalPoly(c.C_p, theta); + + double A_theta = EvalPoly(c.A_theta, theta); + double B_theta = EvalPoly(c.B_theta, theta); + double C_theta = EvalPoly(c.C_theta, theta); + + double A_phi = EvalPoly(c.A_phi, theta); + double B_phi = EvalPoly(c.B_phi, theta); + double C_phi = EvalPoly(c.C_phi, theta); + + // Apply corrections using your exact functional forms + 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); + + // Guard theta,phi divisions + 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 if (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 convention: keep in [0,360) since that is what your Java uses + while (phi_new >= 360.0) { + phi_new -= 360.0; + } + while (phi_new < 0.0) { + phi_new += 360.0; + } + + // Convert back to Cartesian using spherical conversion + double px_out = 0.0, py_out = 0.0, pz_out = 0.0; + SphericalToCartesian(p_new, theta_new, phi_new, px_out, py_out, pz_out); + + return {px_out, py_out, pz_out}; + } + +} // namespace iguana::clas12 \ 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 000000000..a545a2372 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -0,0 +1,74 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" +#include "iguana/algorithms/TypeDefs.h" + +#include +#include +#include +#include + +namespace iguana::clas12 { + + class ProtonEnergyLossCorrection : public Algorithm { + + DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::ProtonEnergyLossCorrection) + + public: + /// @action_function{scalar transformer} + /// Transform a single particle row (px,py,pz) if it is a proton in FD or CD and within an RGA run range. + /// Returns 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: + // Hook API (Start/Run/Stop are final in base class) + void ConfigHook() override; + void StartHook(hipo::banklist& banks) override; + bool RunHook(hipo::banklist& banks) const override; + void StopHook() override; + + // Bank indices + hipo::banklist::size_type b_rec_particle; + hipo::banklist::size_type b_run_config; + + // Internal config structures + struct Poly { + std::vector c; // c0 + c1*x + c2*x^2 + ... + }; + + struct RegionCoeffs { + Poly A_p, B_p, C_p; + Poly A_theta, B_theta, C_theta; + Poly A_phi, B_phi, C_phi; + }; + + struct PeriodDef { + std::vector> run_ranges; + RegionCoeffs fd; + RegionCoeffs cd; + }; + + // period key -> period def + std::map periods_; + + // Helpers + 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 PhiDegLikeJava(double px, double py); + static void SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& pz); + + // Find matching period for run, returns nullptr if none + PeriodDef const* FindPeriod(int run) const; + + }; + +} // namespace iguana::clas12 \ 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 000000000..e69de29bb diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 60a53f0a4..bcccfb61a 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -82,6 +82,13 @@ algo_dict = [ 'has_validator': false, 'test_args': {'banks': [ 'RECFT::Particle' ]}, }, + { + 'name': 'clas12::rga::ProtonEnergyLossCorrection', + 'has_validator': false, + 'test_args': { + 'banks': [ 'REC::Particle', 'RUN::config' ], + }, + }, { 'name': 'clas12::rga::MomentumCorrection', 'has_config': false, From c00e87b8dfecb5399f8854d31f903f502fe73b3d Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 10:08:13 -0500 Subject: [PATCH 02/13] updated GetConfigPath() call in Algorithm.cc --- .../clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index 93a7b81d6..d98a02f91 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -41,7 +41,6 @@ namespace iguana::clas12 { return (180.0 / M_PI) * std::acos(c); } - // This reproduces your Java phi_calculation(x,y) EXACTLY: // phi = toDegrees(atan2(x, y)); // phi = phi - 90; // if (phi < 0) phi = 360 + phi; @@ -100,13 +99,9 @@ namespace iguana::clas12 { void ProtonEnergyLossCorrection::ConfigHook() { periods_.clear(); - // The Algorithm base class provides config path machinery; most iguana algos - // load YAML from the standard algorithm config directory. - // We use the standard YAML-CPP loader on the resolved config file. - // // If your branch uses a helper like GetConfigYAML(), swap this accordingly. // For now we use the conventional path from Algorithm::GetConfigPath(). - std::string cfg_path = GetConfigPath(); + std::string cfg_path = GetConfig()->FindFile("Config.yaml"); YAML::Node root = YAML::LoadFile(cfg_path); From e86b602feec6a872c8d00faafb4b9c2a34b1b128 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 10:34:49 -0500 Subject: [PATCH 03/13] changed registration name of algorithm to clas12::rga::ProtonEnergyLossCorrection --- .../clas12/rga/ProtonEnergyLossCorrection/Algorithm.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h index a545a2372..92ff2c68b 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -12,12 +12,12 @@ namespace iguana::clas12 { class ProtonEnergyLossCorrection : public Algorithm { - DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::ProtonEnergyLossCorrection) + DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection) public: - /// @action_function{scalar transformer} - /// Transform a single particle row (px,py,pz) if it is a proton in FD or CD and within an RGA run range. - /// Returns corrected (px,py,pz). + // @action_function{scalar transformer} + // Transform a single particle row (px,py,pz) if it is an RGA proton + // Returns corrected (px,py,pz). std::tuple Transform( int const pid, int const status, From 63573fa2b15402c709efbb2cad74a8a24f0890c6 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 10:41:05 -0500 Subject: [PATCH 04/13] changed registration name of algorithm to clas12::rga::ProtonEnergyLossCorrection --- .../ProtonEnergyLossCorrection/Algorithm.h | 115 +++++++++--------- 1 file changed, 57 insertions(+), 58 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h index 92ff2c68b..f82106e5d 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -5,70 +5,69 @@ #include #include +#include #include #include -namespace iguana::clas12 { +namespace iguana::clas12::rga { class ProtonEnergyLossCorrection : public Algorithm { - DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection) - public: - // @action_function{scalar transformer} - // Transform a single particle row (px,py,pz) if it is an RGA proton - // Returns 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: - // Hook API (Start/Run/Stop are final in base class) - void ConfigHook() override; - void StartHook(hipo::banklist& banks) override; - bool RunHook(hipo::banklist& banks) const override; - void StopHook() override; - - // Bank indices - hipo::banklist::size_type b_rec_particle; - hipo::banklist::size_type b_run_config; - - // Internal config structures - struct Poly { - std::vector c; // c0 + c1*x + c2*x^2 + ... - }; - - struct RegionCoeffs { - Poly A_p, B_p, C_p; - Poly A_theta, B_theta, C_theta; - Poly A_phi, B_phi, C_phi; - }; - - struct PeriodDef { - std::vector> run_ranges; - RegionCoeffs fd; - RegionCoeffs cd; - }; - - // period key -> period def - std::map periods_; - - // Helpers - 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 PhiDegLikeJava(double px, double py); - static void SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& pz); - - // Find matching period for run, returns nullptr if none - PeriodDef const* FindPeriod(int run) const; - + public: + // @action_function{scalar transformer} + // Transform a single particle row (px,py,pz) if it is an RGA proton. + // Returns 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: + // Hook API (Start/Run/Stop are final in base class) + void ConfigHook() override; + void StartHook(hipo::banklist& banks) override; + bool RunHook(hipo::banklist& banks) const override; + void StopHook() override; + + // Bank indices + hipo::banklist::size_type b_rec_particle; + hipo::banklist::size_type b_run_config; + + // Internal config structures + struct Poly { + std::vector c; // c0 + c1*x + c2*x^2 + ... + }; + + struct RegionCoeffs { + Poly A_p, B_p, C_p; + Poly A_theta, B_theta, C_theta; + Poly A_phi, B_phi, C_phi; + }; + + struct PeriodDef { + std::vector> run_ranges; + RegionCoeffs fd; + RegionCoeffs cd; + }; + + // period key -> period def + std::map periods_; + + // Helpers + 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 PhiDegLikeJava(double px, double py); + static void SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& pz); + + // Find matching period for run, returns nullptr if none + PeriodDef const* FindPeriod(int run) const; }; -} // namespace iguana::clas12 \ No newline at end of file +} // namespace iguana::clas12::rga \ No newline at end of file From 3a858bfae565e5077ddedac6b8bc5374df7559f4 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 10:45:32 -0500 Subject: [PATCH 05/13] changed registration name of algorithm to clas12::rga::ProtonEnergyLossCorrection --- .../ProtonEnergyLossCorrection/Algorithm.cc | 147 ++++++++++-------- 1 file changed, 80 insertions(+), 67 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index d98a02f91..2d878c65b 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -5,37 +5,41 @@ #include -namespace iguana::clas12 { +namespace iguana::clas12::rga { REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); // -------------------------- // Small math helpers // -------------------------- - bool ProtonEnergyLossCorrection::IsFD(int status) { + bool ProtonEnergyLossCorrection::IsFD(int status) + { int s = std::abs(status); return (s >= 2000 && s < 4000); } - bool ProtonEnergyLossCorrection::IsCD(int status) { + bool ProtonEnergyLossCorrection::IsCD(int status) + { int s = std::abs(status); return (s >= 4000 && s < 5000); } - double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) { + double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) + { return std::sqrt(px * px + py * py + pz * pz); } - double ProtonEnergyLossCorrection::ThetaDeg(double px, double py, double pz) { + double ProtonEnergyLossCorrection::ThetaDeg(double px, double py, double pz) + { double r = Pmag(px, py, pz); - if (r <= 0.0) { + if(r <= 0.0) { return 0.0; } double c = pz / r; - if (c > 1.0) { + if(c > 1.0) { c = 1.0; } - if (c < -1.0) { + if(c < -1.0) { c = -1.0; } return (180.0 / M_PI) * std::acos(c); @@ -46,46 +50,50 @@ namespace iguana::clas12 { // if (phi < 0) phi = 360 + phi; // phi = 360 - phi; // Returns phi in [0,360). - double ProtonEnergyLossCorrection::PhiDegLikeJava(double px, double py) { + double ProtonEnergyLossCorrection::PhiDegLikeJava(double px, double py) + { double phi = (180.0 / M_PI) * std::atan2(px, py); - phi = phi - 90.0; - if (phi < 0.0) { + phi = phi - 90.0; + if(phi < 0.0) { phi = 360.0 + phi; } phi = 360.0 - phi; // keep in [0,360) - while (phi >= 360.0) { + while(phi >= 360.0) { phi -= 360.0; } - while (phi < 0.0) { + while(phi < 0.0) { phi += 360.0; } return phi; } - void ProtonEnergyLossCorrection::SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& 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); + 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); } - double ProtonEnergyLossCorrection::EvalPoly(Poly const& p, double x) { + double ProtonEnergyLossCorrection::EvalPoly(Poly const& p, double x) + { double sum = 0.0; - double xn = 1.0; - for (double c : p.c) { + double xn = 1.0; + for(double c : p.c) { sum += c * xn; xn *= x; } return sum; } - ProtonEnergyLossCorrection::PeriodDef const* ProtonEnergyLossCorrection::FindPeriod(int run) const { - for (auto const& kv : periods_) { + ProtonEnergyLossCorrection::PeriodDef const* ProtonEnergyLossCorrection::FindPeriod(int run) const + { + for(auto const& kv : periods_) { auto const& def = kv.second; - for (auto const& rr : def.run_ranges) { - if (run >= rr.first && run <= rr.second) { + for(auto const& rr : def.run_ranges) { + if(run >= rr.first && run <= rr.second) { return &def; } } @@ -96,36 +104,38 @@ namespace iguana::clas12 { // -------------------------- // Hooks // -------------------------- - void ProtonEnergyLossCorrection::ConfigHook() { + void ProtonEnergyLossCorrection::ConfigHook() + { periods_.clear(); - // If your branch uses a helper like GetConfigYAML(), swap this accordingly. - // For now we use the conventional path from Algorithm::GetConfigPath(). + // This uses the algorithm's YAMLReader search paths: + // - ./, any AddDirectory paths, then install prefix. + // It will locate the installed Config.yaml for this algorithm unless overridden. std::string cfg_path = GetConfig()->FindFile("Config.yaml"); YAML::Node root = YAML::LoadFile(cfg_path); - if (!root["clas12"] || !root["clas12"]["ProtonEnergyLossCorrection"]) { + if(!root["clas12"] || !root["clas12"]["ProtonEnergyLossCorrection"]) { throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node clas12:ProtonEnergyLossCorrection"); } YAML::Node node = root["clas12"]["ProtonEnergyLossCorrection"]; - if (!node["periods"]) { + if(!node["periods"]) { throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node periods"); } auto periods_node = node["periods"]; - for (auto it = periods_node.begin(); it != periods_node.end(); ++it) { + for(auto it = periods_node.begin(); it != periods_node.end(); ++it) { std::string key = it->first.as(); YAML::Node pnode = it->second; PeriodDef def; - if (!pnode["run_ranges"]) { + if(!pnode["run_ranges"]) { throw std::runtime_error("ProtonEnergyLossCorrection: missing run_ranges for period " + key); } - for (auto rr : pnode["run_ranges"]) { - if (!rr.IsSequence() || rr.size() != 2) { + for(auto rr : pnode["run_ranges"]) { + if(!rr.IsSequence() || rr.size() != 2) { throw std::runtime_error("ProtonEnergyLossCorrection: run_ranges entry must be [min,max] for period " + key); } def.run_ranges.push_back({rr[0].as(), rr[1].as()}); @@ -133,14 +143,14 @@ namespace iguana::clas12 { auto load_poly = [](YAML::Node const& n, char const* name) -> Poly { Poly p; - if (!n[name]) { + if(!n[name]) { throw std::runtime_error(std::string("ProtonEnergyLossCorrection: missing coeff node ") + name); } YAML::Node a = n[name]; - if (!a.IsSequence()) { + if(!a.IsSequence()) { throw std::runtime_error(std::string("ProtonEnergyLossCorrection: coeff node not sequence: ") + name); } - for (auto v : a) { + for(auto v : a) { p.c.push_back(v.as()); } return p; @@ -148,19 +158,19 @@ namespace iguana::clas12 { 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_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"); + 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; }; - if (!pnode["FD"] || !pnode["CD"]) { + if(!pnode["FD"] || !pnode["CD"]) { throw std::runtime_error("ProtonEnergyLossCorrection: missing FD or CD node for period " + key); } @@ -171,25 +181,27 @@ namespace iguana::clas12 { } } - void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) { + void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) + { b_rec_particle = GetBankIndex(banks, "REC::Particle"); - b_run_config = GetBankIndex(banks, "RUN::config"); + b_run_config = GetBankIndex(banks, "RUN::config"); } - bool ProtonEnergyLossCorrection::RunHook(hipo::banklist& banks) const { + bool ProtonEnergyLossCorrection::RunHook(hipo::banklist& banks) const + { auto& rec = GetBank(banks, b_rec_particle, "REC::Particle"); auto& run = GetBank(banks, b_run_config, "RUN::config"); int runnum = run.getInt("run", 0); - for (auto const& row : rec.getRowList()) { + for(auto const& row : rec.getRowList()) { int pid = rec.getInt("pid", row); - if (pid != 2212) { + if(pid != 2212) { continue; } int status = rec.getInt("status", row); - if (!IsFD(status) && !IsCD(status)) { + if(!IsFD(status) && !IsCD(status)) { continue; } @@ -207,7 +219,8 @@ namespace iguana::clas12 { return true; } - void ProtonEnergyLossCorrection::StopHook() { + void ProtonEnergyLossCorrection::StopHook() + { } // -------------------------- @@ -222,18 +235,18 @@ namespace iguana::clas12 { vector_element_t const pz_in) const { // Only protons - if (pid != 2212) { + if(pid != 2212) { return {px_in, py_in, pz_in}; } bool is_fd = IsFD(status); bool is_cd = IsCD(status); - if (!is_fd && !is_cd) { + if(!is_fd && !is_cd) { return {px_in, py_in, pz_in}; } PeriodDef const* period = FindPeriod(run); - if (!period) { + if(!period) { return {px_in, py_in, pz_in}; } @@ -242,12 +255,12 @@ namespace iguana::clas12 { double pz = pz_in; double p = Pmag(px, py, pz); - if (p <= 0.0) { + if(p <= 0.0) { return {px_in, py_in, pz_in}; } double theta = ThetaDeg(px, py, pz); - double phi = PhiDegLikeJava(px, py); + double phi = PhiDegLikeJava(px, py); RegionCoeffs const& c = is_fd ? period->fd : period->cd; @@ -265,36 +278,36 @@ namespace iguana::clas12 { double C_phi = EvalPoly(c.C_phi, theta); // Apply corrections using your exact functional forms - double p_new = p; + double p_new = p; double theta_new = theta; - double phi_new = phi; + double phi_new = phi; - if (is_fd) { + if(is_fd) { p_new += A_p + B_p / p + C_p / (p * p); // Guard theta,phi divisions - if (theta != 0.0) { + if(theta != 0.0) { theta_new += A_theta + B_theta / theta + C_theta / (theta * theta); } - if (phi != 0.0) { + if(phi != 0.0) { phi_new += A_phi + B_phi / phi + C_phi / (phi * phi); } - } else if (is_cd) { + } else if(is_cd) { p_new += A_p + B_p * p + C_p * p * p; - if (theta != 0.0) { + if(theta != 0.0) { theta_new += A_theta + B_theta / theta + C_theta / (theta * theta); } - if (phi != 0.0) { + if(phi != 0.0) { phi_new += A_phi + B_phi / phi + C_phi / (phi * phi); } } // phi convention: keep in [0,360) since that is what your Java uses - while (phi_new >= 360.0) { + while(phi_new >= 360.0) { phi_new -= 360.0; } - while (phi_new < 0.0) { + while(phi_new < 0.0) { phi_new += 360.0; } @@ -305,4 +318,4 @@ namespace iguana::clas12 { return {px_out, py_out, pz_out}; } -} // namespace iguana::clas12 \ No newline at end of file +} // namespace iguana::clas12::rga \ No newline at end of file From 5adfcaf354ae9def3cf222e63997fa3df37361fc Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 11:49:41 -0500 Subject: [PATCH 06/13] changed registration name of algorithm to clas12::rga::ProtonEnergyLossCorrection --- .../ProtonEnergyLossCorrection/Config.yaml | 127 ++++++++++++++++++ 1 file changed, 127 insertions(+) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml index e69de29bb..4e4a74328 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml +++ 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 From 0773b1d717e47e9ecc1b634a12f7181c408f3842 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 11:54:15 -0500 Subject: [PATCH 07/13] adjusted file look up location in Algorithm.cc for Config.yaml --- .../clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index 2d878c65b..d51b08467 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -111,7 +111,7 @@ namespace iguana::clas12::rga { // This uses the algorithm's YAMLReader search paths: // - ./, any AddDirectory paths, then install prefix. // It will locate the installed Config.yaml for this algorithm unless overridden. - std::string cfg_path = GetConfig()->FindFile("Config.yaml"); + std::string cfg_path = GetConfig()->FindFile("algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml"); YAML::Node root = YAML::LoadFile(cfg_path); From afe7494f96d60c63999352898ba02090265eb5cb Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 12:41:27 -0500 Subject: [PATCH 08/13] wrote validator for proton energy loss corrections, collects protons in bins of theta and prints before and after average momentum --- .../ProtonEnergyLossCorrection/Validator.cc | 173 ++++++++++++++++++ .../ProtonEnergyLossCorrection/Validator.h | 54 ++++++ src/iguana/algorithms/meson.build | 3 +- 3 files changed, 229 insertions(+), 1 deletion(-) create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc create mode 100644 src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h 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 000000000..5e65fb910 --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc @@ -0,0 +1,173 @@ +#include "Validator.h" + +#include +#include +#include +#include + +namespace iguana::clas12::rga { + + REGISTER_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator); + + double ProtonEnergyLossCorrectionValidator::ThetaDegFromPxPyPz(double px, double py, double pz) + { + // theta = atan2(pt, pz) in radians -> degrees + double pt = std::sqrt(px * px + py * py); + double th = std::atan2(pt, pz); + return th * (180.0 / 3.14159265358979323846); + } + + int ProtonEnergyLossCorrectionValidator::ThetaBinIndex(double theta_deg) + { + if(theta_deg < kThetaMinDeg) { + return -1; + } + if(theta_deg > kThetaMaxDeg) { + return -1; + } + + // Include theta == 70 in the last bin + 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) + { + b_particle = GetBankIndex(banks, "REC::Particle"); + b_config = GetBankIndex(banks, "RUN::config"); + + m_algo_seq = std::make_unique(); + m_algo_seq->Add("clas12::rga::ProtonEnergyLossCorrection"); + m_algo_seq->Start(banks); + + // Reset accumulators (safe here: StartHook is non-const) + 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, b_particle, "REC::Particle"); + auto& config = GetBank(banks, b_config, "RUN::config"); + (void)config; + + // Snapshot BEFORE for protons: store (row, bin, p_before) + 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(); + 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 in place + m_algo_seq->Run(banks); + + // AFTER: same rows; accumulate in shared state + 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() + { + 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"; + } + +} // namespace iguana::clas12::rga \ 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 000000000..a537557ff --- /dev/null +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h @@ -0,0 +1,54 @@ +#pragma once + +#include "iguana/algorithms/Validator.h" +#include "iguana/algorithms/AlgorithmSequence.h" + +#include +#include +#include + +namespace iguana::clas12::rga { + + // Validator for clas12::rga::ProtonEnergyLossCorrection + 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: + // Banks + hipo::banklist::size_type b_particle{}; + hipo::banklist::size_type b_config{}; + + // Run algorithm via AlgorithmSequence + std::unique_ptr m_algo_seq; + + // Theta binning: [5,10), [10,15), ..., [65,70] (deg) + 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); // 13 + + struct BinAccum { + long long n = 0; + double sum_p_before = 0.0; + double sum_p_after = 0.0; + }; + + // NOTE: RunHook is const, but validators accumulate state -> mark accumulators mutable. + mutable std::array m_bins{}; + mutable long long m_total_protons_in_range = 0; + mutable long long m_total_protons_all = 0; + + mutable std::mutex m_mutex{}; + + private: + static int ThetaBinIndex(double theta_deg); + static double ThetaDegFromPxPyPz(double px, double py, double pz); + }; + +} // namespace iguana::clas12::rga \ No newline at end of file diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index bcccfb61a..84efb862e 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -84,7 +84,8 @@ algo_dict = [ }, { 'name': 'clas12::rga::ProtonEnergyLossCorrection', - 'has_validator': false, + 'has_validator': true, + 'validator_needs_ROOT': false, 'test_args': { 'banks': [ 'REC::Particle', 'RUN::config' ], }, From df7202f699ca19b4572ca628a2c0156fbefb3d01 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 13:16:55 -0500 Subject: [PATCH 09/13] wrote validator for proton energy loss corrections, collects protons in bins of theta and prints before and after average momentum --- .../ProtonEnergyLossCorrection/Algorithm.cc | 350 ++++++++++++------ .../ProtonEnergyLossCorrection/Algorithm.h | 132 ++++++- .../ProtonEnergyLossCorrection/Validator.cc | 53 ++- .../ProtonEnergyLossCorrection/Validator.h | 110 ++++-- 4 files changed, 453 insertions(+), 192 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index d51b08467..abab6de0a 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -9,9 +9,35 @@ namespace iguana::clas12::rga { REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); - // -------------------------- - // Small math helpers - // -------------------------- + namespace { + // Avoid relying on non-portable M_PI. + static constexpr double kPi = 3.14159265358979323846; + static constexpr double kRadToDeg = 180.0 / kPi; + static constexpr double kDegToRad = kPi / 180.0; + + // Helper: 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; + } + } // namespace + + // ----------------------------------------------------------------------------- + // Detector-region helpers (based on REC::Particle status) + // ----------------------------------------------------------------------------- + // + // Iguana / CLAS12 conventions: + // - Forward Detector tracks typically have status in [2000, 4000). + // - Central Detector tracks typically have status in [4000, 5000). + // + // We use abs(status) because negative values can appear depending on + // reconstruction conventions. bool ProtonEnergyLossCorrection::IsFD(int status) { int s = std::abs(status); @@ -24,11 +50,16 @@ namespace iguana::clas12::rga { return (s >= 4000 && s < 5000); } + // ----------------------------------------------------------------------------- + // Small 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. + // We compute cos(theta) = pz / |p| and clamp to [-1,1] for numerical safety. double ProtonEnergyLossCorrection::ThetaDeg(double px, double py, double pz) { double r = Pmag(px, py, pz); @@ -42,41 +73,45 @@ namespace iguana::clas12::rga { if(c < -1.0) { c = -1.0; } - return (180.0 / M_PI) * std::acos(c); + return kRadToDeg * std::acos(c); } - // phi = toDegrees(atan2(x, y)); + // phi convention matching the Java implementation: + // + // phi = toDegrees(atan2(px, py)); // phi = phi - 90; // if (phi < 0) phi = 360 + phi; // phi = 360 - phi; - // Returns phi in [0,360). + // + // This yields phi in [0,360). double ProtonEnergyLossCorrection::PhiDegLikeJava(double px, double py) { - double phi = (180.0 / M_PI) * std::atan2(px, py); + double phi = kRadToDeg * std::atan2(px, py); phi = phi - 90.0; if(phi < 0.0) { phi = 360.0 + phi; } phi = 360.0 - phi; - // keep in [0,360) - while(phi >= 360.0) { - phi -= 360.0; - } - while(phi < 0.0) { - phi += 360.0; - } - return phi; + return WrapDeg360(phi); } - void ProtonEnergyLossCorrection::SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& pz) + // Convert from spherical-like (p, theta_deg, phi_deg) back to Cartesian. + 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); + double theta = theta_deg * kDegToRad; + double phi = phi_deg * kDegToRad; + + 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; @@ -88,9 +123,15 @@ namespace iguana::clas12::rga { return sum; } + // ----------------------------------------------------------------------------- + // Period lookup + // ----------------------------------------------------------------------------- + // + // Given a run number, find the first PeriodDef 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 : periods_) { + 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) { @@ -101,116 +142,168 @@ namespace iguana::clas12::rga { return nullptr; } - // -------------------------- - // Hooks - // -------------------------- + // ----------------------------------------------------------------------------- + // 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, we load the FD and CD RegionCoeffs, each of which is + // a set of polynomials in theta for A/B/C of p/theta/phi. void ProtonEnergyLossCorrection::ConfigHook() { - periods_.clear(); + m_periods.clear(); - // This uses the algorithm's YAMLReader search paths: - // - ./, any AddDirectory paths, then install prefix. - // It will locate the installed Config.yaml for this algorithm unless overridden. - std::string cfg_path = GetConfig()->FindFile("algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml"); + // Locate the algorithm's installed Config.yaml, unless overridden by search paths. + std::string const cfg_path = + GetConfig()->FindFile("algorithms/clas12/rga/ProtonEnergyLossCorrection/Config.yaml"); YAML::Node root = YAML::LoadFile(cfg_path); - if(!root["clas12"] || !root["clas12"]["ProtonEnergyLossCorrection"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node clas12:ProtonEnergyLossCorrection"); + // Defensive parsing: provide clear error messages if YAML structure changes. + if(!root["clas12"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing top-level key 'clas12'"); + } + if(!root["clas12"]["ProtonEnergyLossCorrection"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing key 'clas12:ProtonEnergyLossCorrection'"); } YAML::Node node = root["clas12"]["ProtonEnergyLossCorrection"]; if(!node["periods"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: missing YAML node periods"); + throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing key 'periods' under clas12:ProtonEnergyLossCorrection"); } - auto periods_node = node["periods"]; + 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 { + if(!region_node[name]) { + throw std::runtime_error(std::string("ProtonEnergyLossCorrection: missing coefficient list '") + name + "'"); + } + YAML::Node a = region_node[name]; + if(!a.IsSequence()) { + throw std::runtime_error(std::string("ProtonEnergyLossCorrection: coefficient '") + name + "' is not a YAML sequence"); + } + + 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; + // run_ranges required for mapping run -> period if(!pnode["run_ranges"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: missing run_ranges for period " + key); + throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing key 'run_ranges'"); } + for(auto rr : pnode["run_ranges"]) { if(!rr.IsSequence() || rr.size() != 2) { - throw std::runtime_error("ProtonEnergyLossCorrection: run_ranges entry must be [min,max] for period " + key); + throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "': each run_ranges entry must be [min,max]"); } def.run_ranges.push_back({rr[0].as(), rr[1].as()}); } - auto load_poly = [](YAML::Node const& n, char const* name) -> Poly { - Poly p; - if(!n[name]) { - throw std::runtime_error(std::string("ProtonEnergyLossCorrection: missing coeff node ") + name); - } - YAML::Node a = n[name]; - if(!a.IsSequence()) { - throw std::runtime_error(std::string("ProtonEnergyLossCorrection: coeff node not sequence: ") + name); - } - for(auto v : a) { - p.c.push_back(v.as()); - } - return p; - }; - - 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; - }; - - if(!pnode["FD"] || !pnode["CD"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: missing FD or CD node for period " + key); + // Both regions must exist. + if(!pnode["FD"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing 'FD' block"); + } + if(!pnode["CD"]) { + throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing 'CD' block"); } def.fd = load_region(pnode["FD"]); def.cd = load_region(pnode["CD"]); - periods_[key] = def; + m_periods[key] = def; } } + // ----------------------------------------------------------------------------- + // StartHook: cache bank indices + // ----------------------------------------------------------------------------- void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) { - b_rec_particle = GetBankIndex(banks, "REC::Particle"); - b_run_config = GetBankIndex(banks, "RUN::config"); + 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, b_rec_particle, "REC::Particle"); - auto& run = GetBank(banks, b_run_config, "RUN::config"); + 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)); @@ -221,98 +314,111 @@ namespace iguana::clas12::rga { void ProtonEnergyLossCorrection::StopHook() { + // No summary output here by default. Validation and monitoring are done + // by the separate Validator class. } - // -------------------------- - // Transform - // -------------------------- - 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 + // ----------------------------------------------------------------------------- + // 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 + // Defensive: 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}; } - double px = px_in; - double py = py_in; - double pz = 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); - double phi = PhiDegLikeJava(px, py); - - RegionCoeffs const& c = is_fd ? period->fd : period->cd; - - // Build A/B/C from theta polynomials - double A_p = EvalPoly(c.A_p, theta); - double B_p = EvalPoly(c.B_p, theta); - double C_p = EvalPoly(c.C_p, theta); - - double A_theta = EvalPoly(c.A_theta, theta); - double B_theta = EvalPoly(c.B_theta, theta); - double C_theta = EvalPoly(c.C_theta, theta); - - double A_phi = EvalPoly(c.A_phi, theta); - double B_phi = EvalPoly(c.B_phi, theta); - double C_phi = EvalPoly(c.C_phi, theta); - - // Apply corrections using your exact functional forms + double theta = ThetaDeg(px, py, pz); // degrees + double phi = PhiDegLikeJava(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 + // theta_new and phi_new use the same inverse-angle form as FD. 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); + p_new += A_p + (B_p / p) + (C_p / (p * p)); - // Guard theta,phi divisions + // Avoid divide-by-zero (theta and phi are in degrees). if(theta != 0.0) { - theta_new += A_theta + B_theta / theta + C_theta / (theta * theta); + 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 += A_phi + (B_phi / phi) + (C_phi / (phi * phi)); } - } else if(is_cd) { - p_new += A_p + B_p * p + C_p * p * p; + } 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); + 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 += A_phi + (B_phi / phi) + (C_phi / (phi * phi)); } } - // phi convention: keep in [0,360) since that is what your Java uses - while(phi_new >= 360.0) { - phi_new -= 360.0; - } - while(phi_new < 0.0) { - phi_new += 360.0; - } + // Keep phi consistent with the Java convention. + phi_new = WrapDeg360(phi_new); - // Convert back to Cartesian using spherical conversion - double px_out = 0.0, py_out = 0.0, pz_out = 0.0; + // 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}; diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h index f82106e5d..a8d641217 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -1,5 +1,50 @@ #pragma once +// ProtonEnergyLossCorrection +// +// Purpose +// ------- +// Apply momentum and angular corrections to reconstructed protons in CLAS12 RGA +// data, reproducing the same functional forms used in a corresponding Java +// implementation. +// +// Where it runs +// ------------- +// This Iguana algorithm 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, +// we 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 +// +// Coordinate conventions +// ---------------------- +// theta is computed from the momentum direction and returned in degrees. +// phi is computed using a convention that matches the Java code: +// +// phi = toDegrees(atan2(px, py)); +// phi = phi - 90; +// if (phi < 0) phi = 360 + phi; +// phi = 360 - phi; +// +// This yields phi in [0,360) degrees. +// +// Notes on correctness / expectations +// ---------------------------------- +// - This algorithm does not try to "fix" events outside known periods; +// if the run is not in any configured run range, it leaves the row unchanged. +// - This algorithm only corrects protons in FD or CD as determined by status. +// - It does not depend on ROOT. + #include "iguana/algorithms/Algorithm.h" #include "iguana/algorithms/TypeDefs.h" @@ -11,13 +56,30 @@ namespace iguana::clas12::rga { - class ProtonEnergyLossCorrection : public Algorithm { + class ProtonEnergyLossCorrection : public Algorithm + { DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection) public: // @action_function{scalar transformer} - // Transform a single particle row (px,py,pz) if it is an RGA proton. - // Returns corrected (px,py,pz). + // + // Transform + // --------- + // Apply the configured correction to a single particle described by: + // pid, status, run, (px,py,pz). + // + // This is written as a scalar action so it can be used by Iguana's action + // system, but it is also used internally by RunHook(). + // + // Inputs: + // pid : PDG ID (we correct only protons: 2212) + // status : REC::Particle status code (used to identify FD vs CD) + // run : run number (selects the period / coefficients) + // px,py,pz : momentum components (GeV) + // + // Output: + // corrected (px,py,pz). If the row should not be corrected, the input + // values are returned unchanged. std::tuple Transform( int const pid, int const status, @@ -27,46 +89,82 @@ namespace iguana::clas12::rga { vector_element_t const pz) const; private: - // Hook API (Start/Run/Stop are final in base class) + // Iguana hook API + // --------------- + // StartHook / RunHook / StopHook are called during processing. + // ConfigHook is called when configuration is loaded/updated. void ConfigHook() override; void StartHook(hipo::banklist& banks) override; bool RunHook(hipo::banklist& banks) const override; void StopHook() override; - // Bank indices - hipo::banklist::size_type b_rec_particle; - hipo::banklist::size_type b_run_config; + // Cached bank indices (resolved once in StartHook) + hipo::banklist::size_type m_b_rec_particle{}; + hipo::banklist::size_type m_b_run_config{}; + + // Internal configuration representation + // ------------------------------------- + // We read YAML coefficients into these structs for fast lookup. - // Internal config structures + // Simple polynomial container: + // c0 + c1*x + c2*x^2 + ... struct Poly { - std::vector c; // c0 + c1*x + c2*x^2 + ... + std::vector c; }; + // Coefficients for a single detector region (FD or CD). + // + // For each of p, theta, phi we 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, B_p, C_p; - Poly A_theta, B_theta, C_theta; - Poly A_phi, B_phi, C_phi; + 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; }; - // period key -> period def - std::map periods_; + // Map period key -> period definition (loaded from YAML) + std::map m_periods{}; - // Helpers + // Helper functions (pure utilities) 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); + + // phi convention that matches the Java implementation (see header comment). static double PhiDegLikeJava(double px, double py); - static void SphericalToCartesian(double p, double theta_deg, double phi_deg, double& px, double& py, double& pz); - // Find matching period for run, returns nullptr if none + // 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 PeriodDef or nullptr if none. PeriodDef const* FindPeriod(int run) const; }; diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc index 5e65fb910..bb99229ac 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc @@ -7,16 +7,31 @@ namespace iguana::clas12::rga { + namespace { + static constexpr double kPi = 3.14159265358979323846; + static constexpr double kRadToDeg = 180.0 / kPi; + } // namespace + REGISTER_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator); + // Compute theta (deg) from momentum components. + // + // theta = atan2(pt, pz) + // where pt = sqrt(px^2 + py^2) + // + // This is robust for pz near 0 and matches common CLAS usage. double ProtonEnergyLossCorrectionValidator::ThetaDegFromPxPyPz(double px, double py, double pz) { - // theta = atan2(pt, pz) in radians -> degrees double pt = std::sqrt(px * px + py * py); double th = std::atan2(pt, pz); - return th * (180.0 / 3.14159265358979323846); + return th * kRadToDeg; } + // 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) { @@ -26,7 +41,7 @@ namespace iguana::clas12::rga { return -1; } - // Include theta == 70 in the last bin + // Include theta == kThetaMaxDeg in the final bin. if(theta_deg >= kThetaMaxDeg) { return kNBins - 1; } @@ -40,14 +55,16 @@ namespace iguana::clas12::rga { void ProtonEnergyLossCorrectionValidator::StartHook(hipo::banklist& banks) { - b_particle = GetBankIndex(banks, "REC::Particle"); - b_config = GetBankIndex(banks, "RUN::config"); + // Cache bank indices for fast access during RunHook(). + m_b_particle = GetBankIndex(banks, "REC::Particle"); + m_b_config = GetBankIndex(banks, "RUN::config"); + // Build and start an AlgorithmSequence containing only the algorithm under test. m_algo_seq = std::make_unique(); m_algo_seq->Add("clas12::rga::ProtonEnergyLossCorrection"); m_algo_seq->Start(banks); - // Reset accumulators (safe here: StartHook is non-const) + // Reset counters/accumulators. for(auto& b : m_bins) { b.n = 0; b.sum_p_before = 0.0; @@ -59,21 +76,25 @@ namespace iguana::clas12::rga { bool ProtonEnergyLossCorrectionValidator::RunHook(hipo::banklist& banks) const { - auto& particle = GetBank(banks, b_particle, "REC::Particle"); - auto& config = GetBank(banks, b_config, "RUN::config"); - (void)config; + 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 this validator summary. - // Snapshot BEFORE for protons: store (row, bin, p_before) + // We must record BEFORE values, then run the algorithm, then compute AFTER. + // Store the particle row index and p_before so we can match exactly the same + // rows after the algorithm modifies the bank in place. struct Entry { - int row = -1; - int bin = -1; - double p_before = 0.0; + 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) { @@ -95,10 +116,11 @@ namespace iguana::clas12::rga { entries.push_back(e); } - // Run algorithm in place + // Run algorithm under test (in place on the same banks). m_algo_seq->Run(banks); - // AFTER: same rows; accumulate in shared state + // Accumulate AFTER values. Use a mutex since multiple events may be processed + // concurrently depending on Iguana configuration. std::scoped_lock lock(m_mutex); for(auto const& e : entries) { @@ -128,6 +150,7 @@ namespace iguana::clas12::rga { void ProtonEnergyLossCorrectionValidator::StopHook() { + // Print a compact, human-readable summary suitable for CI logs. std::cout << "\n"; std::cout << "ProtonEnergyLossCorrectionValidator summary\n"; std::cout << " theta bins: " << kThetaMinDeg << " to " << kThetaMaxDeg diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h index a537557ff..1ec4f065e 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h @@ -1,5 +1,33 @@ #pragma once +// ProtonEnergyLossCorrectionValidator +// +// Purpose +// ------- +// Provide a simple, ROOT-free validation summary for the proton energy loss +// correction algorithm. +// +// What it does +// ------------ +// 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, , , +// +// Threading model +// --------------- +// Iguana may run validators in multi-threaded contexts, so accumulation is +// protected by a mutex. The accumulators are marked "mutable" because Iguana's +// validator API uses RunHook(...) const, but validators naturally accumulate +// per-run state. + #include "iguana/algorithms/Validator.h" #include "iguana/algorithms/AlgorithmSequence.h" @@ -9,46 +37,52 @@ namespace iguana::clas12::rga { - // Validator for clas12::rga::ProtonEnergyLossCorrection 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: - // Banks - hipo::banklist::size_type b_particle{}; - hipo::banklist::size_type b_config{}; - - // Run algorithm via AlgorithmSequence - std::unique_ptr m_algo_seq; - - // Theta binning: [5,10), [10,15), ..., [65,70] (deg) - 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); // 13 - - struct BinAccum { - long long n = 0; - double sum_p_before = 0.0; - double sum_p_after = 0.0; - }; - - // NOTE: RunHook is const, but validators accumulate state -> mark accumulators mutable. - mutable std::array m_bins{}; - mutable long long m_total_protons_in_range = 0; - mutable long long m_total_protons_all = 0; - - mutable std::mutex m_mutex{}; - - private: - static int ThetaBinIndex(double theta_deg); - static double ThetaDegFromPxPyPz(double px, double py, double pz); + 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{}; + + // We run the algorithm under test via an AlgorithmSequence so the validator + // does not need to know algorithm internals. + std::unique_ptr m_algo_seq; + + // Theta binning configuration (degrees). + // Bins: [5,10), [10,15), ..., [65,70] with 70 included in last bin. + 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); // 13 + + // Simple per-bin accumulators. + struct BinAccum { + long long n = 0; + double sum_p_before = 0.0; + double sum_p_after = 0.0; + }; + + // RunHook is const, but we 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); }; } // namespace iguana::clas12::rga \ No newline at end of file From 9f577fd9c90c3a28200a95dc84aa6ddefd3ea0f1 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 13:26:11 -0500 Subject: [PATCH 10/13] numerous documentation/comment updates; renamed phideglikejava to phideg (still rotates) --- .../ProtonEnergyLossCorrection/Algorithm.cc | 4 +- .../ProtonEnergyLossCorrection/Algorithm.h | 48 ++----------------- 2 files changed, 7 insertions(+), 45 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index abab6de0a..f1b040d5a 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -84,7 +84,7 @@ namespace iguana::clas12::rga { // phi = 360 - phi; // // This yields phi in [0,360). - double ProtonEnergyLossCorrection::PhiDegLikeJava(double px, double py) + double ProtonEnergyLossCorrection::PhiDeg(double px, double py) { double phi = kRadToDeg * std::atan2(px, py); phi = phi - 90.0; @@ -358,7 +358,7 @@ namespace iguana::clas12::rga { } double theta = ThetaDeg(px, py, pz); // degrees - double phi = PhiDegLikeJava(px, py); // degrees in [0,360) + double phi = PhiDeg(px, py); // degrees in [0,360) // Choose FD vs CD coefficients. RegionCoeffs const& coeffs = is_fd ? period->fd : period->cd; diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h index a8d641217..330750630 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -2,19 +2,11 @@ // ProtonEnergyLossCorrection // -// Purpose -// ------- // Apply momentum and angular corrections to reconstructed protons in CLAS12 RGA -// data, reproducing the same functional forms used in a corresponding Java -// implementation. -// -// Where it runs -// ------------- -// This Iguana algorithm runs event-by-event and modifies REC::Particle in place. +// 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, -// we compute p, theta, phi from (px,py,pz), -// apply period-dependent corrections, +// compute p, theta, phi from (px,py,pz), apply period-dependent corrections, // then write back corrected (px,py,pz). // // Period dependence @@ -26,24 +18,6 @@ // - FD coefficients // - CD coefficients // -// Coordinate conventions -// ---------------------- -// theta is computed from the momentum direction and returned in degrees. -// phi is computed using a convention that matches the Java code: -// -// phi = toDegrees(atan2(px, py)); -// phi = phi - 90; -// if (phi < 0) phi = 360 + phi; -// phi = 360 - phi; -// -// This yields phi in [0,360) degrees. -// -// Notes on correctness / expectations -// ---------------------------------- -// - This algorithm does not try to "fix" events outside known periods; -// if the run is not in any configured run range, it leaves the row unchanged. -// - This algorithm only corrects protons in FD or CD as determined by status. -// - It does not depend on ROOT. #include "iguana/algorithms/Algorithm.h" #include "iguana/algorithms/TypeDefs.h" @@ -63,14 +37,6 @@ namespace iguana::clas12::rga { public: // @action_function{scalar transformer} // - // Transform - // --------- - // Apply the configured correction to a single particle described by: - // pid, status, run, (px,py,pz). - // - // This is written as a scalar action so it can be used by Iguana's action - // system, but it is also used internally by RunHook(). - // // Inputs: // pid : PDG ID (we correct only protons: 2212) // status : REC::Particle status code (used to identify FD vs CD) @@ -78,8 +44,7 @@ namespace iguana::clas12::rga { // px,py,pz : momentum components (GeV) // // Output: - // corrected (px,py,pz). If the row should not be corrected, the input - // values are returned unchanged. + // corrected (px,py,pz). std::tuple Transform( int const pid, int const status, @@ -91,8 +56,6 @@ namespace iguana::clas12::rga { private: // Iguana hook API // --------------- - // StartHook / RunHook / StopHook are called during processing. - // ConfigHook is called when configuration is loaded/updated. void ConfigHook() override; void StartHook(hipo::banklist& banks) override; bool RunHook(hipo::banklist& banks) const override; @@ -106,8 +69,7 @@ namespace iguana::clas12::rga { // ------------------------------------- // We read YAML coefficients into these structs for fast lookup. - // Simple polynomial container: - // c0 + c1*x + c2*x^2 + ... + // Simple polynomial container: c0 + c1*x + c2*x^2 + ... struct Poly { std::vector c; }; @@ -154,7 +116,7 @@ namespace iguana::clas12::rga { static double ThetaDeg(double px, double py, double pz); // phi convention that matches the Java implementation (see header comment). - static double PhiDegLikeJava(double px, double py); + static double PhiDeg(double px, double py); // Convert (p, theta_deg, phi_deg) back to (px,py,pz). static void SphericalToCartesian(double p, From 44e678b36d8e484ab3cd7b39921c93ce2e8ad4ef Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 13:47:18 -0500 Subject: [PATCH 11/13] further documentation updates for user clarity --- .../ProtonEnergyLossCorrection/Algorithm.cc | 135 ++++++------------ .../ProtonEnergyLossCorrection/Algorithm.h | 14 +- 2 files changed, 45 insertions(+), 104 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index f1b040d5a..2c4b9d3c0 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -10,12 +10,10 @@ namespace iguana::clas12::rga { REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); namespace { - // Avoid relying on non-portable M_PI. - static constexpr double kPi = 3.14159265358979323846; - static constexpr double kRadToDeg = 180.0 / kPi; - static constexpr double kDegToRad = kPi / 180.0; - // Helper: keep an angle in degrees within [0,360). + // Keep an angle in degrees within [0, 360). + // + // We keep this helper because: double WrapDeg360(double x) { while(x >= 360.0) { @@ -26,18 +24,17 @@ namespace iguana::clas12::rga { } return x; } + } // namespace // ----------------------------------------------------------------------------- // Detector-region helpers (based on REC::Particle status) // ----------------------------------------------------------------------------- // - // Iguana / CLAS12 conventions: - // - Forward Detector tracks typically have status in [2000, 4000). - // - Central Detector tracks typically have status in [4000, 5000). + // CLAS12 status definitions: + // - FD tracks: abs(status) in [2000, 4000) + // - CD tracks: abs(status) in [4000, 5000) // - // We use abs(status) because negative values can appear depending on - // reconstruction conventions. bool ProtonEnergyLossCorrection::IsFD(int status) { int s = std::abs(status); @@ -51,7 +48,7 @@ namespace iguana::clas12::rga { } // ----------------------------------------------------------------------------- - // Small math helpers (p, theta, phi) + // math helpers (p, theta, phi) // ----------------------------------------------------------------------------- double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) { @@ -59,13 +56,14 @@ namespace iguana::clas12::rga { } // theta from momentum vector, in degrees. - // We compute cos(theta) = pz / |p| and clamp to [-1,1] for numerical safety. + // We compute: 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; @@ -73,38 +71,33 @@ namespace iguana::clas12::rga { if(c < -1.0) { c = -1.0; } - return kRadToDeg * std::acos(c); + + return (180.0 / M_PI) * std::acos(c); } - // phi convention matching the Java implementation: - // // phi = toDegrees(atan2(px, py)); // phi = phi - 90; // if (phi < 0) phi = 360 + phi; // phi = 360 - phi; // - // This yields phi in [0,360). + // Returns phi in [0,360). double ProtonEnergyLossCorrection::PhiDeg(double px, double py) { - double phi = kRadToDeg * std::atan2(px, 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 spherical-like (p, theta_deg, phi_deg) back to Cartesian. - void ProtonEnergyLossCorrection::SphericalToCartesian(double p, - double theta_deg, - double phi_deg, - double& px, - double& py, - double& pz) - { - double theta = theta_deg * kDegToRad; - double phi = phi_deg * kDegToRad; + // 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); @@ -112,8 +105,7 @@ namespace iguana::clas12::rga { } // Evaluate polynomial c0 + c1*x + c2*x^2 + ... - double ProtonEnergyLossCorrection::EvalPoly(Poly const& p, double x) - { + double ProtonEnergyLossCorrection::EvalPoly(Poly const& p, double x) { double sum = 0.0; double xn = 1.0; for(double c : p.c) { @@ -127,10 +119,10 @@ namespace iguana::clas12::rga { // Period lookup // ----------------------------------------------------------------------------- // - // Given a run number, find the first PeriodDef whose run_ranges contain it. + // 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 - { + 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) { @@ -163,43 +155,22 @@ namespace iguana::clas12::rga { // A_p: [ ... ] // ... // - // For each period, we load the FD and CD RegionCoeffs, each of which is - // a set of polynomials in theta for A/B/C of p/theta/phi. - void ProtonEnergyLossCorrection::ConfigHook() - { + // For each period, load FD and CD RegionCoeffs. Each RegionCoeffs holds + // polynomial parameterizations in theta for A/B/C of p/theta/phi corrections. + void ProtonEnergyLossCorrection::ConfigHook() { m_periods.clear(); - // Locate the algorithm's installed Config.yaml, unless overridden by search paths. + // 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); - - // Defensive parsing: provide clear error messages if YAML structure changes. - if(!root["clas12"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing top-level key 'clas12'"); - } - if(!root["clas12"]["ProtonEnergyLossCorrection"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing key 'clas12:ProtonEnergyLossCorrection'"); - } - YAML::Node node = root["clas12"]["ProtonEnergyLossCorrection"]; - if(!node["periods"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: YAML missing key 'periods' under 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 { - if(!region_node[name]) { - throw std::runtime_error(std::string("ProtonEnergyLossCorrection: missing coefficient list '") + name + "'"); - } YAML::Node a = region_node[name]; - if(!a.IsSequence()) { - throw std::runtime_error(std::string("ProtonEnergyLossCorrection: coefficient '") + name + "' is not a YAML sequence"); - } - Poly p; p.c.reserve((size_t)a.size()); for(auto v : a) { @@ -234,26 +205,10 @@ namespace iguana::clas12::rga { PeriodDef def; - // run_ranges required for mapping run -> period - if(!pnode["run_ranges"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing key 'run_ranges'"); - } - for(auto rr : pnode["run_ranges"]) { - if(!rr.IsSequence() || rr.size() != 2) { - throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "': each run_ranges entry must be [min,max]"); - } def.run_ranges.push_back({rr[0].as(), rr[1].as()}); } - // Both regions must exist. - if(!pnode["FD"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing 'FD' block"); - } - if(!pnode["CD"]) { - throw std::runtime_error("ProtonEnergyLossCorrection: period '" + key + "' missing 'CD' block"); - } - def.fd = load_region(pnode["FD"]); def.cd = load_region(pnode["CD"]); @@ -264,8 +219,7 @@ namespace iguana::clas12::rga { // ----------------------------------------------------------------------------- // StartHook: cache bank indices // ----------------------------------------------------------------------------- - void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) - { + void ProtonEnergyLossCorrection::StartHook(hipo::banklist& banks) { m_b_rec_particle = GetBankIndex(banks, "REC::Particle"); m_b_run_config = GetBankIndex(banks, "RUN::config"); } @@ -273,8 +227,7 @@ namespace iguana::clas12::rga { // ----------------------------------------------------------------------------- // RunHook: apply correction to each matching REC::Particle row // ----------------------------------------------------------------------------- - bool ProtonEnergyLossCorrection::RunHook(hipo::banklist& banks) const - { + 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"); @@ -312,24 +265,17 @@ namespace iguana::clas12::rga { return true; } - void ProtonEnergyLossCorrection::StopHook() - { - // No summary output here by default. Validation and monitoring are done - // by the separate Validator class. + 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 - { - // Defensive: only protons are corrected. + 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}; } @@ -357,8 +303,8 @@ namespace iguana::clas12::rga { return {px_in, py_in, pz_in}; } - double theta = ThetaDeg(px, py, pz); // degrees - double phi = PhiDeg(px, py); // degrees in [0,360) + double theta = ThetaDeg(px, py, pz); // degrees + double phi = PhiDegLikeJava(px, py); // degrees in [0,360) // Choose FD vs CD coefficients. RegionCoeffs const& coeffs = is_fd ? period->fd : period->cd; @@ -385,7 +331,7 @@ namespace iguana::clas12::rga { // // Central Detector (CD): // p_new = p + A_p + B_p*p + C_p*p^2 - // theta_new and phi_new use the same inverse-angle form as FD. + // not all A_, B_ and C_ need be nonzero double p_new = p; double theta_new = theta; double phi_new = phi; @@ -393,7 +339,7 @@ namespace iguana::clas12::rga { if(is_fd) { p_new += A_p + (B_p / p) + (C_p / (p * p)); - // Avoid divide-by-zero (theta and phi are in degrees). + // Avoid divide-by-zero if(theta != 0.0) { theta_new += A_theta + (B_theta / theta) + (C_theta / (theta * theta)); } @@ -412,7 +358,6 @@ namespace iguana::clas12::rga { } } - // Keep phi consistent with the Java convention. phi_new = WrapDeg360(phi_new); // Convert back to Cartesian using the corrected spherical variables. @@ -424,4 +369,4 @@ namespace iguana::clas12::rga { return {px_out, py_out, pz_out}; } -} // namespace iguana::clas12::rga \ No newline at end of file +} \ 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 index 330750630..0b4e65e31 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -119,15 +119,11 @@ namespace iguana::clas12::rga { 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 PeriodDef or nullptr if none. + 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; }; -} // namespace iguana::clas12::rga \ No newline at end of file +} \ No newline at end of file From 6b4907233bbeb03e2124799432654452651f1eb5 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 13:49:39 -0500 Subject: [PATCH 12/13] further documentation updates for user clarity --- .../clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index 2c4b9d3c0..4ffdb7a03 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -304,7 +304,7 @@ namespace iguana::clas12::rga { } double theta = ThetaDeg(px, py, pz); // degrees - double phi = PhiDegLikeJava(px, py); // degrees in [0,360) + double phi = PhiDeg(px, py); // degrees in [0,360) // Choose FD vs CD coefficients. RegionCoeffs const& coeffs = is_fd ? period->fd : period->cd; From c6d3490676d68d1748d671beea10a3b68d68d6d4 Mon Sep 17 00:00:00 2001 From: "Timothy B. Hayward" Date: Wed, 11 Feb 2026 16:04:54 -0500 Subject: [PATCH 13/13] potential final clean up of comments and documentation --- .../ProtonEnergyLossCorrection/Algorithm.cc | 49 ++++++------------ .../ProtonEnergyLossCorrection/Algorithm.h | 26 +++------- .../ProtonEnergyLossCorrection/Validator.cc | 50 ++++++------------- .../ProtonEnergyLossCorrection/Validator.h | 26 ++-------- 4 files changed, 43 insertions(+), 108 deletions(-) diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc index 4ffdb7a03..82d217346 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.cc @@ -10,22 +10,13 @@ namespace iguana::clas12::rga { REGISTER_IGUANA_ALGORITHM(ProtonEnergyLossCorrection); namespace { - // Keep an angle in degrees within [0, 360). - // - // We keep this helper because: - double WrapDeg360(double x) - { - while(x >= 360.0) { - x -= 360.0; - } - while(x < 0.0) { - x += 360.0; - } + double WrapDeg360(double x) { + while(x >= 360.0) { x -= 360.0; } + while(x < 0.0) { x += 360.0; } return x; } - - } // namespace + } // ----------------------------------------------------------------------------- // Detector-region helpers (based on REC::Particle status) @@ -35,14 +26,12 @@ namespace iguana::clas12::rga { // - FD tracks: abs(status) in [2000, 4000) // - CD tracks: abs(status) in [4000, 5000) // - bool ProtonEnergyLossCorrection::IsFD(int status) - { + bool ProtonEnergyLossCorrection::IsFD(int status) { int s = std::abs(status); return (s >= 2000 && s < 4000); } - bool ProtonEnergyLossCorrection::IsCD(int status) - { + bool ProtonEnergyLossCorrection::IsCD(int status) { int s = std::abs(status); return (s >= 4000 && s < 5000); } @@ -50,15 +39,13 @@ namespace iguana::clas12::rga { // ----------------------------------------------------------------------------- // math helpers (p, theta, phi) // ----------------------------------------------------------------------------- - double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) - { + double ProtonEnergyLossCorrection::Pmag(double px, double py, double pz) { return std::sqrt(px * px + py * py + pz * pz); } // theta from momentum vector, in degrees. - // We compute: cos(theta) = pz / |p| - double ProtonEnergyLossCorrection::ThetaDeg(double px, double py, double pz) - { + // 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; @@ -75,16 +62,11 @@ namespace iguana::clas12::rga { return (180.0 / M_PI) * std::acos(c); } - // phi = toDegrees(atan2(px, py)); - // phi = phi - 90; - // if (phi < 0) phi = 360 + phi; - // phi = 360 - phi; - // + // Returns phi in [0,360). - double ProtonEnergyLossCorrection::PhiDeg(double px, double py) - { + double ProtonEnergyLossCorrection::PhiDeg(double px, double py) { double phi = (180.0 / M_PI) * std::atan2(px, py); - phi = phi - 90.0; + phi = phi - 90.0; if(phi < 0.0) { phi = 360.0 + phi; } @@ -157,6 +139,7 @@ namespace iguana::clas12::rga { // // 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(); @@ -302,9 +285,8 @@ namespace iguana::clas12::rga { 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) + 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; @@ -339,7 +321,6 @@ namespace iguana::clas12::rga { if(is_fd) { p_new += A_p + (B_p / p) + (C_p / (p * p)); - // Avoid divide-by-zero if(theta != 0.0) { theta_new += A_theta + (B_theta / theta) + (C_theta / (theta * theta)); } diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h index 0b4e65e31..a9a333f10 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Algorithm.h @@ -5,7 +5,7 @@ // 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, +// - 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). // @@ -30,17 +30,16 @@ namespace iguana::clas12::rga { - class ProtonEnergyLossCorrection : public Algorithm - { + class ProtonEnergyLossCorrection : public Algorithm { DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection) public: // @action_function{scalar transformer} // // Inputs: - // pid : PDG ID (we correct only protons: 2212) - // status : REC::Particle status code (used to identify FD vs CD) - // run : run number (selects the period / coefficients) + // 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: @@ -54,8 +53,6 @@ namespace iguana::clas12::rga { vector_element_t const pz) const; private: - // Iguana hook API - // --------------- void ConfigHook() override; void StartHook(hipo::banklist& banks) override; bool RunHook(hipo::banklist& banks) const override; @@ -65,20 +62,13 @@ namespace iguana::clas12::rga { hipo::banklist::size_type m_b_rec_particle{}; hipo::banklist::size_type m_b_run_config{}; - // Internal configuration representation - // ------------------------------------- - // We read YAML coefficients into these structs for fast lookup. - // 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 we store (A,B,C) polynomials in theta: - // A(theta), B(theta), C(theta) - // + // 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; @@ -106,7 +96,7 @@ namespace iguana::clas12::rga { // Map period key -> period definition (loaded from YAML) std::map m_periods{}; - // Helper functions (pure utilities) + // Helper functions static bool IsFD(int status); static bool IsCD(int status); @@ -114,8 +104,6 @@ namespace iguana::clas12::rga { static double Pmag(double px, double py, double pz); static double ThetaDeg(double px, double py, double pz); - - // phi convention that matches the Java implementation (see header comment). static double PhiDeg(double px, double py); // Convert (p, theta_deg, phi_deg) back to (px,py,pz). diff --git a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc index bb99229ac..a0422e48c 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.cc @@ -7,24 +7,13 @@ namespace iguana::clas12::rga { - namespace { - static constexpr double kPi = 3.14159265358979323846; - static constexpr double kRadToDeg = 180.0 / kPi; - } // namespace - REGISTER_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator); // Compute theta (deg) from momentum components. - // - // theta = atan2(pt, pz) - // where pt = sqrt(px^2 + py^2) - // - // This is robust for pz near 0 and matches common CLAS usage. - double ProtonEnergyLossCorrectionValidator::ThetaDegFromPxPyPz(double px, double py, double pz) - { + 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 * kRadToDeg; + return th * (180.0 / M_PI); } // Convert theta (deg) to a bin index. @@ -32,16 +21,13 @@ namespace iguana::clas12::rga { // Returns: // -1 if theta is outside [kThetaMinDeg, kThetaMaxDeg] // otherwise an index in [0, kNBins-1]. - int ProtonEnergyLossCorrectionValidator::ThetaBinIndex(double theta_deg) - { + int ProtonEnergyLossCorrectionValidator::ThetaBinIndex(double theta_deg) { if(theta_deg < kThetaMinDeg) { return -1; } if(theta_deg > kThetaMaxDeg) { return -1; } - - // Include theta == kThetaMaxDeg in the final bin. if(theta_deg >= kThetaMaxDeg) { return kNBins - 1; } @@ -53,13 +39,12 @@ namespace iguana::clas12::rga { return idx; } - void ProtonEnergyLossCorrectionValidator::StartHook(hipo::banklist& banks) - { - // Cache bank indices for fast access during RunHook(). + 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 containing only the algorithm under test. + // Build and start an AlgorithmSequence m_algo_seq = std::make_unique(); m_algo_seq->Add("clas12::rga::ProtonEnergyLossCorrection"); m_algo_seq->Start(banks); @@ -74,14 +59,13 @@ namespace iguana::clas12::rga { m_total_protons_all = 0; } - bool ProtonEnergyLossCorrectionValidator::RunHook(hipo::banklist& banks) const - { + 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 this validator summary. + (void)config; // run number is not needed for summary of \Delta p. - // We must record BEFORE values, then run the algorithm, then compute AFTER. - // Store the particle row index and p_before so we can match exactly the same + // 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; @@ -94,7 +78,7 @@ namespace iguana::clas12::rga { int const nrows = particle.getRows(); - // Snapshot BEFORE for all protons. + // snapshot before for all protons. for(int i = 0; i < nrows; ++i) { int pid = particle.getInt("pid", i); if(pid != 2212) { @@ -116,11 +100,10 @@ namespace iguana::clas12::rga { entries.push_back(e); } - // Run algorithm under test (in place on the same banks). + // Run algorithm under test m_algo_seq->Run(banks); - // Accumulate AFTER values. Use a mutex since multiple events may be processed - // concurrently depending on Iguana configuration. + // Accumulate after values. std::scoped_lock lock(m_mutex); for(auto const& e : entries) { @@ -148,9 +131,8 @@ namespace iguana::clas12::rga { return true; } - void ProtonEnergyLossCorrectionValidator::StopHook() - { - // Print a compact, human-readable summary suitable for CI logs. + void ProtonEnergyLossCorrectionValidator::StopHook() { + // human-readable summary table std::cout << "\n"; std::cout << "ProtonEnergyLossCorrectionValidator summary\n"; std::cout << " theta bins: " << kThetaMinDeg << " to " << kThetaMaxDeg @@ -193,4 +175,4 @@ namespace iguana::clas12::rga { std::cout << "\n"; } -} // namespace iguana::clas12::rga \ No newline at end of file +} \ 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 index 1ec4f065e..a7cd1ed65 100644 --- a/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h +++ b/src/iguana/algorithms/clas12/rga/ProtonEnergyLossCorrection/Validator.h @@ -2,13 +2,6 @@ // ProtonEnergyLossCorrectionValidator // -// Purpose -// ------- -// Provide a simple, ROOT-free validation summary for the proton energy loss -// correction algorithm. -// -// What it does -// ------------ // For each processed event: // 1) Find all protons (pid==2212) in REC::Particle. // 2) Compute their theta before correction, place them into theta bins. @@ -21,12 +14,6 @@ // In StopHook(), print a table: // theta bin range, N, , , // -// Threading model -// --------------- -// Iguana may run validators in multi-threaded contexts, so accumulation is -// protected by a mutex. The accumulators are marked "mutable" because Iguana's -// validator API uses RunHook(...) const, but validators naturally accumulate -// per-run state. #include "iguana/algorithms/Validator.h" #include "iguana/algorithms/AlgorithmSequence.h" @@ -37,8 +24,7 @@ namespace iguana::clas12::rga { - class ProtonEnergyLossCorrectionValidator : public Validator - { + class ProtonEnergyLossCorrectionValidator : public Validator { DEFINE_IGUANA_VALIDATOR(ProtonEnergyLossCorrectionValidator, clas12::rga::ProtonEnergyLossCorrectionValidator) private: @@ -51,16 +37,14 @@ namespace iguana::clas12::rga { hipo::banklist::size_type m_b_particle{}; hipo::banklist::size_type m_b_config{}; - // We run the algorithm under test via an AlgorithmSequence so the validator - // does not need to know algorithm internals. std::unique_ptr m_algo_seq; // Theta binning configuration (degrees). - // Bins: [5,10), [10,15), ..., [65,70] with 70 included in last bin. + // 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); // 13 + static constexpr int kNBins = (int)((kThetaMaxDeg - kThetaMinDeg) / kThetaStepDeg); // Simple per-bin accumulators. struct BinAccum { @@ -69,7 +53,7 @@ namespace iguana::clas12::rga { double sum_p_after = 0.0; }; - // RunHook is const, but we accumulate over events -> these must be mutable. + // 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; @@ -85,4 +69,4 @@ namespace iguana::clas12::rga { static double ThetaDegFromPxPyPz(double px, double py, double pz); }; -} // namespace iguana::clas12::rga \ No newline at end of file +} \ No newline at end of file