diff --git a/src/algorithms/reco/SecondaryVerticesHelix.cc b/src/algorithms/reco/SecondaryVerticesHelix.cc index 3beff6a1cc..748deb95af 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.cc +++ b/src/algorithms/reco/SecondaryVerticesHelix.cc @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -36,8 +37,8 @@ void SecondaryVerticesHelix::init() {} */ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, const SecondaryVerticesHelix::Output& output) const { - const auto [rcvtx, rcparts] = input; - auto [out_secondary_vertices] = output; + const auto [rcvtx, rcparts_collections] = input; + auto [out_secondary_vertices] = output; auto& particleSvc = algorithms::ParticleSvc::instance(); @@ -60,38 +61,44 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, debug("Primary vertex = ({},{},{})cm \t b field = {} tesla", pVtxPos.x, pVtxPos.y, pVtxPos.z, b_field / dd4hep::tesla); - std::vector hVec; - hVec.clear(); - std::vector indexVec; - indexVec.clear(); - for (unsigned int i = 0; const auto& p : *rcparts) { - if (p.getCharge() == 0) - continue; - Helix h(p, b_field); - double dca = h.distance(pVtxPos) * edm4eic::unit::cm; - if (dca < m_cfg.minDca) - continue; - - hVec.push_back(h); - indexVec.push_back(i); - ++i; - } + struct TrackCandidate { + Helix helix; + edm4eic::ReconstructedParticle particle; + float field; + }; - if (hVec.size() != indexVec.size()) - return; + std::vector candidates; + candidates.clear(); + + for (std::size_t i_coll = 0; i_coll < rcparts_collections.size(); ++i_coll) { + const auto& rcparts = rcparts_collections[i_coll]; + const float trackField = (i_coll == 0) ? b_field : 0; + for (const auto& p : *rcparts) { + if (p.getCharge() == 0) { + continue; + } + Helix h(p, trackField); + double dca = h.distance(pVtxPos) * edm4eic::unit::cm; + if (dca < m_cfg.minDca) { + continue; + } + candidates.push_back({h, p, trackField}); + } + } - debug("\tVector size {}, {}", hVec.size(), indexVec.size()); + debug("\tVector size {}", candidates.size()); - for (unsigned int i1 = 0; i1 < hVec.size(); ++i1) { - for (unsigned int i2 = i1 + 1; i2 < hVec.size(); ++i2) { - const auto& p1 = (*rcparts)[indexVec[i1]]; - const auto& p2 = (*rcparts)[indexVec[i2]]; + for (unsigned int i1 = 0; i1 < candidates.size(); ++i1) { + for (unsigned int i2 = i1 + 1; i2 < candidates.size(); ++i2) { + const auto& p1 = candidates[i1].particle; + const auto& p2 = candidates[i2].particle; - if (!(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0)) + if (!(m_cfg.unlikesign && p1.getCharge() + p2.getCharge() == 0)) { continue; + } - const auto& h1 = hVec[i1]; - const auto& h2 = hVec[i2]; + const auto& h1 = candidates[i1].helix; + const auto& h2 = candidates[i2].helix; // Helix function uses cm unit double dca1 = h1.distance(pVtxPos) * edm4eic::unit::cm; @@ -110,8 +117,8 @@ void SecondaryVerticesHelix::process(const SecondaryVerticesHelix::Input& input, continue; edm4hep::Vector3f pairPos = 0.5 * (h1AtDcaTo2 + h2AtDcaTo1); - edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, b_field); - edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, b_field); + edm4hep::Vector3f h1MomAtDca = h1.momentumAt(ss.first, candidates[i1].field); + edm4hep::Vector3f h2MomAtDca = h2.momentumAt(ss.second, candidates[i2].field); edm4hep::Vector3f pairMom = h1MomAtDca + h2MomAtDca; double e1 = diff --git a/src/algorithms/reco/SecondaryVerticesHelix.h b/src/algorithms/reco/SecondaryVerticesHelix.h index 25427d0231..ba5399c27c 100644 --- a/src/algorithms/reco/SecondaryVerticesHelix.h +++ b/src/algorithms/reco/SecondaryVerticesHelix.h @@ -11,15 +11,17 @@ #include #include // for basic_string #include // for string_view +#include #include "algorithms/interfaces/WithPodConfig.h" #include "algorithms/reco/SecondaryVerticesHelixConfig.h" namespace eicrecon { -using SecondaryVerticesHelixAlgorithm = algorithms::Algorithm< - algorithms::Input, - algorithms::Output>; +using SecondaryVerticesHelixAlgorithm = + algorithms::Algorithm>, + algorithms::Output>; class SecondaryVerticesHelix : public SecondaryVerticesHelixAlgorithm, public WithPodConfig { @@ -28,7 +30,7 @@ class SecondaryVerticesHelix : public SecondaryVerticesHelixAlgorithm, SecondaryVerticesHelix(std::string_view name) : SecondaryVerticesHelixAlgorithm{ name, - {"inputVertices", "inputParticles"}, + {"inputVertices", "inputParticleCollections"}, {"outputSecondaryVertices"}, "Reconstruct secondary vertices in SecondaryVertices collection"} {} diff --git a/src/factories/reco/SecondaryVerticesHelix_factory.h b/src/factories/reco/SecondaryVerticesHelix_factory.h index f66259c4f6..f5537d9788 100644 --- a/src/factories/reco/SecondaryVerticesHelix_factory.h +++ b/src/factories/reco/SecondaryVerticesHelix_factory.h @@ -26,7 +26,7 @@ class SecondaryVerticesHelix_factory std::unique_ptr m_algo; PodioInput m_rc_vertices_input{this}; - PodioInput m_rc_parts_input{this}; + VariadicPodioInput m_rc_parts_input{this}; // Declare outputs PodioOutput m_secondary_vertices_output{this}; @@ -50,8 +50,11 @@ class SecondaryVerticesHelix_factory } void Process(int32_t /* run_number */, uint64_t /* event_number */) { - m_algo->process({m_rc_vertices_input(), m_rc_parts_input()}, - {m_secondary_vertices_output().get()}); + std::vector> rc_parts; + for (const auto* coll : m_rc_parts_input()) { + rc_parts.push_back(gsl::not_null{coll}); + } + m_algo->process({m_rc_vertices_input(), rc_parts}, {m_secondary_vertices_output().get()}); } }; diff --git a/src/global/reco/reco.cc b/src/global/reco/reco.cc index a0506d90fc..e60105e117 100644 --- a/src/global/reco/reco.cc +++ b/src/global/reco/reco.cc @@ -314,5 +314,11 @@ void InitPlugin(JApplication* app) { app->Add(new JOmniFactoryGeneratorT( "SecondaryVerticesHelix", {"PrimaryVertices", "ReconstructedParticles"}, {"SecondaryVerticesHelix"}, {}, app)); + + app->Add(new JOmniFactoryGeneratorT( + "SecondaryVerticesHelixInclFarForward", + {"PrimaryVertices", "ReconstructedParticles", "ForwardRomanPotRecParticles", + "ForwardOffMRecParticles"}, + {"SecondaryVerticesHelixInclFarForward"}, {}, app)); } } // extern "C"