From 4026ad13fdb679d55ad458fb8774d4d0d4b735b9 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Fri, 9 Feb 2024 20:27:26 -0500 Subject: [PATCH 1/9] feat: momentum corrections --- .../algorithms/clas12/MomentumCorrection.cc | 320 ++++++++++++++++++ .../algorithms/clas12/MomentumCorrection.h | 43 +++ src/iguana/algorithms/meson.build | 2 + 3 files changed, 365 insertions(+) create mode 100644 src/iguana/algorithms/clas12/MomentumCorrection.cc create mode 100644 src/iguana/algorithms/clas12/MomentumCorrection.h diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.cc b/src/iguana/algorithms/clas12/MomentumCorrection.cc new file mode 100644 index 00000000..d4c6cc15 --- /dev/null +++ b/src/iguana/algorithms/clas12/MomentumCorrection.cc @@ -0,0 +1,320 @@ +#include "MomentumCorrection.h" +#include + +namespace iguana::clas12 { + + REGISTER_IGUANA_ALGORITHM(MomentumCorrection); + + void MomentumCorrection::Start(hipo::banklist& banks) { + CacheBankIndex(banks, "REC::Particle", b_particle); + CacheBankIndex(banks, "RUN::config", b_config); + } + + + void MomentumCorrection::Run(hipo::banklist& banks) const { + auto& particleBank = GetBank(banks, b_particle, "REC::Particle"); + auto& configBank = GetBank(banks, b_config, "RUN::config"); + ShowBank(particleBank, Logger::Header("INPUT PARTICLES")); + + auto torus = configBank.getFloat("torus", 0); + + for(int row = 0; row < particleBank.getRows(); row++) { + + ////////////////////////////////// + ////////////////////////////////// + int sector = 1; // FIXME: get the true sector + ////////////////////////////////// + ////////////////////////////////// + + auto [px, py, pz] = Transform( + particleBank.getFloat("px", row), + particleBank.getFloat("py", row), + particleBank.getFloat("pz", row), + sector, + particleBank.getInt("pid", row), + torus + ); + particleBank.putFloat("px", row, px); + particleBank.putFloat("py", row, py); + particleBank.putFloat("pz", row, pz); + } + + ShowBank(particleBank, Logger::Header("OUTPUT PARTICLES")); + } + + + std::tuple MomentumCorrection::Transform(float px, float py, float pz, int sec, int pid, float torus) const { + auto cor = torus < 0 + ? CorrectionInbending(px, py, pz, sec, pid) + : CorrectionOutbending(px, py, pz, sec, pid); + return std::tuple{cor*px, cor*py, cor*pz}; + } + + + double MomentumCorrection::CorrectionInbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const { + + // skip the correction if it's not defined + if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 )) + return 1.0; + + // Momentum Magnitude + double pp = sqrt(Px*Px + Py*Py + Pz*Pz); + + // Initializing the correction factor + double dp = 0; + + // Defining Phi Angle + double Phi = (180/3.1415926)*atan2(Py, Px); + + // (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚) + if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){ + Phi += 360; + } + + // Getting Local Phi Angle + double PhiLocal = Phi - (sec - 1)*60; + + // Applying Shift Functions to Phi Angles (local shifted phi = phi) + double phi = PhiLocal; + + // For Electron Shift + if(pid == 11){ + phi = PhiLocal - 30/pp; + } + + // For π+ Pion/Proton Shift + if(pid == 211 || pid == 2212){ + phi = PhiLocal + (32/(pp-0.05)); + } + + // For π- Pion Shift + if(pid == -211){ + phi = PhiLocal - (32/(pp-0.05)); + } + + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //==================================================================================================================================// + //=======================//=======================// Electron Corrections //=======================//=======================// + //==================================================================================================================================// + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + if(pid == 11){ + if(sec == 1){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 1] is: + dp = ((-4.3303e-06)*phi*phi + (1.1006e-04)*phi + (-5.7235e-04))*pp*pp + ((3.2555e-05)*phi*phi + (-0.0014559)*phi + (0.0014878))*pp + ((-1.9577e-05)*phi*phi + (0.0017996)*phi + (0.025963)); + } + if(sec == 2){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 2] is: + dp = ((-9.8045e-07)*phi*phi + (6.7395e-05)*phi + (-4.6757e-05))*pp*pp + ((-1.4958e-05)*phi*phi + (-0.0011191)*phi + (-0.0025143))*pp + ((1.2699e-04)*phi*phi + (0.0033121)*phi + (0.020819)); + } + if(sec == 3){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 3] is: + dp = ((-5.9459e-07)*phi*phi + (-2.8289e-05)*phi + (-4.3541e-04))*pp*pp + ((-1.5025e-05)*phi*phi + (5.7730e-04)*phi + (-0.0077582))*pp + ((7.3348e-05)*phi*phi + (-0.001102)*phi + (0.057052)); + } + if(sec == 4){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 4] is: + dp = ((-2.2714e-06)*phi*phi + (-3.0360e-05)*phi + (-8.9322e-04))*pp*pp + ((2.9737e-05)*phi*phi + (5.1142e-04)*phi + (0.0045641))*pp + ((-1.0582e-04)*phi*phi + (-5.6852e-04)*phi + (0.027506)); + } + if(sec == 5){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 5] is: + dp = ((-1.1490e-06)*phi*phi + (-6.2147e-06)*phi + (-4.7235e-04))*pp*pp + ((3.7039e-06)*phi*phi + (-1.5943e-04)*phi + (-8.5238e-04))*pp + ((4.4069e-05)*phi*phi + (0.0014152)*phi + (0.031933)); + } + if(sec == 6){ + // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 6] is: + dp = ((1.1076e-06)*phi*phi + (4.0156e-05)*phi + (-1.6341e-04))*pp*pp + ((-2.8613e-05)*phi*phi + (-5.1861e-04)*phi + (-0.0056437))*pp + ((1.2419e-04)*phi*phi + (4.9084e-04)*phi + (0.049976)); + } + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //======================//======================// Electron Corrections (End) //======================//======================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //=========================//=========================// π+ Corrections //=========================//=========================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if(pid == 211){ + if(sec == 1){ + dp = ((-5.4904e-07)*phi*phi + (-1.4436e-05)*phi + (3.1534e-04))*pp*pp + ((3.8231e-06)*phi*phi + (3.6582e-04)*phi + (-0.0046759))*pp + ((-5.4913e-06)*phi*phi + (-4.0157e-04)*phi + (0.010767)); + dp = dp + ((6.1103e-07)*phi*phi + (5.5291e-06)*phi + (-1.9120e-04))*pp*pp + ((-3.2300e-06)*phi*phi + (1.5377e-05)*phi + (7.5279e-04))*pp + ((2.1434e-06)*phi*phi + (-6.9572e-06)*phi + (-7.9333e-05)); + dp = dp + ((-1.3049e-06)*phi*phi + (1.1295e-05)*phi + (4.5797e-04))*pp*pp + ((9.3122e-06)*phi*phi + (-5.1074e-05)*phi + (-0.0030757))*pp + ((-1.3102e-05)*phi*phi + (2.2153e-05)*phi + (0.0040938)); + } + if(sec == 2){ + dp = ((-1.0087e-06)*phi*phi + (2.1319e-05)*phi + (7.8641e-04))*pp*pp + ((6.7485e-06)*phi*phi + (7.3716e-05)*phi + (-0.0094591))*pp + ((-1.1820e-05)*phi*phi + (-3.8103e-04)*phi + (0.018936)); + dp = dp + ((8.8155e-07)*phi*phi + (-2.8257e-06)*phi + (-2.6729e-04))*pp*pp + ((-5.4499e-06)*phi*phi + (3.8397e-05)*phi + (0.0015914))*pp + ((6.8926e-06)*phi*phi + (-5.9386e-05)*phi + (-0.0021749)); + dp = dp + ((-2.0147e-07)*phi*phi + (1.1061e-05)*phi + (3.8827e-04))*pp*pp + ((4.9294e-07)*phi*phi + (-6.0257e-05)*phi + (-0.0022087))*pp + ((9.8548e-07)*phi*phi + (5.9047e-05)*phi + (0.0022905)); + } + if(sec == 3){ + dp = ((8.6722e-08)*phi*phi + (-1.7975e-05)*phi + (4.8118e-05))*pp*pp + ((2.6273e-06)*phi*phi + (3.1453e-05)*phi + (-0.0015943))*pp + ((-6.4463e-06)*phi*phi + (-5.8990e-05)*phi + (0.0041703)); + dp = dp + ((9.6317e-07)*phi*phi + (-1.7659e-06)*phi + (-8.8318e-05))*pp*pp + ((-5.1346e-06)*phi*phi + (8.3318e-06)*phi + (3.7723e-04))*pp + ((3.9548e-06)*phi*phi + (-6.9614e-05)*phi + (2.1393e-04)); + dp = dp + ((5.6438e-07)*phi*phi + (8.1678e-06)*phi + (-9.4406e-05))*pp*pp + ((-3.9074e-06)*phi*phi + (-6.5174e-05)*phi + (5.4218e-04))*pp + ((6.3198e-06)*phi*phi + (1.0611e-04)*phi + (-4.5749e-04)); + } + if(sec == 4){ + dp = ((4.3406e-07)*phi*phi + (-4.9036e-06)*phi + (2.3064e-04))*pp*pp + ((1.3624e-06)*phi*phi + (3.2907e-05)*phi + (-0.0034872))*pp + ((-5.1017e-06)*phi*phi + (2.4593e-05)*phi + (0.0092479)); + dp = dp + ((6.0218e-07)*phi*phi + (-1.4383e-05)*phi + (-3.1999e-05))*pp*pp + ((-1.1243e-06)*phi*phi + (9.3884e-05)*phi + (-4.1985e-04))*pp + ((-1.8808e-06)*phi*phi + (-1.2222e-04)*phi + (0.0014037)); + dp = dp + ((-2.5490e-07)*phi*phi + (-8.5120e-07)*phi + (7.9109e-05))*pp*pp + ((2.5879e-06)*phi*phi + (8.6108e-06)*phi + (-5.1533e-04))*pp + ((-4.4521e-06)*phi*phi + (-1.7012e-05)*phi + (7.4848e-04)); + } + if(sec == 5){ + dp = ((2.4292e-07)*phi*phi + (8.8741e-06)*phi + (2.9482e-04))*pp*pp + ((3.7229e-06)*phi*phi + (7.3215e-06)*phi + (-0.0050685))*pp + ((-1.1974e-05)*phi*phi + (-1.3043e-04)*phi + (0.0078836)); + dp = dp + ((1.0867e-06)*phi*phi + (-7.7630e-07)*phi + (-4.4930e-05))*pp*pp + ((-5.6564e-06)*phi*phi + (-1.3417e-05)*phi + (2.5224e-04))*pp + ((6.8460e-06)*phi*phi + (9.0495e-05)*phi + (-4.6587e-04)); + dp = dp + ((8.5720e-07)*phi*phi + (-6.7464e-06)*phi + (-4.0944e-05))*pp*pp + ((-4.7370e-06)*phi*phi + (5.8808e-05)*phi + (1.9047e-04))*pp + ((5.7404e-06)*phi*phi + (-1.1105e-04)*phi + (-1.9392e-04)); + } + if(sec == 6){ + dp = ((2.1191e-06)*phi*phi + (-3.3710e-05)*phi + (2.5741e-04))*pp*pp + ((-1.2915e-05)*phi*phi + (2.3753e-04)*phi + (-2.6882e-04))*pp + ((2.2676e-05)*phi*phi + (-2.3115e-04)*phi + (-0.001283)); + dp = dp + ((6.0270e-07)*phi*phi + (-6.8200e-06)*phi + (1.3103e-04))*pp*pp + ((-1.8745e-06)*phi*phi + (3.8646e-05)*phi + (-8.8056e-04))*pp + ((2.0885e-06)*phi*phi + (-3.4932e-05)*phi + (4.5895e-04)); + dp = dp + ((4.7349e-08)*phi*phi + (-5.7528e-06)*phi + (-3.4097e-06))*pp*pp + ((1.7731e-06)*phi*phi + (3.5865e-05)*phi + (-5.7881e-04))*pp + ((-9.7008e-06)*phi*phi + (-4.1836e-05)*phi + (0.0035403)); + } + + } + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //=========================//=========================// π+ Corrections (End) //=========================//=========================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //==================//==================// π- Corrections (Updated as of 01-13-2023) //==================//==================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if(pid == -211){ + if(sec == 1){ + + dp = (-9.2163E-07*phi*phi + 3.1862E-06*phi + 2.9805E-03)*pp*pp +(1.0435E-05*phi*phi + -8.7298E-05*phi + -1.7730E-02)*pp + -1.5154E-05*phi*phi + -1.3716E-04*phi + 2.2410E-02; + + } + if(sec == 2){ + + dp =(-1.9656E-06*phi*phi + 9.7389E-05*phi + 4.1250E-03)*pp*pp +(1.6439E-06*phi*phi + -4.6007E-04*phi + -1.9809E-02)*pp + 3.5794E-07*phi*phi + 4.8250E-04*phi + 1.7333E-02; + + } + + if(sec == 3){ + + dp =(2.5351E-06*phi*phi + 4.1043E-05*phi + 3.1157E-03)*pp*pp +(-1.3573E-05*phi*phi + -1.7609E-04*phi + -1.6759E-02)*pp + 1.4647E-05*phi*phi + 1.7484E-04*phi + 1.3805E-02; + + } + + if(sec == 4){ + + dp =(2.3500E-06*phi*phi + -7.7894E-05*phi + 4.4837E-03)*pp*pp +(-9.7915E-06*phi*phi + 4.6576E-04*phi + -2.6809E-02)*pp + 1.3819E-05*phi*phi + -5.6017E-04*phi + 3.0320E-02; + + } + + if(sec == 5){ + + dp =(-2.1809E-06*phi*phi + 2.4948E-05*phi + 2.7995E-03)*pp*pp +(6.3908E-06*phi*phi + -6.5122E-05*phi + -1.7571E-02)*pp + -1.9146E-06*phi*phi + -6.3799E-05*phi + 2.0877E-02; + + } + + if(sec == 6){ + + dp =(-9.3043E-06*phi*phi + 6.2678E-05*phi + 5.9660E-03)*pp*pp +(4.0581E-05*phi*phi + -3.0537E-04*phi + -3.1485E-02)*pp + -3.8345E-05*phi*phi + 2.0267E-04*phi + 3.3363E-02; + + + } + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //=======================//=======================// π- Corrections (End) //=======================//=======================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //====================================================================================================================================// + //=======================//=======================// All Proton Corrections //=======================//=======================// + //====================================================================================================================================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if(pid == 2212){ + if(sec == 1){ + dp = ((1 + std::copysign(1, (pp - 1.4)))/2)*((4.4034e-03)*pp + (-0.01703)) + ((1 + std::copysign(1, -(pp - 1.4)))/2)*((-0.10898)*(pp - 1.4)*(pp - 1.4) + (-0.09574)*(pp - 1.4) + ((4.4034e-03)*1.4 + (-0.01703))); + } + if(sec == 2){ + dp = ((1 + std::copysign(1, (pp - 1.5)))/2)*((0.01318)*pp + (-0.03403)) + ((1 + std::copysign(1, -(pp - 1.5)))/2)*((-0.09829)*(pp - 1.5)*(pp - 1.5) + (-0.0986)*(pp - 1.5) + ((0.01318)*1.5 + (-0.03403))); + } + if(sec == 3){ + dp = ((1 + std::copysign(1, (pp - 1.05)))/2)*((-4.7052e-03)*pp + (1.2410e-03)) + ((1 + std::copysign(1, -(pp - 1.05)))/2)*((-0.22721)*(pp - 1.05)*(pp - 1.05) + (-0.09702)*(pp - 1.05) + ((-4.7052e-03)*1.05 + (1.2410e-03))); + } + if(sec == 4){ + dp = ((1 + std::copysign(1, (pp - 1.4)))/2)*((-1.0900e-03)*pp + (-4.0573e-03)) + ((1 + std::copysign(1, -(pp - 1.4)))/2)*((-0.09236)*(pp - 1.4)*(pp - 1.4) + (-0.073)*(pp - 1.4) + ((-1.0900e-03)*1.4 + (-4.0573e-03))); + } + if(sec == 5){ + dp = ((1 + std::copysign(1, (pp - 1.5)))/2)*((7.3965e-03)*pp + (-0.02428)) + ((1 + std::copysign(1, -(pp - 1.5)))/2)*((-0.09539)*(pp - 1.5)*(pp - 1.5) + (-0.09263)*(pp - 1.5) + ((7.3965e-03)*1.5 + (-0.02428))); + } + if(sec == 6){ + dp = ((1 + std::copysign(1, (pp - 1.15)))/2)*((-7.6214e-03)*pp + (8.1014e-03)) + ((1 + std::copysign(1, -(pp - 1.15)))/2)*((-0.12718)*(pp - 1.15)*(pp - 1.15) + (-0.06626)*(pp - 1.15) + ((-7.6214e-03)*1.15 + (8.1014e-03))); + } + } + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=====================================================================================================================================// + //=======================//=======================// End of Proton Corrections //=======================//=======================// + //=====================================================================================================================================// + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + + return dp/pp; + } + + + /* + * INBENDING USAGE + + // The following code is for the Energy Loss Corrections for the proton + double dE_loss = 0; + // Inbending Energy Loss Correction // + if(proth < 27){ + dE_loss = exp(-2.739 - 3.932*pro) + 0.002907; + } + if(proth > 27){ + dE_loss = exp(-1.2 - 4.228*pro) + 0.007502; + } + double feloss = (pro + dE_loss)/pro; + + // Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code: + auto fe = dppC(ex, ey, ez, esec, 0) + 1; + auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1; + auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1; + auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1; + + auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0); + auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957); + auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957); + auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938); + + */ + + + + + double MomentumCorrection::CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const { + + // skip the correction if it's not defined + if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 )) + return 1.0; + + // + // TODO + // + return 1.0; + // + } + + void MomentumCorrection::Stop() { + } + +} diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.h b/src/iguana/algorithms/clas12/MomentumCorrection.h new file mode 100644 index 00000000..c0ef7148 --- /dev/null +++ b/src/iguana/algorithms/clas12/MomentumCorrection.h @@ -0,0 +1,43 @@ +#pragma once + +#include "iguana/algorithms/Algorithm.h" + +namespace iguana::clas12 { + + /// @brief Momentum Corrections + /// + /// Adapted from + class MomentumCorrection : public Algorithm { + + DEFINE_IGUANA_ALGORITHM(MomentumCorrection, clas12::MomentumCorrection) + + public: + + void Start(hipo::banklist& banks) override; + void Run(hipo::banklist& banks) const override; + void Stop() override; + + /// **Action function**: apply the momentum correction + /// @param px @f$p_x@f$ + /// @param py @f$p_y@f$ + /// @param pz @f$p_z@f$ + /// @param sec the sector + /// @param pid the particle PDG + /// @param torus torus setting + /// @returns the transformed momenta + std::tuple Transform(float px, float py, float pz, int sec, int pid, float torus) const; + + protected: + + double CorrectionInbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const; + double CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const; + + private: + + /// `hipo::banklist` index for the particle bank (as an example) + hipo::banklist::size_type b_particle; + hipo::banklist::size_type b_config; + + }; + +} diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 6a75e7b0..1b9766bd 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -7,6 +7,7 @@ algo_sources = [ 'clas12/EventBuilderFilter.cc', 'clas12/ZVertexFilter.cc', 'clas12/LorentzTransformer.cc', + 'clas12/MomentumCorrection.cc', ] # algorithm headers @@ -19,6 +20,7 @@ algo_public_headers = [ 'clas12/EventBuilderFilter.h', 'clas12/ZVertexFilter.h', 'clas12/LorentzTransformer.h', + 'clas12/MomentumCorrection.h', ] algo_config_files = [ From 8ba2e6196a733021c6ee0480b1e5f00bf39425cb Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 12 Feb 2024 15:49:20 -0500 Subject: [PATCH 2/9] feat: outbending corrections --- src/iguana/algorithms/TypeDefs.h | 2 + .../algorithms/clas12/MomentumCorrection.cc | 158 +++++++++++++++++- 2 files changed, 155 insertions(+), 5 deletions(-) diff --git a/src/iguana/algorithms/TypeDefs.h b/src/iguana/algorithms/TypeDefs.h index 25082fb2..e527e1da 100644 --- a/src/iguana/algorithms/TypeDefs.h +++ b/src/iguana/algorithms/TypeDefs.h @@ -11,4 +11,6 @@ namespace iguana { /// Generic Lorentz vector container type using lorentz_vector_t = std::tuple; + // TODO: enum for particle PDG + } diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.cc b/src/iguana/algorithms/clas12/MomentumCorrection.cc index d4c6cc15..62835817 100644 --- a/src/iguana/algorithms/clas12/MomentumCorrection.cc +++ b/src/iguana/algorithms/clas12/MomentumCorrection.cc @@ -307,13 +307,161 @@ namespace iguana::clas12 { if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 )) return 1.0; - // - // TODO - // - return 1.0; - // + // Momentum Magnitude + double pp = sqrt(Px*Px + Py*Py + Pz*Pz); + + // Initializing the correction factor + double dp = 0; + + // Defining Phi Angle + double Phi = (180/3.1415926)*atan2(Py, Px); + + // (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚) + if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){ + Phi += 360; + } + + // Getting Local Phi Angle + double PhiLocal = Phi - (sec - 1)*60; + + // Applying Shift Functions to Phi Angles (local shifted phi = phi) + double phi = PhiLocal; + // For Electron Shift + if(pid == 11){ + phi = PhiLocal - 30/pp; + } + // For π+ Pion/Proton Shift + if(pid == 211 || pid == 2212){ + phi = PhiLocal + (32/(pp-0.05)); + } + // For π- Pion Shift + if(pid == -211){ + phi = PhiLocal - (32/(pp-0.05)); + } + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=======================//=======================// Electron Corrections //=======================//=======================// + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if(pid == 11){ + if(sec == 1){ + dp = ((1.3189e-06)*phi*phi + (4.26057e-05)*phi + (-0.002322628))*pp*pp + ((-1.1409e-05)*phi*phi + (2.2188e-05)*phi + (0.02878927))*pp + ((2.4950e-05)*phi*phi + (1.6170e-06)*phi + (-0.061816275)); + } + if(sec == 2){ + dp = ((-2.9240e-07)*phi*phi + (3.2448e-07)*phi + (-0.001848308))*pp*pp + ((4.4500e-07)*phi*phi + (4.76324e-04)*phi + (0.02219469))*pp + ((6.9220e-06)*phi*phi + (-0.00153517)*phi + (-0.0479058)); + } + if(sec == 3){ + dp = ((2.71911e-06)*phi*phi + (1.657148e-05)*phi + (-0.001822211))*pp*pp + ((-4.96814e-05)*phi*phi + (-3.761117e-04)*phi + (0.02564148))*pp + ((1.97748e-04)*phi*phi + (9.58259e-04)*phi + (-0.05818292)); + } + if(sec == 4){ + dp = ((1.90966e-06)*phi*phi + (-2.4761e-05)*phi + (-0.00231562))*pp*pp + ((-2.3927e-05)*phi*phi + (2.25262e-04)*phi + (0.0291831))*pp + ((8.0515e-05)*phi*phi + (-6.42098e-04)*phi + (-0.06159197)); + } + if(sec == 5){ + dp = ((-3.6760323e-06)*phi*phi + (4.04398e-05)*phi + (-0.0021967515))*pp*pp + ((4.90857e-05)*phi*phi + (-4.37437e-04)*phi + (0.02494339))*pp + ((-1.08257e-04)*phi*phi + (0.00146111)*phi + (-0.0648485)); + } + if(sec == 6){ + dp = ((-6.2488e-08)*phi*phi + (2.23173e-05)*phi + (-0.00227522))*pp*pp + ((1.8372e-05)*phi*phi + (-7.5227e-05)*phi + (0.032636))*pp + ((-6.6566e-05)*phi*phi + (-2.4450e-04)*phi + (-0.072293)); + } + } + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //======================//======================// Electron Corrections (End) //======================//======================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=========================//=========================// π+ Corrections //=========================//=========================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if(pid == 211){ + if(sec == 1){ + dp = ((-1.7334e-06)*phi*phi + (1.45112e-05)*phi + (0.00150721))*pp*pp + ((6.6234e-06)*phi*phi + (-4.81191e-04)*phi + (-0.0138695))*pp + ((-3.23625e-06)*phi*phi + (2.79751e-04)*phi + (0.027726)); + } + if(sec == 2){ + dp = ((-4.475464e-06)*phi*phi + (-4.11573e-05)*phi + (0.00204557))*pp*pp + ((2.468278e-05)*phi*phi + (9.3590e-05)*phi + (-0.015399))*pp + ((-1.61547e-05)*phi*phi + (-2.4206e-04)*phi + (0.0231743)); + } + if(sec == 3){ + dp = ((-8.0374e-07)*phi*phi + (2.8728e-06)*phi + (0.00152163))*pp*pp + ((5.1347e-06)*phi*phi + (3.71709e-04)*phi + (-0.0165735))*pp + ((4.0105e-06)*phi*phi + (-5.289869e-04)*phi + (0.02175395)); + } + if(sec == 4){ + dp = ((-3.8790e-07)*phi*phi + (-4.78445e-05)*phi + (0.002324725))*pp*pp + ((6.80543e-06)*phi*phi + (5.69358e-04)*phi + (-0.0199162))*pp + ((-1.30264e-05)*phi*phi + (-5.91606e-04)*phi + (0.03202088)); + } + if(sec == 5){ + dp = ((2.198518e-06)*phi*phi + (-1.52535e-05)*phi + (0.001187761))*pp*pp + ((-1.000264e-05)*phi*phi + (1.63976e-04)*phi + (-0.01429673))*pp + ((9.4962e-06)*phi*phi + (-3.86691e-04)*phi + (0.0303695)); + } + if(sec == 6){ + dp = ((-3.92944e-07)*phi*phi + (1.45848e-05)*phi + (0.00120668))*pp*pp + ((3.7899e-06)*phi*phi + (-1.98219e-04)*phi + (-0.0131312))*pp + ((-3.9961e-06)*phi*phi + (-1.32883e-04)*phi + (0.0294497)); + } + } + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=======================//=======================// π+ Corrections (End) //=======================//=======================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=======================//=======================// π- Corrections (Start) //=======================//=======================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + if(pid == -211){ + if(sec == 1){ + dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02; + + + } + if(sec == 2){ + dp =(-4.6611E-06*phi*phi + -8.1637E-05*phi + 7.5013E-03)*pp*pp +(1.7616E-05*phi*phi + 3.5439E-04*phi + -3.7122E-02)*pp + -1.6286E-05*phi*phi + -2.6545E-04*phi + 4.5659E-02; + + } + + if(sec == 3){ + dp =(4.5270E-06*phi*phi + 2.2578E-04*phi + 5.9214E-03)*pp*pp +(-1.6419E-05*phi*phi + -8.1776E-04*phi + -3.2776E-02)*pp + 1.3734E-05*phi*phi + 6.6125E-04*phi + 4.5784E-02; + } + + if(sec == 4){ + dp =(-1.3141E-06*phi*phi + 1.9648E-04*phi + 7.6109E-03-0.006)*pp*pp +(8.0912E-06*phi*phi + -8.2672E-04*phi + -4.0495E-02+0.03)*pp + -3.1380E-06*phi*phi + 6.2211E-04*phi + 5.3361E-02-0.04; + + } + + if(sec == 5){ + dp =(-5.4065E-06*phi*phi + -1.6325E-05*phi + 1.2269E-02-0.002)*pp*pp +(1.9512E-05*phi*phi + 1.0228E-04*phi + -6.2351E-02+0.01)*pp + -9.5023E-06*phi*phi + -3.7997E-05*phi + 7.1061E-02-0.02; + + } + + if(sec == 6){ + dp =(-1.1882E-05*phi*phi + 2.0101E-04*phi + 1.1635E-02-0.01)*pp*pp +(5.8488E-05*phi*phi + -6.4709E-04*phi + -5.3833E-02+0.05)*pp + -4.4462E-05*phi*phi + 3.7529E-04*phi + 6.2130E-02-0.06; + + } + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //=======================//=======================// π- Corrections (End) //=======================//=======================// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + return dp/pp; } + + /* + * OUTBENDING USAGE + +// The following code is for the Energy Loss Corrections for the proton +double dE_loss = 0; +// Outbending Energy Loss Correction // +if(proth > 27){ + dE_loss = exp(-1.871 - 3.063*pro) + 0.007517; +} +double feloss = (pro + dE_loss)/pro; + +// Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code: +auto fe = dppC(ex, ey, ez, esec, 0) + 1; +auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1; +auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1; +auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1; + +auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0); +auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957); +auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957); +auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938); + +*/ + + void MomentumCorrection::Stop() { } From 9d4726bb3d111c934f4d15c8ce32a10b921d73fb Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 12 Feb 2024 15:53:38 -0500 Subject: [PATCH 3/9] doc: correction functions --- .../algorithms/clas12/MomentumCorrection.h | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.h b/src/iguana/algorithms/clas12/MomentumCorrection.h index c0ef7148..d732d5b5 100644 --- a/src/iguana/algorithms/clas12/MomentumCorrection.h +++ b/src/iguana/algorithms/clas12/MomentumCorrection.h @@ -27,15 +27,29 @@ namespace iguana::clas12 { /// @returns the transformed momenta std::tuple Transform(float px, float py, float pz, int sec, int pid, float torus) const; - protected: - + /// Calculate the correction factor for inbending data + /// @param Px @f$p_x@f$ + /// @param Py @f$p_y@f$ + /// @param Pz @f$p_z@f$ + /// @param sec the sector + /// @param pid the particle PDG + /// @returns the correction factor double CorrectionInbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const; + + /// Calculate the correction factor for outbending data + /// @param Px @f$p_x@f$ + /// @param Py @f$p_y@f$ + /// @param Pz @f$p_z@f$ + /// @param sec the sector + /// @param pid the particle PDG + /// @returns the correction factor double CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const; private: - /// `hipo::banklist` index for the particle bank (as an example) + /// `hipo::banklist` index for the particle bank hipo::banklist::size_type b_particle; + /// `hipo::banklist` index for the config bank hipo::banklist::size_type b_config; }; From ee5479a7bb3740b8263165cc1a7e516baa5078d3 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 12 Feb 2024 16:49:02 -0500 Subject: [PATCH 4/9] feat: energy loss corrections for proton --- .../algorithms/clas12/MomentumCorrection.cc | 199 ++++++------------ .../algorithms/clas12/MomentumCorrection.h | 16 ++ 2 files changed, 84 insertions(+), 131 deletions(-) diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.cc b/src/iguana/algorithms/clas12/MomentumCorrection.cc index 62835817..d74609be 100644 --- a/src/iguana/algorithms/clas12/MomentumCorrection.cc +++ b/src/iguana/algorithms/clas12/MomentumCorrection.cc @@ -44,10 +44,20 @@ namespace iguana::clas12 { std::tuple MomentumCorrection::Transform(float px, float py, float pz, int sec, int pid, float torus) const { - auto cor = torus < 0 - ? CorrectionInbending(px, py, pz, sec, pid) - : CorrectionOutbending(px, py, pz, sec, pid); - return std::tuple{cor*px, cor*py, cor*pz}; + // energy loss correction + auto e_cor = torus < 0 + ? EnergyLossInbending(px, py, pz, pid) + : EnergyLossOutbending(px, py, pz, pid); + // momentum correction + auto p_cor = torus < 0 + ? CorrectionInbending(e_cor * px, e_cor * py, e_cor * pz, sec, pid) + : CorrectionOutbending(e_cor * px, e_cor * py, e_cor * pz, sec, pid); + // return the corrected momentum + return std::tuple{ + e_cor * p_cor * px, + e_cor * p_cor * py, + e_cor * p_cor * pz + }; } @@ -64,7 +74,7 @@ namespace iguana::clas12 { double dp = 0; // Defining Phi Angle - double Phi = (180/3.1415926)*atan2(Py, Px); + double Phi = (180/M_PI)*atan2(Py, Px); // (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚) if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){ @@ -92,13 +102,11 @@ namespace iguana::clas12 { phi = PhiLocal - (32/(pp-0.05)); } - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //==================================================================================================================================// //=======================//=======================// Electron Corrections //=======================//=======================// //==================================================================================================================================// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 11){ if(sec == 1){ // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 1] is: @@ -126,14 +134,6 @@ namespace iguana::clas12 { } } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //====================================================================================================================================// - //======================//======================// Electron Corrections (End) //======================//======================// - //====================================================================================================================================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //====================================================================================================================================// //=========================//=========================// π+ Corrections //=========================//=========================// @@ -173,15 +173,6 @@ namespace iguana::clas12 { } - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //====================================================================================================================================// - //=========================//=========================// π+ Corrections (End) //=========================//=========================// - //====================================================================================================================================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //====================================================================================================================================// //==================//==================// π- Corrections (Updated as of 01-13-2023) //==================//==================// @@ -225,14 +216,6 @@ namespace iguana::clas12 { } } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //====================================================================================================================================// - //=======================//=======================// π- Corrections (End) //=======================//=======================// - //====================================================================================================================================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //====================================================================================================================================// //=======================//=======================// All Proton Corrections //=======================//=======================// @@ -259,52 +242,14 @@ namespace iguana::clas12 { } } - ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //=====================================================================================================================================// - //=======================//=======================// End of Proton Corrections //=======================//=======================// - //=====================================================================================================================================// - ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - - - return dp/pp; + return dp/pp + 1; } - /* - * INBENDING USAGE - - // The following code is for the Energy Loss Corrections for the proton - double dE_loss = 0; - // Inbending Energy Loss Correction // - if(proth < 27){ - dE_loss = exp(-2.739 - 3.932*pro) + 0.002907; - } - if(proth > 27){ - dE_loss = exp(-1.2 - 4.228*pro) + 0.007502; - } - double feloss = (pro + dE_loss)/pro; - - // Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code: - auto fe = dppC(ex, ey, ez, esec, 0) + 1; - auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1; - auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1; - auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1; - - auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0); - auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957); - auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957); - auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938); - - */ - - - - double MomentumCorrection::CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const { // skip the correction if it's not defined - if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 )) + if(!( pid == 11 || std::abs(pid) == 211)) return 1.0; // Momentum Magnitude @@ -314,7 +259,7 @@ namespace iguana::clas12 { double dp = 0; // Defining Phi Angle - double Phi = (180/3.1415926)*atan2(Py, Px); + double Phi = (180/M_PI)*atan2(Py, Px); // (Initial) Shift of the Phi Angle (done to realign sectors whose data is separated when plotted from ±180˚) if(((sec == 4 || sec == 3) && Phi < 0) || (sec > 4 && Phi < 90)){ @@ -362,9 +307,6 @@ namespace iguana::clas12 { dp = ((-6.2488e-08)*phi*phi + (2.23173e-05)*phi + (-0.00227522))*pp*pp + ((1.8372e-05)*phi*phi + (-7.5227e-05)*phi + (0.032636))*pp + ((-6.6566e-05)*phi*phi + (-2.4450e-04)*phi + (-0.072293)); } } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //======================//======================// Electron Corrections (End) //======================//======================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //=========================//=========================// π+ Corrections //=========================//=========================// @@ -389,77 +331,72 @@ namespace iguana::clas12 { dp = ((-3.92944e-07)*phi*phi + (1.45848e-05)*phi + (0.00120668))*pp*pp + ((3.7899e-06)*phi*phi + (-1.98219e-04)*phi + (-0.0131312))*pp + ((-3.9961e-06)*phi*phi + (-1.32883e-04)*phi + (0.0294497)); } } - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //=======================//=======================// π+ Corrections (End) //=======================//=======================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //=======================//=======================// π- Corrections (Start) //=======================//=======================// + //=======================//=======================// π- Corrections //=======================//=======================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == -211){ - if(sec == 1){ - dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02; - - - } - if(sec == 2){ - dp =(-4.6611E-06*phi*phi + -8.1637E-05*phi + 7.5013E-03)*pp*pp +(1.7616E-05*phi*phi + 3.5439E-04*phi + -3.7122E-02)*pp + -1.6286E-05*phi*phi + -2.6545E-04*phi + 4.5659E-02; - - } - - if(sec == 3){ - dp =(4.5270E-06*phi*phi + 2.2578E-04*phi + 5.9214E-03)*pp*pp +(-1.6419E-05*phi*phi + -8.1776E-04*phi + -3.2776E-02)*pp + 1.3734E-05*phi*phi + 6.6125E-04*phi + 4.5784E-02; - } - - if(sec == 4){ - dp =(-1.3141E-06*phi*phi + 1.9648E-04*phi + 7.6109E-03-0.006)*pp*pp +(8.0912E-06*phi*phi + -8.2672E-04*phi + -4.0495E-02+0.03)*pp + -3.1380E-06*phi*phi + 6.2211E-04*phi + 5.3361E-02-0.04; - - } + if(sec == 1){ + dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02; + } + if(sec == 2){ + dp =(-4.6611E-06*phi*phi + -8.1637E-05*phi + 7.5013E-03)*pp*pp +(1.7616E-05*phi*phi + 3.5439E-04*phi + -3.7122E-02)*pp + -1.6286E-05*phi*phi + -2.6545E-04*phi + 4.5659E-02; + } + if(sec == 3){ + dp =(4.5270E-06*phi*phi + 2.2578E-04*phi + 5.9214E-03)*pp*pp +(-1.6419E-05*phi*phi + -8.1776E-04*phi + -3.2776E-02)*pp + 1.3734E-05*phi*phi + 6.6125E-04*phi + 4.5784E-02; + } + if(sec == 4){ + dp =(-1.3141E-06*phi*phi + 1.9648E-04*phi + 7.6109E-03-0.006)*pp*pp +(8.0912E-06*phi*phi + -8.2672E-04*phi + -4.0495E-02+0.03)*pp + -3.1380E-06*phi*phi + 6.2211E-04*phi + 5.3361E-02-0.04; + } + if(sec == 5){ + dp =(-5.4065E-06*phi*phi + -1.6325E-05*phi + 1.2269E-02-0.002)*pp*pp +(1.9512E-05*phi*phi + 1.0228E-04*phi + -6.2351E-02+0.01)*pp + -9.5023E-06*phi*phi + -3.7997E-05*phi + 7.1061E-02-0.02; + } + if(sec == 6){ + dp =(-1.1882E-05*phi*phi + 2.0101E-04*phi + 1.1635E-02-0.01)*pp*pp +(5.8488E-05*phi*phi + -6.4709E-04*phi + -5.3833E-02+0.05)*pp + -4.4462E-05*phi*phi + 3.7529E-04*phi + 6.2130E-02-0.06; + } + } - if(sec == 5){ - dp =(-5.4065E-06*phi*phi + -1.6325E-05*phi + 1.2269E-02-0.002)*pp*pp +(1.9512E-05*phi*phi + 1.0228E-04*phi + -6.2351E-02+0.01)*pp + -9.5023E-06*phi*phi + -3.7997E-05*phi + 7.1061E-02-0.02; + return dp/pp + 1; + } - } - if(sec == 6){ - dp =(-1.1882E-05*phi*phi + 2.0101E-04*phi + 1.1635E-02-0.01)*pp*pp +(5.8488E-05*phi*phi + -6.4709E-04*phi + -5.3833E-02+0.05)*pp + -4.4462E-05*phi*phi + 3.7529E-04*phi + 6.2130E-02-0.06; + double MomentumCorrection::EnergyLossInbending(const float Px, const float Py, const float Pz, const int pid) const { - } - } + // The following code is for the Energy Loss Corrections for the proton + if(pid != 2212) + return 1.0; - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - //=======================//=======================// π- Corrections (End) //=======================//=======================// - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + double dE_loss = 0; + auto pro = sqrt(Px*Px + Py*Py + Pz*Pz); + auto proth = atan2(sqrt(Px*Px + Py*Py), Pz)*(180/M_PI); - return dp/pp; + // Inbending Energy Loss Correction // + if(proth < 27){ + dE_loss = exp(-2.739 - 3.932*pro) + 0.002907; + } + else { + dE_loss = exp(-1.2 - 4.228*pro) + 0.007502; + } + return (pro + dE_loss)/pro; } - /* - * OUTBENDING USAGE + double MomentumCorrection::EnergyLossOutbending(const float Px, const float Py, const float Pz, const int pid) const { -// The following code is for the Energy Loss Corrections for the proton -double dE_loss = 0; -// Outbending Energy Loss Correction // -if(proth > 27){ - dE_loss = exp(-1.871 - 3.063*pro) + 0.007517; -} -double feloss = (pro + dE_loss)/pro; - -// Below shows how the corrections are to be applied using the ROOT momentum 4-vector using the above code: -auto fe = dppC(ex, ey, ez, esec, 0) + 1; -auto fpip = dppC(pipx, pipy, pipz, pipsec, 1) + 1; -auto fpim = dppC(pimx, pimy, pimz, pimsec, 2) + 1; -auto fpro = dppC(prox*feloss, proy*feloss, proz*feloss, prosec, 3) + 1; + // The following code is for the Energy Loss Corrections for the proton + if(pid != 2212) + return 1.0; -auto eleC = ROOT::Math::PxPyPzMVector(ex*fe, ey*fe, ez*fe, 0); -auto pipC = ROOT::Math::PxPyPzMVector(pipx*fpip, pipy*fpip, pipz*fpip, 0.13957); -auto pimC = ROOT::Math::PxPyPzMVector(pimx*fpim, pimy*fpim, pimz*fpim, 0.13957); -auto proC = ROOT::Math::PxPyPzMVector(prox*feloss*fpro, proy*feloss*fpro, proz*feloss*fpro, 0.938); + double dE_loss = 0; + auto pro = sqrt(Px*Px + Py*Py + Pz*Pz); + auto proth = atan2(sqrt(Px*Px + Py*Py), Pz)*(180/M_PI); -*/ + // Outbending Energy Loss Correction // + if(proth > 27){ + dE_loss = exp(-1.871 - 3.063*pro) + 0.007517; + } + return (pro + dE_loss)/pro; + } void MomentumCorrection::Stop() { diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.h b/src/iguana/algorithms/clas12/MomentumCorrection.h index d732d5b5..e256834f 100644 --- a/src/iguana/algorithms/clas12/MomentumCorrection.h +++ b/src/iguana/algorithms/clas12/MomentumCorrection.h @@ -45,6 +45,22 @@ namespace iguana::clas12 { /// @returns the correction factor double CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const; + /// Energy loss correction for inbending data + /// @param Px @f$p_x@f$ + /// @param Py @f$p_y@f$ + /// @param Pz @f$p_z@f$ + /// @param pid the particle PDG + /// @returns the correction factor + double EnergyLossInbending(const float Px, const float Py, const float Pz, const int pid) const; + + /// Energy loss correction for outbending data + /// @param Px @f$p_x@f$ + /// @param Py @f$p_y@f$ + /// @param Pz @f$p_z@f$ + /// @param pid the particle PDG + /// @returns the correction factor + double EnergyLossOutbending(const float Px, const float Py, const float Pz, const int pid) const; + private: /// `hipo::banklist` index for the particle bank From 60dc540984d1794d43c4e1b429699b6b4946965d Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 12 Feb 2024 17:02:10 -0500 Subject: [PATCH 5/9] feat: PDG codes --- src/iguana/algorithms/TypeDefs.h | 13 ++++++- .../algorithms/clas12/MomentumCorrection.cc | 35 ++++++++++--------- 2 files changed, 30 insertions(+), 18 deletions(-) diff --git a/src/iguana/algorithms/TypeDefs.h b/src/iguana/algorithms/TypeDefs.h index e527e1da..aa29c4da 100644 --- a/src/iguana/algorithms/TypeDefs.h +++ b/src/iguana/algorithms/TypeDefs.h @@ -11,6 +11,17 @@ namespace iguana { /// Generic Lorentz vector container type using lorentz_vector_t = std::tuple; - // TODO: enum for particle PDG + /// Light-weight namespace for particle constants + namespace particle { + + /// PDG codes + enum PDG { + electron = 11, + pi_plus = 211, + pi_minus = -211, + proton = 2212 + }; + + } } diff --git a/src/iguana/algorithms/clas12/MomentumCorrection.cc b/src/iguana/algorithms/clas12/MomentumCorrection.cc index d74609be..f9ade0bd 100644 --- a/src/iguana/algorithms/clas12/MomentumCorrection.cc +++ b/src/iguana/algorithms/clas12/MomentumCorrection.cc @@ -1,4 +1,5 @@ #include "MomentumCorrection.h" +#include "iguana/algorithms/TypeDefs.h" #include namespace iguana::clas12 { @@ -64,7 +65,7 @@ namespace iguana::clas12 { double MomentumCorrection::CorrectionInbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const { // skip the correction if it's not defined - if(!( pid == 11 || std::abs(pid) == 211 || pid == 2212 )) + if(!( pid == particle::electron || pid == particle::pi_plus || pid == particle::pi_minus || pid == particle::proton )) return 1.0; // Momentum Magnitude @@ -88,17 +89,17 @@ namespace iguana::clas12 { double phi = PhiLocal; // For Electron Shift - if(pid == 11){ + if(pid == particle::electron){ phi = PhiLocal - 30/pp; } // For π+ Pion/Proton Shift - if(pid == 211 || pid == 2212){ + if(pid == particle::pi_plus || pid == particle::proton){ phi = PhiLocal + (32/(pp-0.05)); } // For π- Pion Shift - if(pid == -211){ + if(pid == particle::pi_minus){ phi = PhiLocal - (32/(pp-0.05)); } @@ -107,7 +108,7 @@ namespace iguana::clas12 { //=======================//=======================// Electron Corrections //=======================//=======================// //==================================================================================================================================// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 11){ + if(pid == particle::electron){ if(sec == 1){ // The CONTINUOUS QUADRATIC function predicted for ∆p_{El} for [Cor = Uncorrected][Sector 1] is: dp = ((-4.3303e-06)*phi*phi + (1.1006e-04)*phi + (-5.7235e-04))*pp*pp + ((3.2555e-05)*phi*phi + (-0.0014559)*phi + (0.0014878))*pp + ((-1.9577e-05)*phi*phi + (0.0017996)*phi + (0.025963)); @@ -139,7 +140,7 @@ namespace iguana::clas12 { //=========================//=========================// π+ Corrections //=========================//=========================// //====================================================================================================================================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 211){ + if(pid == particle::pi_plus){ if(sec == 1){ dp = ((-5.4904e-07)*phi*phi + (-1.4436e-05)*phi + (3.1534e-04))*pp*pp + ((3.8231e-06)*phi*phi + (3.6582e-04)*phi + (-0.0046759))*pp + ((-5.4913e-06)*phi*phi + (-4.0157e-04)*phi + (0.010767)); dp = dp + ((6.1103e-07)*phi*phi + (5.5291e-06)*phi + (-1.9120e-04))*pp*pp + ((-3.2300e-06)*phi*phi + (1.5377e-05)*phi + (7.5279e-04))*pp + ((2.1434e-06)*phi*phi + (-6.9572e-06)*phi + (-7.9333e-05)); @@ -178,7 +179,7 @@ namespace iguana::clas12 { //==================//==================// π- Corrections (Updated as of 01-13-2023) //==================//==================// //====================================================================================================================================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == -211){ + if(pid == particle::pi_minus){ if(sec == 1){ dp = (-9.2163E-07*phi*phi + 3.1862E-06*phi + 2.9805E-03)*pp*pp +(1.0435E-05*phi*phi + -8.7298E-05*phi + -1.7730E-02)*pp + -1.5154E-05*phi*phi + -1.3716E-04*phi + 2.2410E-02; @@ -221,7 +222,7 @@ namespace iguana::clas12 { //=======================//=======================// All Proton Corrections //=======================//=======================// //====================================================================================================================================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 2212){ + if(pid == particle::proton){ if(sec == 1){ dp = ((1 + std::copysign(1, (pp - 1.4)))/2)*((4.4034e-03)*pp + (-0.01703)) + ((1 + std::copysign(1, -(pp - 1.4)))/2)*((-0.10898)*(pp - 1.4)*(pp - 1.4) + (-0.09574)*(pp - 1.4) + ((4.4034e-03)*1.4 + (-0.01703))); } @@ -249,7 +250,7 @@ namespace iguana::clas12 { double MomentumCorrection::CorrectionOutbending(const float Px, const float Py, const float Pz, const int sec, const int pid) const { // skip the correction if it's not defined - if(!( pid == 11 || std::abs(pid) == 211)) + if(!( pid == particle::electron || pid == particle::pi_plus || pid == particle::pi_minus )) return 1.0; // Momentum Magnitude @@ -272,22 +273,22 @@ namespace iguana::clas12 { // Applying Shift Functions to Phi Angles (local shifted phi = phi) double phi = PhiLocal; // For Electron Shift - if(pid == 11){ + if(pid == particle::electron){ phi = PhiLocal - 30/pp; } // For π+ Pion/Proton Shift - if(pid == 211 || pid == 2212){ + if(pid == particle::pi_plus || pid == particle::proton){ phi = PhiLocal + (32/(pp-0.05)); } // For π- Pion Shift - if(pid == -211){ + if(pid == particle::pi_minus){ phi = PhiLocal - (32/(pp-0.05)); } ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //=======================//=======================// Electron Corrections //=======================//=======================// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 11){ + if(pid == particle::electron){ if(sec == 1){ dp = ((1.3189e-06)*phi*phi + (4.26057e-05)*phi + (-0.002322628))*pp*pp + ((-1.1409e-05)*phi*phi + (2.2188e-05)*phi + (0.02878927))*pp + ((2.4950e-05)*phi*phi + (1.6170e-06)*phi + (-0.061816275)); } @@ -311,7 +312,7 @@ namespace iguana::clas12 { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //=========================//=========================// π+ Corrections //=========================//=========================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == 211){ + if(pid == particle::pi_plus){ if(sec == 1){ dp = ((-1.7334e-06)*phi*phi + (1.45112e-05)*phi + (0.00150721))*pp*pp + ((6.6234e-06)*phi*phi + (-4.81191e-04)*phi + (-0.0138695))*pp + ((-3.23625e-06)*phi*phi + (2.79751e-04)*phi + (0.027726)); } @@ -335,7 +336,7 @@ namespace iguana::clas12 { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //=======================//=======================// π- Corrections //=======================//=======================// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - if(pid == -211){ + if(pid == particle::pi_minus){ if(sec == 1){ dp = (7.8044E-06*phi*phi + -9.4703E-05*phi + 4.6696E-03)*pp*pp +(-3.4668E-05*phi*phi + 6.2280E-04*phi + -2.4273E-02)*pp + 2.3566E-05*phi*phi + -5.8519E-04*phi + 3.9226E-02; } @@ -363,7 +364,7 @@ namespace iguana::clas12 { double MomentumCorrection::EnergyLossInbending(const float Px, const float Py, const float Pz, const int pid) const { // The following code is for the Energy Loss Corrections for the proton - if(pid != 2212) + if(pid != particle::proton) return 1.0; double dE_loss = 0; @@ -384,7 +385,7 @@ namespace iguana::clas12 { double MomentumCorrection::EnergyLossOutbending(const float Px, const float Py, const float Pz, const int pid) const { // The following code is for the Energy Loss Corrections for the proton - if(pid != 2212) + if(pid != particle::proton) return 1.0; double dE_loss = 0; From d5af10b6a87883ccfd4284839535a38a1e9c738a Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Tue, 20 Feb 2024 17:22:53 -0500 Subject: [PATCH 6/9] chore: codeowners --- CODEOWNERS | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CODEOWNERS b/CODEOWNERS index e3bed0fc..02dac547 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -2,6 +2,7 @@ # GitHub Handle Name # ============= ==== # @c-dilks Christopher Dilks +# @RichCap Richard Capobianco # @rtysonCLAS12 Richard Tyson ################################################################################# @@ -9,6 +10,8 @@ examples/config/ex2.yaml @rtysonCLAS12 examples/iguana-example-02-YAMLReader.cc @rtysonCLAS12 examples/iguana-example-03-zvertex-filter.cc @rtysonCLAS12 +src/iguana/algorithms/clas12/MomentumCorrection.cc @RichCap @c-dilks +src/iguana/algorithms/clas12/MomentumCorrection.h @RichCap @c-dilks src/iguana/algorithms/clas12/ZVertexFilter.cc @rtysonCLAS12 src/iguana/algorithms/clas12/ZVertexFilter.h @rtysonCLAS12 src/iguana/algorithms/clas12/ZVertexFilter.yaml @rtysonCLAS12 From 92c33305e0587c0f8146bdb4f2ed07c8a150fe1d Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Tue, 20 Feb 2024 17:41:54 -0500 Subject: [PATCH 7/9] fix: add to unit tests --- src/iguana/algorithms/meson.build | 1 + 1 file changed, 1 insertion(+) diff --git a/src/iguana/algorithms/meson.build b/src/iguana/algorithms/meson.build index 1b9766bd..99609a7b 100644 --- a/src/iguana/algorithms/meson.build +++ b/src/iguana/algorithms/meson.build @@ -33,6 +33,7 @@ algos_and_banks_for_unit_testing = { 'clas12::EventBuilderFilter': ['REC::Particle'], 'clas12::ZVertexFilter': ['REC::Particle'], 'clas12::LorentzTransformer': ['REC::Particle'], + 'clas12::MomentumCorrection': ['REC::Particle', 'RUN::config'], } algo_lib = shared_library( From 358610a1d8bacee7b89ce0ce71bdcfd72693f73e Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Tue, 20 Feb 2024 18:40:18 -0500 Subject: [PATCH 8/9] ci: UBSAN suppression --- .github/UBSAN.supp | 1 + .github/workflows/ci.yml | 3 +++ 2 files changed, 4 insertions(+) create mode 100644 .github/UBSAN.supp diff --git a/.github/UBSAN.supp b/.github/UBSAN.supp new file mode 100644 index 00000000..fdcac37e --- /dev/null +++ b/.github/UBSAN.supp @@ -0,0 +1 @@ +alignment:hipo::structure::getFloatAt diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 79bd41ac..07a8ac7b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,6 +42,9 @@ defaults: env: hipo_version: f40da676bbd1745398e9fdf233ff213ff98798f1 num_events: 10 + UBSAN_OPTIONS: halt_on_error=1:abort_on_error=1:print_summary=1:print_stacktrace=1:suppressions=${{ github.workspace }}/.github/UBSAN.supp + ASAN_OPTIONS: halt_on_error=1:abort_on_error=1:print_summary=1 + jobs: From 1590cf608483bb4d7e2dcee9558ce1d81f972672 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Tue, 20 Feb 2024 19:02:16 -0500 Subject: [PATCH 9/9] fix: use relative path for UBSAN suppressions --- .github/workflows/ci.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 07a8ac7b..d2dcf882 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,7 +42,8 @@ defaults: env: hipo_version: f40da676bbd1745398e9fdf233ff213ff98798f1 num_events: 10 - UBSAN_OPTIONS: halt_on_error=1:abort_on_error=1:print_summary=1:print_stacktrace=1:suppressions=${{ github.workspace }}/.github/UBSAN.supp + # sanitizer options; NOTE: suppression path assumes the build directory is a subdirectory of the top-level repo directory + UBSAN_OPTIONS: halt_on_error=1:abort_on_error=1:print_summary=1:print_stacktrace=1:suppressions=../.github/UBSAN.supp ASAN_OPTIONS: halt_on_error=1:abort_on_error=1:print_summary=1