From dd4ec490ae6741eb2930e7673999984067e31965 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Sat, 27 Jul 2024 16:39:03 -0400 Subject: [PATCH] feat: use InclusiveKinematics algorithms from EICrecon --- .../components/InclusiveKinematicsTruth.cpp | 161 +---------- .../src/components/InclusiveKinematicsDA.cpp | 241 +---------------- .../InclusiveKinematicsElectron.cpp | 249 +----------------- .../src/components/InclusiveKinematicsJB.cpp | 231 +--------------- .../components/InclusiveKinematicsSigma.cpp | 234 +--------------- .../components/InclusiveKinematicseSigma.cpp | 243 +---------------- 6 files changed, 36 insertions(+), 1323 deletions(-) diff --git a/JugFast/src/components/InclusiveKinematicsTruth.cpp b/JugFast/src/components/InclusiveKinematicsTruth.cpp index 822e8c1..720b93b 100644 --- a/JugFast/src/components/InclusiveKinematicsTruth.cpp +++ b/JugFast/src/components/InclusiveKinematicsTruth.cpp @@ -1,158 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" -#include -#include +#include +#include -#include "JugBase/IParticleSvc.h" -#include +namespace algorithms { -#include "JugBase/Utilities/Beam.h" +const std::shared_ptr ParticleSvc::kParticleMap = + std::make_shared(ParticleSvc::ParticleMap{ + { 0, { 0, 0, 0.0 , "unknown" }}, + }); -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -#include "edm4hep/utils/vector_utils.h" - -namespace Jug::Fast { - -class InclusiveKinematicsTruth : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - SmartIF m_pidSvc; - double m_proton{0}; - double m_neutron{0}; - double m_electron{0}; - -public: - InclusiveKinematicsTruth(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsTruth"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) { - return StatusCode::FAILURE; - } - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // Loop over generated particles to get incoming electron and proton beams - // and the scattered electron. In the presence of QED radition on the incoming - // or outgoing electron line, the vertex kinematics will be different than the - // kinematics calculated using the scattered electron as done here. - // Also need to update for CC events. - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ei_p = ei_coll[0].getMomentum(); - const auto ei_p_mag = edm4hep::utils::magnitude(ei_p); - const auto ei_mass = m_electron; - const PxPyPzEVector ei(ei_p.x, ei_p.y, ei_p.z, std::hypot(ei_p_mag, ei_mass)); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto pi_p = pi_coll[0].getMomentum(); - const auto pi_p_mag = edm4hep::utils::magnitude(pi_p); - const auto pi_mass = pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron; - const PxPyPzEVector pi(pi_p.x, pi_p.y, pi_p.z, std::hypot(pi_p_mag, pi_mass)); - - // Get first scattered electron - // Scattered electron. Currently taken as first status==1 electron in HEPMC record, - // which seems to be correct based on a cursory glance at the Pythia8 output. In the future, - // it may be better to trace back each final-state electron and see which one originates from - // the beam. - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_p = ef_coll[0].getMomentum(); - const auto ef_p_mag = edm4hep::utils::magnitude(ef_p); - const auto ef_mass = m_electron; - const PxPyPzEVector ef(ef_p.x, ef_p.y, ef_p.z, std::hypot(ef_p_mag, ef_mass)); - - // DIS kinematics calculations - const auto q = ei - ef; - const auto q_dot_pi = q.Dot(pi); - const auto Q2 = -q.Dot(q); - const auto y = q_dot_pi / ei.Dot(pi); - const auto nu = q_dot_pi / pi_mass; - const auto x = Q2 / (2.*q_dot_pi); - const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2); - const auto kin = out_kinematics.create(x, Q2, W, y, nu); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "ef = " << ef << endmsg; - debug() << "q = " << q << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +} // namespace algorithms // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicsTruth) - -} // namespace Jug::Fast +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicsTruth, eicrecon::InclusiveKinematicsTruth, Jug::Fast) diff --git a/JugReco/src/components/InclusiveKinematicsDA.cpp b/JugReco/src/components/InclusiveKinematicsDA.cpp index af63a5c..e05460a 100644 --- a/JugReco/src/components/InclusiveKinematicsDA.cpp +++ b/JugReco/src/components/InclusiveKinematicsDA.cpp @@ -1,238 +1,17 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" -#include -#include +#include +#include -#include "JugBase/IParticleSvc.h" -#include +namespace algorithms { -#include "JugBase/Utilities/Beam.h" -#include "JugBase/Utilities/Boost.h" +const std::shared_ptr ParticleSvc::kParticleMap = + std::make_shared(ParticleSvc::ParticleMap{ + { 0, { 0, 0, 0.0 , "unknown" }}, + }); -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -namespace Jug::Reco { - -class InclusiveKinematicsDA : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleCollection{ - "inputReconstructedParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleAssociation{ - "inputParticleAssociations", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - Gaudi::Property m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian}; - - SmartIF m_pidSvc; - double m_proton{0}, m_neutron{0}, m_electron{0}; - -public: - InclusiveKinematicsDA(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles"); - declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsDA"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) - return StatusCode::FAILURE; - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - const auto& rcparts = *(m_inputParticleCollection.get()); - const auto& rcassoc = *(m_inputParticleAssociation.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector ei( - Jug::Base::Beam::round_beam_four_momentum( - ei_coll[0].getMomentum(), - m_electron, - {-5.0, -10.0, -18.0}, - 0.0) - ); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector pi( - Jug::Base::Beam::round_beam_four_momentum( - pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, - {41.0, 100.0, 275.0}, - m_crossingAngle) - ); - - // Get first scattered electron - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc.begin(), - // rcassoc.end(), - // [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; }); - auto ef_assoc = rcassoc.begin(); - for (; ef_assoc != rcassoc.end(); ++ef_assoc) { - if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) { - break; - } - } - if (!(ef_assoc != rcassoc.end())) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Truth scattered electron not in reconstructed particles" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to caculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - //Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - double theta_e = 0; - - // Get boost to colinear frame - auto boost = Jug::Base::Boost::determine_boost(ei, pi); - - for(const auto& p : rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab); - - theta_e = e_boosted.Theta(); - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } - - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto ptsum = sqrt(pxsum*pxsum + pysum*pysum); - auto theta_h = 2.*atan(sigma_h/ptsum); - - // If no scattered electron was found - if (sigma_h <= 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - - // Calculate kinematic variables - const auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) ); - const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) ); - const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da); - const auto nu_da = Q2_da / (2.*m_proton*x_da); - const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da); - auto kin = out_kinematics.create(x_da, Q2_da, W_da, y_da, nu_da); - kin.setScat(ef_rc); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +} // namespace algorithms // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicsDA) - -} // namespace Jug::Reco +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicsDA, eicrecon::InclusiveKinematicsDA, Jug::Reco) diff --git a/JugReco/src/components/InclusiveKinematicsElectron.cpp b/JugReco/src/components/InclusiveKinematicsElectron.cpp index 928a5ae..4ee82ab 100644 --- a/JugReco/src/components/InclusiveKinematicsElectron.cpp +++ b/JugReco/src/components/InclusiveKinematicsElectron.cpp @@ -1,249 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" -#include -#include - -#include "JugBase/IParticleSvc.h" -#include - -#include "JugBase/Utilities/Beam.h" -#include "JugBase/Utilities/Boost.h" - -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -namespace Jug::Reco { - -class InclusiveKinematicsElectron : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleCollection{ - "inputReconstructedParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleAssociation{ - "inputParticleAssociations", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - Gaudi::Property m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian}; - - SmartIF m_pidSvc; - double m_proton{0}, m_neutron{0}, m_electron{0}; - -public: - InclusiveKinematicsElectron(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles"); - declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsElectron"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) - return StatusCode::FAILURE; - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - const auto& rcparts = *(m_inputParticleCollection.get()); - const auto& rcassoc = *(m_inputParticleAssociation.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // 1. find_if - //const auto mc_first_electron = std::find_if( - // mcparts.begin(), - // mcparts.end(), - // [](const auto& p){ return p.getPDG() == 11; }); - - // 2a. simple loop over iterator (post-increment) - //auto mc_first_electron = mcparts.end(); - //for (auto p = mcparts.begin(); p != mcparts.end(); p++) { - // if (p->getPDG() == 11) { - // mc_first_electron = p; - // break; - // } - //} - // 2b. simple loop over iterator (pre-increment) - //auto mc_first_electron = mcparts.end(); - //for (auto p = mcparts.begin(); p != mcparts.end(); ++p) { - // if (p->getPDG() == 11) { - // mc_first_electron = p; - // break; - // } - //} - - // 3. pre-initialized simple loop - //auto mc_first_electron = mcparts.begin(); - //for (; mc_first_electron != mcparts.end(); ++mc_first_electron) { - // if (mc_first_electron->getPDG() == 11) { - // break; - // } - //} - - // 4a. iterator equality - //if (mc_first_electron == mcparts.end()) { - // debug() << "No electron found" << endmsg; - // return StatusCode::FAILURE; - //} - // 4b. iterator inequality - //if (!(mc_first_electron != mcparts.end())) { - // debug() << "No electron found" << endmsg; - // return StatusCode::FAILURE; - //} - - // 5. ranges and views - //auto is_electron = [](const auto& p){ return p.getPDG() == 11; }; - //for (const auto& e: mcparts | std::views::filter(is_electron)) { - // break; - //} - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector ei( - Jug::Base::Beam::round_beam_four_momentum( - ei_coll[0].getMomentum(), - m_electron, - {-5.0, -10.0, -18.0}, - 0.0) - ); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector pi( - Jug::Base::Beam::round_beam_four_momentum( - pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, - {41.0, 100.0, 275.0}, - m_crossingAngle) - ); - - // Get first scattered electron - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc.begin(), - // rcassoc.end(), - // [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; }); - auto ef_assoc = rcassoc.begin(); - for (; ef_assoc != rcassoc.end(); ++ef_assoc) { - if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) { - break; - } - } - if (!(ef_assoc != rcassoc.end())) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Truth scattered electron not in reconstructed particles" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get outgoing scattered electron - // Use the true scattered electron from the MC information - std::vector electrons; - for (const auto& p: rcparts) { - if (p.getObjectID().index == ef_rc_id) { - electrons.emplace_back(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - break; - } - } - - // If no scattered electron was found - if (electrons.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - - // DIS kinematics calculations - const auto ef = electrons.front(); - const auto q = ei - ef; - const auto q_dot_pi = q.Dot(pi); - const auto Q2 = -q.Dot(q); - const auto y = q_dot_pi / ei.Dot(pi); - const auto nu = q_dot_pi / m_proton; - const auto x = Q2 / (2. * q_dot_pi); - const auto W = sqrt( + 2.*q_dot_pi - Q2); - auto kin = out_kinematics.create(x, Q2, W, y, nu); - kin.setScat(ef_rc); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "ef = " << ef << endmsg; - debug() << "q = " << q << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicsElectron) - -} // namespace Jug::Reco +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicsElectron, eicrecon::InclusiveKinematicsElectron, Jug::Reco) diff --git a/JugReco/src/components/InclusiveKinematicsJB.cpp b/JugReco/src/components/InclusiveKinematicsJB.cpp index f0f64c7..f794a68 100644 --- a/JugReco/src/components/InclusiveKinematicsJB.cpp +++ b/JugReco/src/components/InclusiveKinematicsJB.cpp @@ -1,231 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/PhysicalConstants.h" -#include "GaudiKernel/RndmGenerators.h" -#include -#include - -#include "JugBase/IParticleSvc.h" -#include - -#include "JugBase/Utilities/Beam.h" -#include "JugBase/Utilities/Boost.h" - -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -namespace Jug::Reco { - -class InclusiveKinematicsJB : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleCollection{ - "inputReconstructedParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleAssociation{ - "inputParticleAssociations", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - Gaudi::Property m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian}; - - SmartIF m_pidSvc; - double m_proton{0}, m_neutron{0}, m_electron{0}; - -public: - InclusiveKinematicsJB(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles"); - declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsJB"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) - return StatusCode::FAILURE; - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - const auto& rcparts = *(m_inputParticleCollection.get()); - const auto& rcassoc = *(m_inputParticleAssociation.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector ei( - Jug::Base::Beam::round_beam_four_momentum( - ei_coll[0].getMomentum(), - m_electron, - {-5.0, -10.0, -18.0}, - 0.0) - ); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector pi( - Jug::Base::Beam::round_beam_four_momentum( - pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, - {41.0, 100.0, 275.0}, - m_crossingAngle) - ); - - // Get first scattered electron - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc.begin(), - // rcassoc.end(), - // [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; }); - auto ef_assoc = rcassoc.begin(); - for (; ef_assoc != rcassoc.end(); ++ef_assoc) { - if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) { - break; - } - } - if (!(ef_assoc != rcassoc.end())) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Truth scattered electron not in reconstructed particles" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles other than the scattered electron - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to caculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - // Sums in colinear frame - double pxsum = 0; - double pysum = 0; - double pzsum = 0; - double Esum = 0; - - // Get boost to colinear frame - auto boost = Jug::Base::Boost::determine_boost(ei, pi); - - for (const auto& p: rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - - // Sum over all particles other than scattered electron - } else { - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab); - - pxsum += hf_boosted.Px(); - pysum += hf_boosted.Py(); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } - - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto ptsum = sqrt(pxsum*pxsum + pysum*pysum); - - // Sigma zero or negative - if (sigma_h <= 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Sigma zero or negative" << endmsg; - } - return StatusCode::SUCCESS; - } - - // Calculate kinematic variables - const auto y_jb = sigma_h / (2.*ei.energy()); - const auto Q2_jb = ptsum*ptsum / (1. - y_jb); - const auto x_jb = Q2_jb / (4.*ei.energy()*pi.energy()*y_jb); - const auto nu_jb = Q2_jb / (2.*m_proton*x_jb); - const auto W_jb = sqrt(m_proton*m_proton + 2*m_proton*nu_jb - Q2_jb); - auto kin = out_kinematics.create(x_jb, Q2_jb, W_jb, y_jb, nu_jb); - kin.setScat(ef_rc); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicsJB) - -} // namespace Jug::Reco +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicsJB, eicrecon::InclusiveKinematicsJB, Jug::Reco) diff --git a/JugReco/src/components/InclusiveKinematicsSigma.cpp b/JugReco/src/components/InclusiveKinematicsSigma.cpp index 59634a8..5244211 100644 --- a/JugReco/src/components/InclusiveKinematicsSigma.cpp +++ b/JugReco/src/components/InclusiveKinematicsSigma.cpp @@ -1,234 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck, Barak Schmookler +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" -#include -#include - -#include "JugBase/IParticleSvc.h" -#include - -#include "JugBase/Utilities/Beam.h" -#include "JugBase/Utilities/Boost.h" - -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -namespace Jug::Reco { - -class InclusiveKinematicsSigma : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleCollection{ - "inputReconstructedParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleAssociation{ - "inputParticleAssociations", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - Gaudi::Property m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian}; - - SmartIF m_pidSvc; - double m_proton{0}, m_neutron{0}, m_electron{0}; - -public: - InclusiveKinematicsSigma(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles"); - declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsSigma"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) - return StatusCode::FAILURE; - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - const auto& rcparts = *(m_inputParticleCollection.get()); - const auto& rcassoc = *(m_inputParticleAssociation.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector ei( - Jug::Base::Beam::round_beam_four_momentum( - ei_coll[0].getMomentum(), - m_electron, - {-5.0, -10.0, -18.0}, - 0.0) - ); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector pi( - Jug::Base::Beam::round_beam_four_momentum( - pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, - {41.0, 100.0, 275.0}, - m_crossingAngle) - ); - - // Get first scattered electron - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc.begin(), - // rcassoc.end(), - // [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; }); - auto ef_assoc = rcassoc.begin(); - for (; ef_assoc != rcassoc.end(); ++ef_assoc) { - if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) { - break; - } - } - if (!(ef_assoc != rcassoc.end())) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Truth scattered electron not in reconstructed particles" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to caculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - // Sums in colinear frame - double pzsum = 0; - double Esum = 0; - - double pt_e = 0; - double sigma_e = 0; - - // Get boost to colinear frame - auto boost = Jug::Base::Boost::determine_boost(ei, pi); - - for(const auto& p: rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab); - - pt_e = e_boosted.Pt(); - sigma_e = e_boosted.E() - e_boosted.Pz(); - - // Sum over all particles other than scattered electron - } else{ - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab); - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } - - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto sigma_tot = sigma_e + sigma_h; - - if (sigma_h <= 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No scattered electron found or sigma zero or negative" << endmsg; - } - return StatusCode::SUCCESS; - } - - // Calculate kinematic variables - const auto y_sig = sigma_h / sigma_tot; - const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig); - const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig); - const auto nu_sig = Q2_sig / (2.*m_proton*x_sig); - const auto W_sig = sqrt(m_proton*m_proton + 2*m_proton*nu_sig - Q2_sig); - auto kin = out_kinematics.create(x_sig, Q2_sig, W_sig, y_sig, nu_sig); - kin.setScat(ef_rc); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicsSigma) - -} // namespace Jug::Reco +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicsSigma, eicrecon::InclusiveKinematicsSigma, Jug::Reco) diff --git a/JugReco/src/components/InclusiveKinematicseSigma.cpp b/JugReco/src/components/InclusiveKinematicseSigma.cpp index 5d674fc..87d3ae9 100644 --- a/JugReco/src/components/InclusiveKinematicseSigma.cpp +++ b/JugReco/src/components/InclusiveKinematicseSigma.cpp @@ -1,243 +1,8 @@ // SPDX-License-Identifier: LGPL-3.0-or-later -// Copyright (C) 2022 Wouter Deconinck, Barak Schmookler +// Copyright (C) 2024 Wouter Deconinck -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Producer.h" -#include "GaudiAlg/Transformer.h" -#include "GaudiKernel/RndmGenerators.h" -#include "GaudiKernel/PhysicalConstants.h" -#include -#include - -#include "JugBase/IParticleSvc.h" -#include - -#include "JugBase/Utilities/Beam.h" -#include "JugBase/Utilities/Boost.h" - -#include "Math/Vector4D.h" -using ROOT::Math::PxPyPzEVector; - -// Event Model related classes -#include "edm4hep/MCParticleCollection.h" -#include "edm4eic/MCRecoParticleAssociationCollection.h" -#include "edm4eic/ReconstructedParticleCollection.h" -#include "edm4eic/InclusiveKinematicsCollection.h" - -namespace Jug::Reco { - -class InclusiveKinematicseSigma : public GaudiAlgorithm { -private: - DataHandle m_inputMCParticleCollection{ - "inputMCParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleCollection{ - "inputReconstructedParticles", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_inputParticleAssociation{ - "inputParticleAssociations", - Gaudi::DataHandle::Reader, - this}; - DataHandle m_outputInclusiveKinematicsCollection{ - "outputInclusiveKinematics", - Gaudi::DataHandle::Writer, - this}; - - Gaudi::Property m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian}; - - SmartIF m_pidSvc; - double m_proton{0}, m_neutron{0}, m_electron{0}; - -public: - InclusiveKinematicseSigma(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc) { - declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles"); - declareProperty("inputReconstructedParticles", m_inputParticleCollection, "ReconstructedParticles"); - declareProperty("inputParticleAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation"); - declareProperty("outputInclusiveKinematics", m_outputInclusiveKinematicsCollection, "InclusiveKinematicseSigma"); - } - - StatusCode initialize() override { - if (GaudiAlgorithm::initialize().isFailure()) - return StatusCode::FAILURE; - - m_pidSvc = service("ParticleSvc"); - if (!m_pidSvc) { - error() << "Unable to locate Particle Service. " - << "Make sure you have ParticleSvc in the configuration." - << endmsg; - return StatusCode::FAILURE; - } - m_proton = m_pidSvc->particle(2212).mass; - m_neutron = m_pidSvc->particle(2112).mass; - m_electron = m_pidSvc->particle(11).mass; - - return StatusCode::SUCCESS; - } - - StatusCode execute() override { - // input collections - const auto& mcparts = *(m_inputMCParticleCollection.get()); - const auto& rcparts = *(m_inputParticleCollection.get()); - const auto& rcassoc = *(m_inputParticleAssociation.get()); - // output collection - auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut()); - - // Get incoming electron beam - const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts); - if (ei_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector ei( - Jug::Base::Beam::round_beam_four_momentum( - ei_coll[0].getMomentum(), - m_electron, - {-5.0, -10.0, -18.0}, - 0.0) - ); - - // Get incoming hadron beam - const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts); - if (pi_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No beam hadron found" << endmsg; - } - return StatusCode::SUCCESS; - } - const PxPyPzEVector pi( - Jug::Base::Beam::round_beam_four_momentum( - pi_coll[0].getMomentum(), - pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron, - {41.0, 100.0, 275.0}, - m_crossingAngle) - ); - - // Get first scattered electron - const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts); - if (ef_coll.size() == 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No truth scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - // Associate first scattered electron with reconstructed electrons - //const auto ef_assoc = std::find_if( - // rcassoc.begin(), - // rcassoc.end(), - // [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].getObjectID().index; }); - auto ef_assoc = rcassoc.begin(); - for (; ef_assoc != rcassoc.end(); ++ef_assoc) { - if (ef_assoc->getSimID() == (unsigned) ef_coll[0].getObjectID().index) { - break; - } - } - if (!(ef_assoc != rcassoc.end())) { - if (msgLevel(MSG::DEBUG)) { - debug() << "Truth scattered electron not in reconstructed particles" << endmsg; - } - return StatusCode::SUCCESS; - } - const auto ef_rc{ef_assoc->getRec()}; - const auto ef_rc_id{ef_rc.getObjectID().index}; - - // Loop over reconstructed particles to get all outgoing particles - // ----------------------------------------------------------------- - // Right now, everything is taken from Reconstructed particles branches. - // - // This means the tracking detector is used for charged particles to caculate the momentum, - // and the magnitude of this momentum plus the true PID to calculate the energy. - // No requirement is made that these particles produce a hit in any other detector - // - // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the - // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the - // magnitude of the momentum. - // ----------------------------------------------------------------- - - //Sums in colinear frame - double pzsum = 0; - double Esum = 0; - - double pt_e = 0; - double sigma_e = 0; - - // Get boost to colinear frame - auto boost = Jug::Base::Boost::determine_boost(ei, pi); - - for (const auto& p: rcparts) { - // Get the scattered electron index and angle - if (p.getObjectID().index == ef_rc_id) { - // Lorentz vector in lab frame - PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab); - - pt_e = e_boosted.Pt(); - sigma_e = e_boosted.E() - e_boosted.Pz(); - - // Sum over all particles other than scattered electron - } else{ - // Lorentz vector in lab frame - PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); - // Boost to colinear frame - PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab); - - pzsum += hf_boosted.Pz(); - Esum += hf_boosted.E(); - } - } - - // DIS kinematics calculations - auto sigma_h = Esum - pzsum; - auto sigma_tot = sigma_e + sigma_h; - - // If no scattered electron was found - if (sigma_h <= 0) { - if (msgLevel(MSG::DEBUG)) { - debug() << "No scattered electron found" << endmsg; - } - return StatusCode::SUCCESS; - } - - // Calculate kinematic variables - const auto y_e = 1. - sigma_e / (2.*ei.energy()); - const auto Q2_e = (pt_e*pt_e) / (1. - y_e); - - const auto y_sig = sigma_h / sigma_tot; - const auto Q2_sig = (pt_e*pt_e) / (1. - y_sig); - const auto x_sig = Q2_sig / (4.*ei.energy()*pi.energy()*y_sig); - - const auto Q2_esig = Q2_e; - const auto x_esig = x_sig; - const auto y_esig = Q2_esig / (4.*ei.energy()*pi.energy()*x_esig); //equivalent to (2*ei.energy() / sigma_tot)*y_sig - const auto nu_esig = Q2_esig / (2.*m_proton*x_esig); - const auto W_esig = sqrt(m_proton*m_proton + 2*m_proton*nu_esig - Q2_esig); - auto kin = out_kinematics.create(x_esig, Q2_esig, W_esig, y_esig, nu_esig); - kin.setScat(ef_rc); - - // Debugging output - if (msgLevel(MSG::DEBUG)) { - debug() << "pi = " << pi << endmsg; - debug() << "ei = " << ei << endmsg; - debug() << "x,Q2,W,y,nu = " - << kin.getX() << "," - << kin.getQ2() << "," - << kin.getW() << "," - << kin.getY() << "," - << kin.getNu() - << endmsg; - } - - return StatusCode::SUCCESS; - } -}; +#include +#include // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) -DECLARE_COMPONENT(InclusiveKinematicseSigma) - -} // namespace Jug::Reco +JUGALGO_DEFINE_ALGORITHM(InclusiveKinematicseSigma, eicrecon::InclusiveKinematicseSigma, Jug::Reco)