diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc new file mode 100644 index 0000000000..96de76822e --- /dev/null +++ b/src/algorithms/reco/HadronicFinalState.cc @@ -0,0 +1,154 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Tyler Kutz + +#include +#if EDM4EIC_VERSION_MAJOR >= 6 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Beam.h" +#include "Boost.h" +#include "HadronicFinalState.h" + +using ROOT::Math::PxPyPzEVector; + +namespace eicrecon { + + void HadronicFinalState::init(std::shared_ptr& logger) { + m_log = logger; + // m_pidSvc = service("ParticleSvc"); + // if (!m_pidSvc) { + // m_log->debug("Unable to locate Particle Service. " + // "Make sure you have ParticleSvc in the configuration."); + // } + } + + void HadronicFinalState::process( + const HadronicFinalState::Input& input, + const HadronicFinalState::Output& output) const { + + const auto [mcparts, rcparts, rcassoc] = input; + auto [hadronicfinalstate] = output; + + // Get incoming electron beam + const auto ei_coll = find_first_beam_electron(mcparts); + if (ei_coll.size() == 0) { + m_log->debug("No beam electron found"); + return; + } + const PxPyPzEVector ei( + 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 = find_first_beam_hadron(mcparts); + if (pi_coll.size() == 0) { + m_log->debug("No beam hadron found"); + return; + } + const PxPyPzEVector pi( + 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 = find_first_scattered_electron(mcparts); + if (ef_coll.size() == 0) { + m_log->debug("No truth scattered electron found"); + return; + } + // 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.getSim().getObjectID() == ef_coll[0].getObjectID(); }); + auto ef_assoc = rcassoc->begin(); + for (; ef_assoc != rcassoc->end(); ++ef_assoc) { + if (ef_assoc->getSim().getObjectID() == ef_coll[0].getObjectID()) { + break; + } + } + if (!(ef_assoc != rcassoc->end())) { + m_log->debug("Truth scattered electron not in reconstructed particles"); + return; + } + const auto ef_rc{ef_assoc->getRec()}; + const auto ef_rc_id{ef_rc.getObjectID().index}; + + // Sums in colinear frame + double pxsum = 0; + double pysum = 0; + double pzsum = 0; + double Esum = 0; + + // Get boost to colinear frame + auto boost = determine_boost(ei, pi); + + auto hfs = hadronicfinalstate->create(0., 0., 0.); + + for (const auto& p: *rcparts) { + + bool isHadron = true; + // Check if it's the scattered electron + if (p.getObjectID().index == ef_rc_id) isHadron = false; + // Check for non-hadron PDG codes + if (p.getPDG() == 11) isHadron = false; + if (p.getPDG() == 22) isHadron = false; + if (p.getPDG() == 13) isHadron = false; + // If it's the scattered electron or not a hadron, skip + if(!isHadron) continue; + + // 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 = apply_boost(boost, hf_lab); + + pxsum += hf_boosted.Px(); + pysum += hf_boosted.Py(); + pzsum += hf_boosted.Pz(); + Esum += hf_boosted.E(); + + hfs.addToHadrons(p); + + } + + // Hadronic final state calculations + auto sigma = Esum - pzsum; + auto pT = sqrt(pxsum*pxsum + pysum*pysum); + auto gamma = (pT*pT - sigma*sigma)/(pT*pT + sigma*sigma); + + hfs.setSigma(sigma); + hfs.setPT(pT); + hfs.setGamma(gamma); + + // Sigma zero or negative + if (sigma <= 0) { + m_log->debug("Sigma zero or negative"); + return; + } + + m_log->debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma()); + + } + +} // namespace Jug::Reco +#endif diff --git a/src/algorithms/reco/HadronicFinalState.h b/src/algorithms/reco/HadronicFinalState.h new file mode 100644 index 0000000000..e4de62ebc9 --- /dev/null +++ b/src/algorithms/reco/HadronicFinalState.h @@ -0,0 +1,48 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Tyler Kutz + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace eicrecon { + + using HadronicFinalStateAlgorithm = algorithms::Algorithm< + algorithms::Input< + edm4hep::MCParticleCollection, + edm4eic::ReconstructedParticleCollection, + edm4eic::MCRecoParticleAssociationCollection + >, + algorithms::Output< + edm4eic::HadronicFinalStateCollection + > + >; + + class HadronicFinalState + : public HadronicFinalStateAlgorithm { + + public: + HadronicFinalState(std::string_view name) + : HadronicFinalStateAlgorithm{name, + {"MCParticles", "inputParticles", "inputAssociations"}, + {"hadronicFinalState"}, + "Calculate summed quantities of the hadronic final state."} {} + + void init(std::shared_ptr& logger); + void process(const Input&, const Output&) const final; + + private: + std::shared_ptr m_log; + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; + }; + +} // namespace eicrecon diff --git a/src/factories/reco/HadronicFinalState_factory.h b/src/factories/reco/HadronicFinalState_factory.h new file mode 100644 index 0000000000..30b0dea850 --- /dev/null +++ b/src/factories/reco/HadronicFinalState_factory.h @@ -0,0 +1,47 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Tyler Kutz + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "extensions/jana/JOmniFactory.h" + +namespace eicrecon { + +template +class HadronicFinalState_factory : + public JOmniFactory> { + +public: + using FactoryT = JOmniFactory>; + +private: + std::unique_ptr m_algo; + + typename FactoryT::template PodioInput m_mc_particles_input {this}; + typename FactoryT::template PodioInput m_rc_particles_input {this}; + typename FactoryT::template PodioInput m_rc_particles_assoc_input {this}; + typename FactoryT::template PodioOutput m_hadronic_final_state_output {this}; + +public: + void Configure() { + m_algo = std::make_unique(this->GetPrefix()); + m_algo->init(this->logger()); + } + + void ChangeRun(int64_t run_number) { + } + + void Process(int64_t run_number, uint64_t event_number) { + m_algo->process({m_mc_particles_input(), m_rc_particles_input(), m_rc_particles_assoc_input()}, + {m_hadronic_final_state_output().get()}); + } +}; + +} // eicrecon diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index b4502b1212..b2b7fae889 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -19,6 +20,9 @@ #include "algorithms/reco/InclusiveKinematicsJB.h" #include "algorithms/reco/InclusiveKinematicsSigma.h" #include "algorithms/reco/InclusiveKinematicseSigma.h" +#if EDM4EIC_VERSION_MAJOR >= 6 +#include "algorithms/reco/HadronicFinalState.h" +#endif #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/meta/CollectionCollector_factory.h" #include "factories/meta/FilterMatching_factory.h" @@ -27,6 +31,9 @@ #include "factories/reco/InclusiveKinematicsTruth_factory.h" #include "factories/reco/JetReconstruction_factory.h" #include "factories/reco/TransformBreitFrame_factory.h" +#if EDM4EIC_VERSION_MAJOR >= 6 +#include "factories/reco/HadronicFinalState_factory.h" +#endif #include "global/reco/ChargedReconstructedParticleSelector_factory.h" #include "global/reco/MC2SmearedParticle_factory.h" #include "global/reco/MatchClusters_factory.h" @@ -277,5 +284,20 @@ void InitPlugin(JApplication *app) { app )); +#if EDM4EIC_VERSION_MAJOR >= 6 + app->Add(new JOmniFactoryGeneratorT>( + "HadronicFinalState", + { + "MCParticles", + "ReconstructedParticles", + "ReconstructedParticleAssociations" + }, + { + "HadronicFinalState" + }, + app + )); +#endif + } } // extern "C" diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index ea562e2ba5..ac94074cbc 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -1,6 +1,8 @@ #include "JEventProcessorPODIO.h" +#include + #include #include #include @@ -186,6 +188,9 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "ReconstructedElectrons", "ScatteredElectronsTruth", "ScatteredElectronsEMinusPz", +#if EDM4EIC_VERSION_MAJOR >= 6 + "HadronicFinalState", +#endif // Track projections "CalorimeterTrackProjections",