Algorithm.h

#pragma once
 
// ProtonEnergyLossCorrection
//
// Apply momentum and angular corrections to reconstructed protons in CLAS12 RGA
// data; runs event-by-event and modifies REC::Particle in place.
// For each row in REC::Particle:
//   - If pid == 2212 (proton) AND status indicates FD or CD (should always),
//     compute p, theta, phi from (px,py,pz), apply period-dependent corrections,
//     then write back corrected (px,py,pz).
//
// Period dependence
// -----------------
// The coefficients used for the corrections depend on the run number.
// The mapping from run -> (coefficients) is defined in Config.yaml.
// Config.yaml provides one or more "periods" blocks, each with:
//   - run_ranges: list of [min,max] run-number ranges
//   - FD coefficients
//   - CD coefficients
//
 
#include "iguana/algorithms/Algorithm.h"
#include "iguana/algorithms/TypeDefs.h"
 
#include <map>
#include <string>
#include <utility>
#include <vector>
 
namespace iguana::clas12::rga {

 
  class ProtonEnergyLossCorrection : public Algorithm {
    DEFINE_IGUANA_ALGORITHM(ProtonEnergyLossCorrection, clas12::rga::ProtonEnergyLossCorrection)
 
  public:
    Momentum3 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;

    bool Run(hipo::bank& rec, hipo::bank const& run) const;
 
  private:
    void ConfigHook() override;
    void StartHook(hipo::banklist& banks) override;
    bool RunHook(hipo::banklist& banks) const override;
    void StopHook() override;
 
    hipo::banklist::size_type m_b_rec_particle{};
    hipo::banklist::size_type m_b_run_config{};
 
    struct Poly {
      std::vector<double> c;
    };
 
    struct RegionCoeffs {
      Poly A_p;
      Poly B_p;
      Poly C_p;
 
      Poly A_theta;
      Poly B_theta;
      Poly C_theta;
 
      Poly A_phi;
      Poly B_phi;
      Poly C_phi;
    };
 
    struct PeriodDef {
      std::vector<std::pair<int, int>> run_ranges;
      RegionCoeffs fd;
      RegionCoeffs cd;
    };
 
    std::map<std::string, PeriodDef> m_periods{};
 
    static bool IsFD(int status);
    static bool IsCD(int status);
 
    static double EvalPoly(Poly const& p, double x);
 
    static double Pmag(double px, double py, double pz);
    static double ThetaDeg(double px, double py, double pz);
    static double PhiDeg(double px, double py);
 
    static void SphericalToCartesian(double p, double theta_deg, double phi_deg,
      double& px, double& py, double& pz);
 
    PeriodDef const* FindPeriod(int run) const;
  };
 
}