From 8bd7f12a0431b2ac5c4b73c676d3e72f3a39c5bb Mon Sep 17 00:00:00 2001 From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch> Date: Fri, 15 Oct 2021 11:57:18 +0200 Subject: [PATCH 01/26] Utilities for the ThOr combiner --- Phys/SelKernel/include/SelKernel/Utilities.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index ff73ce1a562..4d72084e099 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -194,6 +194,13 @@ namespace Sel::Utils { constexpr bool has_tracklike_API = details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::v3::Tracks, T>::value; + template <typename T> + constexpr bool has_neutrals_API = + details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value; + + template <typename T> + constexpr bool has_basic_particle_API = (has_tracklike_API<T> || has_neutrals_API<T>); + using StateLocation = LHCb::Event::v3::Tracks::StateLocation; template <typename Tracklike> auto get_kinematic_state( Tracklike track ) { -- GitLab From fde79c3533d7ab54c244d1f0cae4a3099c484ad0 Mon Sep 17 00:00:00 2001 From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch> Date: Fri, 15 Oct 2021 11:58:53 +0200 Subject: [PATCH 02/26] lb-format --- Phys/SelKernel/include/SelKernel/Utilities.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index 4d72084e099..6a6040c72af 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -196,10 +196,10 @@ namespace Sel::Utils { template <typename T> constexpr bool has_neutrals_API = - details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value; + details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value; template <typename T> - constexpr bool has_basic_particle_API = (has_tracklike_API<T> || has_neutrals_API<T>); + constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> ); using StateLocation = LHCb::Event::v3::Tracks::StateLocation; template <typename Tracklike> -- GitLab From a735b9da9ba9f61745ff1a59e1395b1085c747f1 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 22 Oct 2021 10:11:44 +0100 Subject: [PATCH 03/26] renaming --- Phys/SelKernel/include/SelKernel/Utilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index a105e7dbe37..904378150f6 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -197,7 +197,7 @@ namespace Sel::Utils { template <typename T> constexpr bool has_neutrals_API = - details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::HypothesesWithOrigin, T>::value; + details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithOrigin, T>::value; template <typename T> constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> ); -- GitLab From 448ddfc8aa09139fd7de048c7ba2fd9017de53ed Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 5 Nov 2021 10:48:01 +0000 Subject: [PATCH 04/26] Add empty function to return a state vector from neutral objects --- Phys/SelTools/include/SelTools/Utilities.h | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h index 294caba6621..53fa5b80be1 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -82,4 +82,19 @@ namespace Sel { } return s; } -} // namespace Sel \ No newline at end of file + + /** @fn stateVectorFromOther + @brief Convert a particle (4-momentum + position) to 4-state (at some z). + + TODO: This must be properly filled + */ + template <typename particle_t> + auto stateVectorFromOther( particle_t const& p ) { + + using float_v = std::decay_t<decltype( p.pt() )>; + + State4<float_v> s{}; + + return s; + } +} // namespace Sel -- GitLab From c87bbbd6af3c1107a5010d7b51d09cd970b6819f Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Mon, 22 Nov 2021 12:37:51 +0000 Subject: [PATCH 05/26] Rename state accessor for neutrals --- Phys/SelTools/include/SelTools/Utilities.h | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h index 53fa5b80be1..5f2be8da784 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -83,17 +83,28 @@ namespace Sel { return s; } - /** @fn stateVectorFromOther + /** @fn stateVectorFromNeutral @brief Convert a particle (4-momentum + position) to 4-state (at some z). TODO: This must be properly filled */ template <typename particle_t> - auto stateVectorFromOther( particle_t const& p ) { + auto stateVectorFromNeutral( particle_t const& p ) { using float_v = std::decay_t<decltype( p.pt() )>; + auto const invpz = 1.f / p.pz(); + State4<float_v> s{}; + s.x() = p.hypothesis_position_x(); + s.y() = p.hypothesis_position_y(); + s.z() = p.hypothesis_position_z(); + s.tx() = p.px() * invpz; + s.ty() = p.py() * invpz; + + s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); + s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); + s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); return s; } -- GitLab From 995f7ce19ad67d4f1a62ac4bbbbf38b0143593ad Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Wed, 24 Nov 2021 12:06:51 +0000 Subject: [PATCH 06/26] Fix setting the state as empty --- Phys/SelTools/include/SelTools/Utilities.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h index 5f2be8da784..e06e304c388 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -102,9 +102,9 @@ namespace Sel { s.tx() = p.px() * invpz; s.ty() = p.py() * invpz; - s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); - s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); - s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 3>>(); + s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>(); + s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>(); + s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 2>>(); return s; } -- GitLab From c338b482a11633d007d1dde8309fabf17839338f Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 26 Nov 2021 11:44:34 +0000 Subject: [PATCH 07/26] Change the name of the container --- Phys/SelKernel/include/SelKernel/Utilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index 3c7a6c48ed3..60ba84228ae 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -201,7 +201,7 @@ namespace Sel::Utils { template <typename T> constexpr bool has_neutrals_API = - details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithOrigin, T>::value; + details::is_proxy<T>::value&& details::unpack_and_contains<LHCb::Event::CaloHypothesesWithDirection, T>::value; template <typename T> constexpr bool has_basic_particle_API = ( has_tracklike_API<T> || has_neutrals_API<T> ); -- GitLab From e8603da46a3ddb1cb27fb74159f9954b4f3c22f1 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Thu, 2 Dec 2021 19:35:05 +0000 Subject: [PATCH 08/26] Add some interfaces for neutrals --- Phys/SelAlgorithms/src/MakeZip.cpp | 3 +++ Phys/SelAlgorithms/src/Monitor.cpp | 2 ++ Phys/SelTools/include/SelTools/DistanceCalculator.h | 2 ++ 3 files changed, 7 insertions(+) diff --git a/Phys/SelAlgorithms/src/MakeZip.cpp b/Phys/SelAlgorithms/src/MakeZip.cpp index 0daf3654532..e5dca6a180b 100644 --- a/Phys/SelAlgorithms/src/MakeZip.cpp +++ b/Phys/SelAlgorithms/src/MakeZip.cpp @@ -63,6 +63,9 @@ struct MakeZip final : public Helper<OutputType, Inputs...>::Base { using MakeZip__ChargedBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::ChargedBasics>; DECLARE_COMPONENT_WITH_ID( MakeZip__ChargedBasics_as_Particles, "MakeZip__ChargedBasics_as_Particles" ) +using MakeZip__NeutralBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::NeutralBasics>; +DECLARE_COMPONENT_WITH_ID( MakeZip__NeutralBasics_as_Particles, "MakeZip__NeutralBasics_as_Particles" ) + using MakeZip__PrFittedForwardTracks__BestVertexRelations = MakeZip<deduce_return_type_t, LHCb::Pr::Fitted::Forward::Tracks, BestVertexRelations>; DECLARE_COMPONENT_WITH_ID( MakeZip__PrFittedForwardTracks__BestVertexRelations, diff --git a/Phys/SelAlgorithms/src/Monitor.cpp b/Phys/SelAlgorithms/src/Monitor.cpp index 8154407d841..59230867ba7 100644 --- a/Phys/SelAlgorithms/src/Monitor.cpp +++ b/Phys/SelAlgorithms/src/Monitor.cpp @@ -31,6 +31,8 @@ namespace SelAlgorithms { DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::ChargedBasics>, "Monitor__ChargedBasics" ) + DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::NeutralBasics>, "Monitor__NeutralBasics" ) + DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Event::Composites>, "Monitor__Composites" ) DECLARE_COMPONENT_WITH_ID( Monitor<LHCb::Particle::Range>, "Monitor__ParticleRange" ) diff --git a/Phys/SelTools/include/SelTools/DistanceCalculator.h b/Phys/SelTools/include/SelTools/DistanceCalculator.h index 4328e55e9f3..4b0ef410d79 100644 --- a/Phys/SelTools/include/SelTools/DistanceCalculator.h +++ b/Phys/SelTools/include/SelTools/DistanceCalculator.h @@ -45,6 +45,8 @@ namespace Sel { return p.state( StateLocation::ClosestToBeam ); } else if constexpr ( Sel::type_traits::has_closestToBeamState_v<Particle> ) { return p.closestToBeamState(); + } else if constexpr ( Sel::Utils::has_neutrals_API<Particle> ) { + return stateVectorFromNeutral( p ); } else { // Make a state vector from a particle (position + 4-momentum) return stateVectorFromComposite( p ); -- GitLab From b0103c4af653c1632b02e0e77968cb8a1dd1ecd0 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 3 Dec 2021 10:47:11 +0000 Subject: [PATCH 09/26] Remove the MakeZip__NeutralBasics_as_Particles algorithm --- Phys/SelAlgorithms/src/MakeZip.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/Phys/SelAlgorithms/src/MakeZip.cpp b/Phys/SelAlgorithms/src/MakeZip.cpp index e5dca6a180b..0daf3654532 100644 --- a/Phys/SelAlgorithms/src/MakeZip.cpp +++ b/Phys/SelAlgorithms/src/MakeZip.cpp @@ -63,9 +63,6 @@ struct MakeZip final : public Helper<OutputType, Inputs...>::Base { using MakeZip__ChargedBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::ChargedBasics>; DECLARE_COMPONENT_WITH_ID( MakeZip__ChargedBasics_as_Particles, "MakeZip__ChargedBasics_as_Particles" ) -using MakeZip__NeutralBasics_as_Particles = MakeZip<LHCb::Event::Particles, LHCb::Event::NeutralBasics>; -DECLARE_COMPONENT_WITH_ID( MakeZip__NeutralBasics_as_Particles, "MakeZip__NeutralBasics_as_Particles" ) - using MakeZip__PrFittedForwardTracks__BestVertexRelations = MakeZip<deduce_return_type_t, LHCb::Pr::Fitted::Forward::Tracks, BestVertexRelations>; DECLARE_COMPONENT_WITH_ID( MakeZip__PrFittedForwardTracks__BestVertexRelations, -- GitLab From b15ab2cceefa2bdfab4d8d128f9c7cc5eaa6d0df Mon Sep 17 00:00:00 2001 From: Miguel Ramos Pernas <miguel.ramos.pernas@cern.ch> Date: Fri, 15 Oct 2021 11:56:42 +0200 Subject: [PATCH 10/26] Extend ThOr combiner to use neutrals --- Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h index 88cc1324eac..628568c515d 100644 --- a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h +++ b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h @@ -165,11 +165,11 @@ namespace ThOr::detail::Combiner { template <typename T0, typename T1> auto have_overlap( T0 const& p0, T1 const& p1 ) { - using Sel::Utils::has_tracklike_API; - if constexpr ( Sel::Utils::has_tracklike_API<T0> && Sel::Utils::has_tracklike_API<T1> ) { + using Sel::Utils::has_basic_particle_API; + if constexpr ( Sel::Utils::has_basic_particle_API<T0> && Sel::Utils::has_basic_particle_API<T1> ) { // track-like/track-like return p0.unique_id() == p1.unique_id(); - } else if constexpr ( Sel::Utils::has_tracklike_API<T0> && !Sel::Utils::has_tracklike_API<T1> ) { + } else if constexpr ( Sel::Utils::has_basic_particle_API<T0> && !Sel::Utils::has_basic_particle_API<T1> ) { // track-like/composite using mask_t = decltype( p0.unique_id() == p0.unique_id() ); mask_t mask{false}; @@ -177,7 +177,7 @@ namespace ThOr::detail::Combiner { mask = mask | ( p1.descendant_unique_id( i ) == p0.unique_id() ); } return mask; - } else if constexpr ( !Sel::Utils::has_tracklike_API<T0> && Sel::Utils::has_tracklike_API<T1> ) + } else if constexpr ( !Sel::Utils::has_basic_particle_API<T0> && Sel::Utils::has_basic_particle_API<T1> ) // composite/track-like return have_overlap( p1, p0 ); else { -- GitLab From 8cd40d564cb311afab71c1c39b987e6f0a91d2a6 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 5 Nov 2021 10:49:01 +0000 Subject: [PATCH 11/26] Add support for CALO objects in the vertex fitter --- .../include/VertexFit/ParticleVertexFitter.h | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index a2178deded6..1e2da8abae4 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -118,12 +118,17 @@ namespace Sel::Fitters { enum struct ChildType : unsigned char { TrackWithVelo = 0, TrackWithoutVelo, + Track, Composite, Resonance, Other, NumTypes }; + template <ChildType C> + static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance || + C == ChildType::Other ); + /** Indirectly stolen from TrackVertexUtils::poca */ template <typename float_v, typename mask_v> @@ -180,12 +185,12 @@ namespace Sel::Fitters { } // Class for daughters with a straight-line trajectory - template <typename particle_t, bool represents_composite> + template <typename particle_t, ChildType child_type> struct VertexTraj final { using simd_t = typename particle_t::simd_t; using int_v = typename simd_t::int_v; using float_v = typename simd_t::float_v; - using state_t = std::conditional_t<represents_composite, State4<float_v>, MiniState<float_v>>; + using state_t = std::conditional_t<has_full_state_v<child_type>, State4<float_v>, MiniState<float_v>>; private: particle_t m_particle; // TODO maybe just cache some information from this? @@ -252,10 +257,11 @@ namespace Sel::Fitters { // ParticleVertexFitter this was done with a full specialisation, but // now it would have to be a partial specialisation and it all gets a // bit messy. - if constexpr ( represents_composite ) { + if constexpr ( child_type == ChildType::Composite ) { addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); - } else { + } else if constexpr ( child_type == ChildType::Track ) { addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); + } else { } } @@ -456,8 +462,8 @@ namespace Sel::Fitters { child_t const& child ) { assert( type == ChildType::TrackWithVelo || type == ChildType::TrackWithoutVelo ); using VertexTrack = - VertexTraj<child_t, false /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; + VertexTraj<child_t, ChildType::Track /* flag this represents a track not a weakly-decaying composite */>; + return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )}; } }; @@ -510,7 +516,8 @@ namespace Sel::Fitters { child_t const& child ) { assert( type == ChildType::Composite || type == ChildType::Resonance ); using VertexComposite = - VertexTraj<child_t, true /* flag this represents a weakly-decaying composite and not a track */>; + VertexTraj<child_t, + ChildType::Composite /* flag this represents a weakly-decaying composite and not a track */>; return VertexComposite{child, stateVectorFromComposite( child )}; } }; @@ -533,9 +540,11 @@ namespace Sel::Fitters { template <typename mask_v, typename float_v> static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, - [[maybe_unused]] child_t const& child ) { + child_t const& child ) { assert( type == ChildType::Other ); - return bool{}; + using VertexOther = + VertexTraj<child_t, ChildType::Other /* flag this represents a track not a weakly-decaying composite */>; + return VertexOther{child, stateVectorFromOther( child )}; } }; -- GitLab From 0933d2cbed568b2534f7e29e4445cf0caa285269 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Wed, 24 Nov 2021 11:39:54 +0000 Subject: [PATCH 12/26] Include Neutral as a ChildType in the vertex fitter --- .../include/VertexFit/ParticleVertexFitter.h | 107 +++++++++++++++--- 1 file changed, 92 insertions(+), 15 deletions(-) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index 1e2da8abae4..66719c52a4b 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -118,16 +118,16 @@ namespace Sel::Fitters { enum struct ChildType : unsigned char { TrackWithVelo = 0, TrackWithoutVelo, - Track, Composite, Resonance, + Neutral, Other, NumTypes }; template <ChildType C> static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance || - C == ChildType::Other ); + C == ChildType::Neutral || C == ChildType::Other ); /** Indirectly stolen from TrackVertexUtils::poca */ @@ -259,7 +259,7 @@ namespace Sel::Fitters { // bit messy. if constexpr ( child_type == ChildType::Composite ) { addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); - } else if constexpr ( child_type == ChildType::Track ) { + } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) { addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); } else { } @@ -423,7 +423,7 @@ namespace Sel::Fitters { * @todo Handle TrackWithoutVelo. */ template <typename child_t> - struct TrackChildHelper { + struct TrackWithVeloChildHelper { /** Determine if the given `child_t` type represents a [chunk of] tracks. * @todo Make this more sophisticated and follow changes in the event * model. @@ -460,9 +460,60 @@ namespace Sel::Fitters { static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, child_t const& child ) { - assert( type == ChildType::TrackWithVelo || type == ChildType::TrackWithoutVelo ); + assert( type == ChildType::TrackWithVelo ); using VertexTrack = - VertexTraj<child_t, ChildType::Track /* flag this represents a track not a weakly-decaying composite */>; + VertexTraj<child_t, + ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>; + return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )}; + } + }; + + /** Helper type for identifying and handling TrackWithVelo and + * TrackWithoutVelo children. + * @todo Handle TrackWithoutVelo. + */ + template <typename child_t> + struct TrackWithoutVeloChildHelper { + /** Determine if the given `child_t` type represents a [chunk of] tracks. + * @todo Make this more sophisticated and follow changes in the event + * model. + */ + constexpr static bool value = Sel::Utils::has_tracklike_API<child_t>; + + /** Given that `value` was `true`, derive (at runtime) the `ChildType` of + * `child_t`. This basically means "is it a long or downstream track". + * @todo Implement check and maybe return `ChildType::TrackWithoutVelo`. + * @todo Clearly state what is supported when `child` is a vector type. + * For example, if it is a proxy into a container that supports + * heterogeneous row types (e.g. some rows are long, some rows are + * downstream) then there is not a single correct return value we + * can give here. Is that a runtime error? Should the return type + * be changed so that can be supported? + */ + static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) { + return ChildType::TrackWithoutVelo; + } + + /** Given that `value` was `true`, and that `type` was returned by our + * `deriveType` function, return the narrowest-possible std::variant of + * types implementing the API used in the vertex fit. For tracks there + * should just be one return type (i.e. a variant of size 1), but there + * will be some runtime branching on the `type` to figure out how to get + * the relevant state. + * @todo Use `vertex_position` and its validity mask + * `vertex_position_mask` to determine if a sophisticated + * extrapolation is needed. + * @todo Implement this logic, don't just blindly call the + * `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]]. + */ + template <typename mask_v, typename float_v> + static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, + child_t const& child ) { + assert( type == ChildType::TrackWithoutVelo ); + using VertexTrack = + VertexTraj<child_t, + ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>; return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )}; } }; @@ -522,7 +573,33 @@ namespace Sel::Fitters { } }; - /** Helper type for the default `Other` case. + /** Helper type for the default `Neutral` case. + */ + template <typename child_t> + struct NeutralChildHelper { + /** This helper just tries generic logic, it should be tried last and it + * accepts everything. If it doesn't work, a compile error is the + * correct outcome. + */ + constexpr static bool value = Sel::Utils::has_neutrals_API<child_t>; + /** Catch-all, don't expect other logic here. + */ + static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; } + /** @todo Need to implement the helper class that handles generic + * children and return an instance of it. As a 1-member variant. + */ + template <typename mask_v, typename float_v> + static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, + child_t const& child ) { + assert( type == ChildType::Neutral ); + using VertexNeutral = + VertexTraj<child_t, ChildType::Neutral /* flag this represents a track not a weakly-decaying composite */>; + return VertexNeutral{child, stateVectorFromNeutral( child )}; + } + }; + + /** Helper type for the default `Neutral` case. */ template <typename child_t> struct LeftoverHelper { @@ -533,7 +610,7 @@ namespace Sel::Fitters { constexpr static bool value = true; /** Catch-all, don't expect other logic here. */ - static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Other; } + static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; } /** @todo Need to implement the helper class that handles generic * children and return an instance of it. As a 1-member variant. */ @@ -542,9 +619,7 @@ namespace Sel::Fitters { [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, child_t const& child ) { assert( type == ChildType::Other ); - using VertexOther = - VertexTraj<child_t, ChildType::Other /* flag this represents a track not a weakly-decaying composite */>; - return VertexOther{child, stateVectorFromOther( child )}; + return bool{}; } }; @@ -553,7 +628,8 @@ namespace Sel::Fitters { */ template <typename child_t> using all_helper_types = - boost::mp11::mp_list<TrackChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>; + boost::mp11::mp_list<TrackWithVeloChildHelper<child_t>, TrackWithoutVeloChildHelper<child_t>, + NeutralChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>; /** Helper types that match (::value == true) for a given child type. */ @@ -721,7 +797,7 @@ namespace Sel::Fitters { oss << "Could not initialise Sel::Fitters::ParticleVertex: #TrackWithVelo = " << count( ChildType::TrackWithVelo ) << ", #TrackWithoutVelo = " << count( ChildType::TrackWithoutVelo ) << ", #Composite = " << count( ChildType::Composite ) << ", #Resonance = " << count( ChildType::Resonance ) - << ", #Other = " << count( ChildType::Other ); + << ", #Neutral = " << count( ChildType::Neutral ) << ", #Other = " << count( ChildType::Other ); throw exception( oss.str() ); } } @@ -757,7 +833,8 @@ namespace Sel::Fitters { // - TrackWithoutVelo: tracks without velo hits (only downstream, hopefully) // - Resonance: composites with ctau < 1micron // - Composite: other composites - // - Other: everything else, e.g. photons, pi0, jets + // - Neutral: neutral objects, e.g. photons, pi0 + // - Other: everything else, e.g. jets // Not all of these classifications can be made at compile time based on // the template pack `child_ts...`. If the input children are themselves // variants then `deriveType` is used as a visitor, increasing the amount @@ -780,7 +857,7 @@ namespace Sel::Fitters { // Now estimate everything that goes into the vertex, this array will be // filled with variants; the different variant types encode the different - // handling required for TrackWith[out]Velo/Composite/Resonance/Other. + // handling required for TrackWith[out]Velo/Composite/Resonance/Neutral/Other. std::tuple extrapolated_children{ prepareChild( vertex_position, prelim_fit_mask, types[IChild], children.template get<IChild>() )...}; -- GitLab From f6b9e8a4f4b897971a761d68f767f740ff866976 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 26 Nov 2021 11:46:41 +0000 Subject: [PATCH 13/26] Add a NeutralBasics particle maker --- Phys/ParticleMaker/src/NeutralMakers.cpp | 43 ++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 43d39bdea1a..346f16657a8 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -14,6 +14,7 @@ #include "GaudiAlg/Transformer.h" #include "Event/Particle.h" +#include "Event/Particle_v2.h" #include "Event/ProtoParticle.h" #include "Event/RecVertex.h" @@ -634,4 +635,46 @@ namespace LHCb::Phys::ParticleMakers { debug() << " Skipped : " << m_SkipPi0sCounter << endmsg; debug() << "--------------------" << endmsg; } + + using neutral_basics_maker_output_t = + std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>; + + /** @class NeutralBasicsMaker + * + * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices + * + * + */ + class NeutralBasicsMaker : public Gaudi::Functional::Transformer<neutral_basics_maker_output_t( + LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&, + LHCb::Event::v2::RecVertices const& )> { + + public: + NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, + KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, + KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, + {KeyValue{"OutputParticles", ""}} ) {} + + neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const& unique_id_gen, + LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, + LHCb::Event::v2::RecVertices const& primary_vertices ) const override { + + auto zn = Zipping::generateZipIdentifier(); + auto calo_hypos_with_direction = std::make_unique<LHCb::Event::CaloHypothesesWithDirection>( unique_id_gen, zn ); + + calo_hypos_with_direction->reserve( calo_hypos.size() * primary_vertices.size() ); + + // here we are creating one hypothesis per primary vertex + for ( auto const& ch : calo_hypos ) + for ( auto const& vtx : primary_vertices ) calo_hypos_with_direction->emplace_back( unique_id_gen, ch, vtx ); + + return {LHCb::Pr::make_zip( std::as_const( *calo_hypos_with_direction ) ), + std::move( calo_hypos_with_direction )}; + } + }; + + DECLARE_COMPONENT( NeutralBasicsMaker ) + } // namespace LHCb::Phys::ParticleMakers -- GitLab From a1020cc4760634d9e4709ea8ebb911786090f232 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Mon, 29 Nov 2021 14:03:25 +0000 Subject: [PATCH 14/26] Fix NeutralBasicsMaker --- Phys/ParticleMaker/src/NeutralMakers.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 346f16657a8..e918e242897 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -655,7 +655,7 @@ namespace LHCb::Phys::ParticleMakers { {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, - {KeyValue{"OutputParticles", ""}} ) {} + {KeyValue{"Particles", ""}} ) {} neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const& unique_id_gen, LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, @@ -675,6 +675,6 @@ namespace LHCb::Phys::ParticleMakers { } }; - DECLARE_COMPONENT( NeutralBasicsMaker ) + DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker ) } // namespace LHCb::Phys::ParticleMakers -- GitLab From 53ca06d58c97e7040f04a8afebd42c8555048259 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Wed, 1 Dec 2021 16:31:43 +0000 Subject: [PATCH 15/26] Define a NeutralBasics filter using Pr::Filter --- Phys/ParticleCombiners/src/Filter.cpp | 2 ++ Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h | 7 +++++++ 2 files changed, 9 insertions(+) diff --git a/Phys/ParticleCombiners/src/Filter.cpp b/Phys/ParticleCombiners/src/Filter.cpp index d87b7eee794..db12706628b 100644 --- a/Phys/ParticleCombiners/src/Filter.cpp +++ b/Phys/ParticleCombiners/src/Filter.cpp @@ -28,4 +28,6 @@ namespace LHCb { DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Event::ChargedBasics>, "ChargedBasicsFilter" ) +DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Event::NeutralBasics>, "NeutralBasicsFilter" ) + DECLARE_COMPONENT_WITH_ID( Pr::Filter<LHCb::Particle::Range>, "ParticleRangeFilter" ) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index 66719c52a4b..dedc95c348e 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -261,10 +261,17 @@ namespace Sel::Fitters { addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) { addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); + } else if constexpr ( child_type == ChildType::Neutral ) { + addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix ); } else { } } + /** @todo FIXME correctly process the errors + */ + void addToFourVector_Neutral( VecN<3, float_v> const&, VecN<4, float_v>&, SymNxN<4, float_v>&, + Matrix<4, 3, float_v>& ) const {} + /** @todo FIXME implement logic for bremsstrahlung photons. */ void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, -- GitLab From fefdc0dbd03d69e57b19b5bd61bfdb1d4eac6986 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Thu, 2 Dec 2021 19:34:39 +0000 Subject: [PATCH 16/26] Add a basic composite + neutral particle combiner --- .../src/ThOrCombiner_Composites.cpp | 4 ++++ Phys/ParticleMaker/src/NeutralMakers.cpp | 14 +++++++------- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp index b90e62625d7..57808bdc750 100644 --- a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp +++ b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp @@ -19,6 +19,10 @@ using ThOrCombiner__Composites3ChargedBasics = ThOr::CombinerBest<LHCb::Event::Composites, LHCb::Event::ChargedBasics, LHCb::Event::ChargedBasics, LHCb::Event::ChargedBasics>; +using ThOrCombiner__CompositesNeutralBasics = ThOr::CombinerBest<LHCb::Event::Composites, LHCb::Event::NeutralBasics>; + DECLARE_COMPONENT_WITH_ID( ThOrCombiner__CompositesChargedBasics, "ThOrCombiner__CompositesChargedBasics" ) DECLARE_COMPONENT_WITH_ID( ThOrCombiner__Composites2ChargedBasics, "ThOrCombiner__Composites2ChargedBasics" ) DECLARE_COMPONENT_WITH_ID( ThOrCombiner__Composites3ChargedBasics, "ThOrCombiner__Composites3ChargedBasics" ) + +DECLARE_COMPONENT_WITH_ID( ThOrCombiner__CompositesNeutralBasics, "ThOrCombiner__CompositesNeutralBasics" ) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index e918e242897..3d95931b5f1 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -637,7 +637,7 @@ namespace LHCb::Phys::ParticleMakers { } using neutral_basics_maker_output_t = - std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>; + std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>> ; /** @class NeutralBasicsMaker * @@ -645,17 +645,17 @@ namespace LHCb::Phys::ParticleMakers { * * */ - class NeutralBasicsMaker : public Gaudi::Functional::Transformer<neutral_basics_maker_output_t( + class NeutralBasicsMaker : public Gaudi::Functional::MultiTransformer<neutral_basics_maker_output_t( LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&, LHCb::Event::v2::RecVertices const& )> { public: NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator ) - : Transformer( name, pSvcLocator, - {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, - KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, - KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, - {KeyValue{"Particles", ""}} ) {} + : MultiTransformer( name, pSvcLocator, + {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, + KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, + KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, + {KeyValue{"Particles", ""}, KeyValue{"CaloHypotheses", ""}} ) {} neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const& unique_id_gen, LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, -- GitLab From cef33116009660d349867ee3f7f35b6da678c81e Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 3 Dec 2021 10:30:45 +0000 Subject: [PATCH 17/26] Use CaloHypothesesWithDirection as NeutralBasics --- Phys/ParticleMaker/src/NeutralMakers.cpp | 32 +++++++++++------------- 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 3d95931b5f1..80558d58f7e 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -636,42 +636,38 @@ namespace LHCb::Phys::ParticleMakers { debug() << "--------------------" << endmsg; } - using neutral_basics_maker_output_t = - std::tuple<LHCb::Event::NeutralBasics, std::unique_ptr<LHCb::Event::CaloHypothesesWithDirection>>> ; - /** @class NeutralBasicsMaker * * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices * * */ - class NeutralBasicsMaker : public Gaudi::Functional::MultiTransformer<neutral_basics_maker_output_t( + class NeutralBasicsMaker : public Gaudi::Functional::Transformer<LHCb::Event::NeutralBasics( LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&, LHCb::Event::v2::RecVertices const& )> { public: NeutralBasicsMaker( const std::string& name, ISvcLocator* pSvcLocator ) - : MultiTransformer( name, pSvcLocator, - {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, - KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, - KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, - {KeyValue{"Particles", ""}, KeyValue{"CaloHypotheses", ""}} ) {} + : Transformer( name, pSvcLocator, + {KeyValue{"InputUniqueIDGenerator", LHCb::UniqueIDGeneratorLocation::Default}, + KeyValue{"InputCaloObjects", LHCb::Event::Calo::v2::HypothesesLocation::Default}, + KeyValue{"InputPrimaryVertices", LHCb::Event::v2::RecVertexLocation::Primary}}, + {KeyValue{"Particles", ""}} ) {} - neutral_basics_maker_output_t operator()( LHCb::UniqueIDGenerator const& unique_id_gen, - LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, - LHCb::Event::v2::RecVertices const& primary_vertices ) const override { + LHCb::Event::NeutralBasics operator()( LHCb::UniqueIDGenerator const& unique_id_gen, + LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, + LHCb::Event::v2::RecVertices const& primary_vertices ) const override { - auto zn = Zipping::generateZipIdentifier(); - auto calo_hypos_with_direction = std::make_unique<LHCb::Event::CaloHypothesesWithDirection>( unique_id_gen, zn ); + auto zn = Zipping::generateZipIdentifier(); + LHCb::Event::NeutralBasics neutrals( unique_id_gen, zn ); - calo_hypos_with_direction->reserve( calo_hypos.size() * primary_vertices.size() ); + neutrals.reserve( calo_hypos.size() * primary_vertices.size() ); // here we are creating one hypothesis per primary vertex for ( auto const& ch : calo_hypos ) - for ( auto const& vtx : primary_vertices ) calo_hypos_with_direction->emplace_back( unique_id_gen, ch, vtx ); + for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx ); - return {LHCb::Pr::make_zip( std::as_const( *calo_hypos_with_direction ) ), - std::move( calo_hypos_with_direction )}; + return neutrals; } }; -- GitLab From 837fe1bf77951c9f36b74d9087bd2cac43d96a39 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Mon, 6 Dec 2021 11:18:25 +0000 Subject: [PATCH 18/26] Allow to choose whether to use all the vertices or the average to build the neutrals --- Phys/ParticleMaker/src/NeutralMakers.cpp | 31 +++++++++++++++++++----- 1 file changed, 25 insertions(+), 6 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 80558d58f7e..62243130782 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -17,6 +17,7 @@ #include "Event/Particle_v2.h" #include "Event/ProtoParticle.h" #include "Event/RecVertex.h" +#include "Event/RecVertex_v2.h" #include "CaloFutureUtils/CaloMomentum.h" #include "Kernel/IParticlePropertySvc.h" @@ -658,17 +659,35 @@ namespace LHCb::Phys::ParticleMakers { LHCb::Event::Calo::v2::Hypotheses const& calo_hypos, LHCb::Event::v2::RecVertices const& primary_vertices ) const override { - auto zn = Zipping::generateZipIdentifier(); - LHCb::Event::NeutralBasics neutrals( unique_id_gen, zn ); + LHCb::Event::NeutralBasics neutrals( unique_id_gen ); - neutrals.reserve( calo_hypos.size() * primary_vertices.size() ); + if ( m_AveragePVs.value() ) { + // Calculate the average the primary vertex positions + neutrals.reserve( calo_hypos.size() ); - // here we are creating one hypothesis per primary vertex - for ( auto const& ch : calo_hypos ) - for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx ); + auto avg = std::accumulate( primary_vertices.cbegin(), primary_vertices.cend(), Gaudi::XYZPoint{0.f, 0.f, 0.f}, + [&primary_vertices]( auto const& avg, auto const& vtx ) -> Gaudi::XYZPoint { + return {avg.x() + vtx.position().x() / primary_vertices.size(), + avg.y() + vtx.position().y() / primary_vertices.size(), + avg.z() + vtx.position().z() / primary_vertices.size()}; + } ); + + for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, avg ); + } else { + // Create one particle per primary vertex and ECAL hypothesis + neutrals.reserve( calo_hypos.size() * primary_vertices.size() ); + + for ( auto const& ch : calo_hypos ) + for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx.position() ); + } return neutrals; } + + private: + Gaudi::Property<bool> m_AveragePVs{this, "AveragePVs", true, + "If set to \"true\"; uses the average of the primary vertices in order to " + "determine the pointing direction of the neutral particles"}; }; DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker ) -- GitLab From 3d2312d646c1230e40c1eff670345f0c15445035 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Tue, 7 Dec 2021 16:58:01 +0000 Subject: [PATCH 19/26] Add the four-momenta of the neutral objects corrected by the vertex position to the result of the combination --- .../include/VertexFit/ParticleVertexFitter.h | 24 ++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index dedc95c348e..0502f8e744d 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -267,10 +267,28 @@ namespace Sel::Fitters { } } - /** @todo FIXME correctly process the errors + /** Add the momenta of a neutral particle to the result of a combination + + Update the information of the momenta of the neutral objects based on the vertex of + the decaying particle. It is crucial that we add the neutral objects after + determining the decay vertex (i.e. at the end of the combination). + + @todo FIXME correctly process the covariance matrix */ - void addToFourVector_Neutral( VecN<3, float_v> const&, VecN<4, float_v>&, SymNxN<4, float_v>&, - Matrix<4, 3, float_v>& ) const {} + void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&, + Matrix<4, 3, float_v>& ) const { + // TODO: make this a global function that modifies the neutral objects at the + // persistency stage + auto dx = m_state.x() - vertexpos( 0 ); + auto dy = m_state.y() - vertexpos( 1 ); + auto dz = m_state.z() - vertexpos( 2 ); + + auto d = sqrt( dx * dx + dy * dy + dz * dz ); + + auto p = m_particle.p(); + + p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()}; + } /** @todo FIXME implement logic for bremsstrahlung photons. */ -- GitLab From fd29320e8c877f16e0e055e722b94ff098762a90 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Mon, 17 Jan 2022 12:24:29 +0000 Subject: [PATCH 20/26] Fixes due to the new zipping infrastructure --- Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index 0502f8e744d..6ca5b104046 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -489,7 +489,7 @@ namespace Sel::Fitters { using VertexTrack = VertexTraj<child_t, ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )}; + return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; } }; @@ -539,7 +539,7 @@ namespace Sel::Fitters { using VertexTrack = VertexTraj<child_t, ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( child_t::StateLocation::ClosestToBeam )}; + return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; } }; @@ -788,8 +788,8 @@ namespace Sel::Fitters { velotrajs[nvelotrajs++] = StateVector4<float_v>{child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; } else { - throw exception( - "TrackWithVelo type didn't have a state( child_t::StateLocation::ClosestToBeam)" ); + throw exception( "TrackWithVelo type didn't have a state( " + "LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam)" ); } } else if ( types[i] == ChildType::Composite ) { if constexpr ( Sel::type_traits::has_endVertex_v<child_t> && -- GitLab From e63f366248a8ef18e9121d50e215e7d30dcf29eb Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Tue, 18 Jan 2022 14:53:35 +0000 Subject: [PATCH 21/26] Implement missing member functions for ParticleCombination --- .../SelKernel/include/SelKernel/ParticleCombination.h | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/Phys/SelKernel/include/SelKernel/ParticleCombination.h b/Phys/SelKernel/include/SelKernel/ParticleCombination.h index 3bd96f8bafe..3b8bad9e038 100644 --- a/Phys/SelKernel/include/SelKernel/ParticleCombination.h +++ b/Phys/SelKernel/include/SelKernel/ParticleCombination.h @@ -120,10 +120,19 @@ namespace Sel::detail { return transform_reduce( []( auto const& p ) { return p.charge(); }, std::plus<>{} ); } + auto mass2() const { + auto mom = this->momentum(); + auto px = mom( 0 ); + auto py = mom( 1 ); + auto pz = mom( 2 ); + auto e = mom( 3 ); + return e * e - px * px + py * py + pz * pz; + } + auto mass() const { using Sel::Utils::mass2; using std::sqrt; - return sqrt( mass2( this->fourMomentum() ) ); + return sqrt( mass2( *this ) ); } /** Unfortunately this is needed in order for explicit mask functionality -- GitLab From b203bf2cabee1b977c8ed381f9512b466c5ab71d Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Thu, 31 Mar 2022 10:52:32 +0100 Subject: [PATCH 22/26] Fix compilation issues by defining new members in proxy types --- .../include/CombKernel/ThOrCombiner.h | 8 +- Phys/SelKernel/include/SelKernel/Utilities.h | 6 + Phys/SelTools/include/SelTools/Utilities.h | 2 + .../include/VertexFit/ParticleVertexFitter.h | 1116 ++++++++--------- 4 files changed, 508 insertions(+), 624 deletions(-) diff --git a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h index 17af73e8ca9..252db5d0eff 100644 --- a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h +++ b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h @@ -169,11 +169,10 @@ namespace ThOr::detail::Combiner { template <typename T0, typename T1> auto have_overlap( T0 const& p0, T1 const& p1 ) { - if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<T0> && Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) { + if constexpr ( Sel::Utils::isBasicParticle_v<T0> && Sel::Utils::isBasicParticle_v<T1> ) { // track-like/track-like return p0.unique_id() == p1.unique_id(); - } else if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<T0> && - !Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) { + } else if constexpr ( Sel::Utils::isBasicParticle_v<T0> && !Sel::Utils::isBasicParticle_v<T1> ) { // track-like/composite using mask_t = decltype( p0.unique_id() == p0.unique_id() ); mask_t mask{false}; @@ -181,8 +180,7 @@ namespace ThOr::detail::Combiner { mask = mask | ( p1.descendant_unique_id( i ) == p0.unique_id() ); } return mask; - } else if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<T0> && - Sel::Utils::canBeExtrapolatedDownstream_v<T1> ) { + } else if constexpr ( !Sel::Utils::isBasicParticle_v<T0> && Sel::Utils::isBasicParticle_v<T1> ) { // composite/track-like return have_overlap( p1, p0 ); } else { diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index f1c30e4f6ce..7422842d3a9 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -294,6 +294,12 @@ namespace Sel::Utils { template <typename T> constexpr auto canBeExtrapolatedDownstream_v = std::remove_pointer_t<std::decay_t<T>>::canBeExtrapolatedDownstream; + template <typename T> + constexpr auto isBasicParticle_v = std::remove_pointer_t<std::decay_t<T>>::isBasicParticle; + + template <typename T> + constexpr auto hasTrack_v = std::remove_pointer_t<std::decay_t<T>>::hasTrack; + template <typename T> auto threeMomPosCovMatrix( T const& item ) -> decltype( item.threeMomPosCovMatrix() ) { return item.threeMomPosCovMatrix(); diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h index 6075bac5262..0dd1ebe108a 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -9,10 +9,12 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ #pragma once +#include "Event/Particle_v2.h" #include "LHCbMath/MatVec.h" #include "SelKernel/State4.h" namespace Sel { + /** @fn stateVectorFromComposite * @brief Convert a particle (3-momentum + position) to 4-state (at some z). * diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index b2e88434e79..400448403c8 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -22,6 +22,473 @@ #include "Gaudi/Algorithm.h" +namespace { + + template <std::size_t N, typename T> + using VecN = LHCb::LinAlg::Vec<T, N>; + + template <std::size_t N, typename T> + using SymNxN = LHCb::LinAlg::MatSym<T, N>; + + template <std::size_t N, std::size_t M, typename T> + using Matrix = LHCb::LinAlg::Mat<T, N, M>; + + template <std::size_t N, typename T> + using MatrixNxN = LHCb::LinAlg::Mat<T, N>; + + // Fast utility function for inverting a 2x2 matrix. This just comes + // from ROOT/Math/Dinv, but I took out the test on the determinant + // since we do not need it. + template <typename float_v> + SymNxN<2, float_v> Invert2x2( SymNxN<2, float_v> const& rhs ) { + SymNxN<2, float_v> ret; + auto const det = rhs( 0, 0 ) * rhs( 1, 1 ) - rhs( 0, 1 ) * rhs( 0, 1 ); + auto const invdet = 1.f / det; + // if (det == T(0.)) { return false; } + ret( 0, 0 ) = rhs( 1, 1 ) * invdet; + ret( 1, 1 ) = rhs( 0, 0 ) * invdet; + ret( 0, 1 ) = -1.f * rhs( 0, 1 ) * invdet; + return ret; + } + + // Different daughter types + // TODO look at using the meta enum thing + enum struct ChildType : unsigned char { + TrackWithVelo = 0, + TrackWithoutVelo, + Composite, + Resonance, + Neutral, + Other, + NumTypes + }; + + template <ChildType C> + static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance || + C == ChildType::Neutral || C == ChildType::Other ); + + /* TODO look at merging this with the stuff in State4.h + */ + template <typename F> + struct MiniState { + SymNxN<5, F> cov{}; + auto& x() { return m_x; } + auto const& x() const { return m_x; } + auto& y() { return m_y; } + auto const& y() const { return m_y; } + auto& z() { return m_z; } + auto const& z() const { return m_z; } + auto& qOverP() { return m_qOverP; } + auto const& qOverP() const { return m_qOverP; } + auto& tx() { return m_tx; } + auto const& tx() const { return m_tx; } + auto& ty() { return m_ty; } + auto const& ty() const { return m_ty; } + + template <typename T> + MiniState( T const& other ) + : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} { + static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> ); + auto const& other_cov = other.covariance(); + LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { + LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } ); + } ); + } + + void linearTransportTo( F const& new_z ) { + const auto dz = new_z - m_z; + const auto dz2 = dz * dz; + m_x += dz * m_tx; + m_y += dz * m_ty; + cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 ); + cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) ); + cov( 2, 0 ) += dz * cov( 2, 2 ); + cov( 3, 0 ) += dz * cov( 3, 2 ); + cov( 4, 0 ) += dz * cov( 4, 2 ); + cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 ); + cov( 2, 1 ) += dz * cov( 3, 2 ); + cov( 3, 1 ) += dz * cov( 3, 3 ); + cov( 4, 1 ) += dz * cov( 4, 3 ); + m_z = new_z; + } + + SymNxN<2, F> covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); } + MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); } + + private: + F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{}; + }; + + // Class for daughters with a straight-line trajectory + template <typename particle_t, ChildType child_type> + struct VertexTraj final { + using simd_t = typename particle_t::simd_t; + using int_v = typename simd_t::int_v; + using float_v = typename simd_t::float_v; + using state_t = std::conditional_t<has_full_state_v<child_type>, Sel::State4<float_v>, MiniState<float_v>>; + + private: + particle_t m_particle; // TODO maybe just cache some information from this? + state_t m_state; // TODO evaluate what is really needed + VecN<2, float_v> m_q; // predicted/fitted slope (tx,ty) + SymNxN<2, float_v> m_G; // weight matrix of (x,y) of state + Matrix<2, 3, float_v> m_A; // projection matrix for vertex position + + public: + VertexTraj( particle_t p, state_t state ) + : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {} + + void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX, SymNxN<3, float_v>& halfD2ChisqDX2, + float_v& chi2, int_v& ndof ) { + // move the state (not sure we should do it this way: maybe better to just cache the z.) + m_state.linearTransportTo( vertexpos( 2 ) ); + + // compute the weight matrix + // TODO Just use invChol? + m_G = Invert2x2( m_state.covXX() ); + + // compute residual + VecN<2, float_v> residual{vertexpos( 0 ) - m_state.x(), vertexpos( 1 ) - m_state.y()}; + + // fill the projection matrix: use the fitted momentum! + m_A( 0, 0 ) = m_A( 1, 1 ) = 1.f; + m_A( 0, 2 ) = -1.f * m_q( 0 ); + m_A( 1, 2 ) = -1.f * m_q( 1 ); + + // I tried to make this faster by writing it out, but H does + // not contain sufficiently many zeroes. Better to + // parallelise. + + // I am not sure that the compiler realizes that A^T*G is used + // more than once here, so I'll substitute the code from + // similarity. Note that similarity does not return an expression. + // halfD2ChisqDX2 += ROOT::Math::Similarity( ROOT::Math::Transpose(m_A), m_G ) ; + // halfDChisqDX += ( ROOT::Math::Transpose(m_A) * m_G) * res ; + auto const ATG = m_A.transpose() * m_G.cast_to_mat(); + // updating relevant fit results + halfD2ChisqDX2 = halfD2ChisqDX2 + ( ATG * m_A ).cast_to_sym(); + halfDChisqDX = halfDChisqDX + ATG * residual; + chi2 = chi2 + LHCb::LinAlg::similarity( residual, m_G ); + ndof = ndof + 2; + } + + void updateSlopes( VecN<3, float_v> const& vertexpos ) { + // first update the residual. (note the subtle difference with that in project!) + float_v const dz = vertexpos( 2 ) - m_state.z(); + VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_state.tx() * dz ), + vertexpos( 1 ) - ( m_state.y() + m_state.ty() * dz )}; + + // get the matrix that is the correlation of (x,y) and (tx,ty,qop) + // ROOT::Math::SMatrix<double,2,2> Vba = m_state.covariance().template + // Sub<ROOT::Math::SMatrix<double,2,2>>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit, + // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix<double,2,2> K = Vba*m_G ; compute the momentum vector + m_q = VecN<2, float_v>{{m_state.tx(), m_state.ty()}} + m_state.covXT().transpose() * m_G.cast_to_mat() * residual; + } + + void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, + Matrix<4, 3, float_v>& gainmatrix ) const { + // Despatch is a bit crude, but it'll do...in the original + // ParticleVertexFitter this was done with a full specialisation, but + // now it would have to be a partial specialisation and it all gets a + // bit messy. + if constexpr ( child_type == ChildType::Composite ) { + addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); + } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) { + addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); + } else if constexpr ( child_type == ChildType::Neutral ) { + addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix ); + } else { + } + } + + /** Add the momenta of a neutral particle to the result of a combination + + Update the information of the momenta of the neutral objects based on the vertex of + the decaying particle. It is crucial that we add the neutral objects after + determining the decay vertex (i.e. at the end of the combination). + + @todo FIXME correctly process the covariance matrix + */ + void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&, + Matrix<4, 3, float_v>& ) const { + // TODO: make this a global function that modifies the neutral objects at the + // persistency stage + auto dx = m_state.x() - vertexpos( 0 ); + auto dy = m_state.y() - vertexpos( 1 ); + auto dz = m_state.z() - vertexpos( 2 ); + + auto d = sqrt( dx * dx + dy * dy + dz * dz ); + + auto p = m_particle.p(); + + p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()}; + } + + /** @todo FIXME implement logic for bremsstrahlung photons. + */ + void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, + Matrix<4, 3, float_v>& gain_matrix ) const { + // This specialisation is for TrackWithVelo and TrackWithoutVelo + // we first need to update q/p. since we also need the full gain matrix, we could as well redo that part. + float_v const dz = vertexpos( 2 ) - m_state.z(); + VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_q( 0 ) * dz ), + vertexpos( 1 ) - ( m_state.y() + m_q( 1 ) * dz )}; + + auto const Vba = m_state.cov.template sub<Matrix<3, 2, float_v>, 2, 0>(); + auto K = Vba * m_G.cast_to_mat(); + auto tx = m_state.tx() + ( K * residual )( 0 ); + auto ty = m_state.ty() + ( K * residual )( 1 ); + auto qOverP = m_state.qOverP() + ( K * residual )( 2 ); + // update 4-momentum + // TODO can we just store this value in the constructor? what else is needed? bremsstrahlung info? + auto const mass = m_particle.mass(); + auto const p = 1.f / abs( qOverP ); + auto const n2 = 1.f + tx * tx + ty * ty; + auto const n = sqrt( n2 ); + auto const pz = p / n; + auto const px = tx * pz; + auto const py = ty * pz; + auto const E = sqrt( p * p + mass * mass ); + p4 = p4 + VecN<4, float_v>{{px, py, pz, E}}; // TODO check whether this gives right results + // it was E px py pz before + + Matrix<4, 3, float_v> dp4dmom{}; // jacobi matrix of del(p4)/del(tx,ty,qOverP) + auto const n3 = n2 * n; + auto const invn3 = 1.f / n3; + dp4dmom( 0, 0 ) = p * ( 1.f + ty * ty ) * invn3; // dpx/dtx + dp4dmom( 0, 1 ) = p * tx * -ty * invn3; // dpx/dty + dp4dmom( 0, 2 ) = -px / qOverP; // dpx/dqop + + dp4dmom( 1, 0 ) = p * ty * -tx * invn3; // dpy/dtx + dp4dmom( 1, 1 ) = p * ( 1 + tx * tx ) * invn3; // dpy/dty + dp4dmom( 1, 2 ) = -py / qOverP; // dpy/dqop + + dp4dmom( 2, 0 ) = pz * -tx / n2; // dpz/dtx + dp4dmom( 2, 1 ) = pz * -ty / n2; // dpz/dtx + dp4dmom( 2, 2 ) = -pz / qOverP; // dpz/dqop + + dp4dmom( 3, 0 ) = 0.0; // dE/dtx + dp4dmom( 3, 1 ) = 0.0; // dE/dty + dp4dmom( 3, 2 ) = p / E * -p / qOverP; // dE/dqop + + auto FK = dp4dmom * K; + + p4cov = p4cov + LHCb::LinAlg::similarity( dp4dmom, m_state.cov.template sub<SymNxN<3, float_v>, 2, 2>() ) - + LHCb::LinAlg::similarity( FK, m_state.cov.template sub<SymNxN<2, float_v>, 0, 0>() ); + + gain_matrix = gain_matrix + FK * m_A; + // // For brem-recovered electrons, we need to do something + // // special. So, electrons are an ugly exception here. We could + // // also add to q/p instead, which is cheaper, but perhaps even + // // more ugly. + // if ( particle().particleID().abspid() == 11 ) { + // const double absqop = std::abs( m_state.qOverP() ); + // const double momentumFromParticle = particle().momentum().P(); + // // is 1% a reasonable threshold? + // if ( momentumFromParticle * absqop > 1.01 ) { + // const double momentumFromTrack = 1 / absqop; + // const double momentumError2FromTrack = m_state.covariance()( 4, 4 ) * std::pow( + // momentumFromTrack, 4 ); const double momentumError2FromParticle = + // particle().momCovMatrix()( 3, 3 ) * std::pow( particle().momentum().E() / + // momentumFromParticle, 2 ); + // const double bremEnergyCov = momentumError2FromParticle - momentumError2FromTrack; + // const double bremEnergy = momentumFromParticle - momentumFromTrack; + // // if the correction is unphysical, ignore it. + // if ( bremEnergyCov > 0 ) { + // const auto tx = q( 0 ); + // const auto ty = q( 1 ); + // auto t = std::sqrt( 1 + tx * tx + ty * ty ); + // // we could also 'scale' the momentum, but we anyway need the components for the + // jacobian Gaudi::LorentzVector p4brem{bremEnergy * tx / t, bremEnergy * ty / t, + // bremEnergy / t, bremEnergy}; p4 += p4brem; ROOT::Math::SMatrix<double, 4, 1> J; + // J( 0, 0 ) = tx / t; + // J( 1, 0 ) = ty / t; + // J( 2, 0 ) = 1 / t; + // J( 3, 0 ) = 1; + // p4cov += ROOT::Math::Similarity( J, Gaudi::SymMatrix1x1{bremEnergyCov} ); + // } + // } + // } + } + + void addToFourVector_Composite( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, + Matrix<4, 3, float_v>& gainmatrix ) const { + // first need to 'update' the momentum vector. but for that we + // first need to 'transport'. + using Sel::Utils::fourMomentum; + auto mom = fourMomentum( m_particle ); // 4-momentum + float_v const pz = Z( mom ); + float_v const tx = m_q( 0 ); + float_v const ty = m_q( 1 ); + + // first update the residual + using namespace LHCb::LinAlg; + using Sel::Utils::referencePoint; + auto particlepos = referencePoint( m_particle ); + float_v const dz = Z( vertexpos ) - particlepos.Z(); + VecN<2, float_v> res{X( vertexpos ) - ( particlepos.X() + tx * dz ), + Y( vertexpos ) - ( particlepos.Y() + ty * dz )}; + auto const& R = m_state.covXX(); + auto const& Rinv = m_G; + using Sel::Utils::momCovMatrix; + auto const momCov = momCovMatrix( m_particle ); + using Sel::Utils::momPosCovMatrix; + auto const momPosCov = momPosCovMatrix( m_particle ); + + // To do this right we need THREE projection matrices for the residual: + // r = A*x_vertex + B*mom_dau + C*x_dau + // Note that C=-A. + // Lets call the matrix F^T = (B^T C^T). then the uncertainty on the residual is + // R = F V77 F^T where V77 is the 7x7 covariance matrix of the daughter + // We will need R every iteration + // The correlation matrix between residual and the 7 parameters of the daughter is + // Vr_xmom = V77 F^T + // Consequently, the gain matrix for the momentum-decayvertex is + // K72 = V77 F^T R^-1 + // Now, we do not need the updated parameters for the decay + // vertex. If we just need the momentum part, then this reduces to + // K42 = ( V43 * C^T + V44 * B^T ) R^-1 + + // The easiest implementation is this ... + // Gaudi::Matrix3x2 CT{ ROOT::Math::Transpose(m_A) } ; + // ROOT::Math::SMatrix<double,4,2> BT ; + // BT(0,0) = +dz/pz; + // BT(2,0) = -dz/pz*tx ; + // BT(1,1) = +dz/pz ; + // BT(2,1) = -dz/pz*ty ; + // ROOT::Math::SMatrix<double,4,2> K42 = ( momPosCov * CT + momCov * BT ) * Rinv ; + // but since B and C are largely empty, we better write it out: + float_v const dzopz = dz / pz; + Matrix<4, 2, float_v> K42part; // uninitialised here + for ( int i = 0; i < 4; ++i ) { + K42part( i, 0 ) = + momPosCov( i, 0 ) - momPosCov( i, 2 ) * tx + momCov( i, 0 ) * dzopz - momCov( i, 2 ) * ( dzopz * tx ); + K42part( i, 1 ) = + momPosCov( i, 1 ) - momPosCov( i, 2 ) * ty + momCov( i, 1 ) * dzopz - momCov( i, 2 ) * ( dzopz * ty ); + } + Matrix<4, 2, float_v> const K42 = K42part * Rinv.cast_to_mat(); + + // update the fourvector + auto const deltap4 = K42 * res; + p4 = p4 + mom + deltap4; + // to understand this, compare to what is done for the track + p4cov = p4cov + momCov - similarity( K42, R ); + gainmatrix = gainmatrix + K42 * m_A; + } + }; + + template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition> + struct child_type_for_conditions; + + template <> + struct child_type_for_conditions<true, true, true> { + static constexpr auto value = ChildType::TrackWithVelo; + }; + + template <> + struct child_type_for_conditions<true, false, false> { + static constexpr auto value = ChildType::Neutral; + }; + + template <> + struct child_type_for_conditions<false, false, false> { + static constexpr auto value = ChildType::Composite; + }; + + template <class Proxy> + struct child_type_for_proxy + : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>, + Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {}; + + template <class Proxy> + static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value; + + template <ChildType> + struct ChildHelper; + + /** Helper type for identifying and handling TrackWithVelo and + * TrackWithoutVelo children. + * @todo Handle TrackWithoutVelo. + */ + template <> + struct ChildHelper<ChildType::TrackWithVelo> { + + static constexpr auto child_type = ChildType::TrackWithVelo; + + /** Given that `value` was `true`, and that `type` was returned by our + * `deriveType` function, return the narrowest-possible std::variant of + * types implementing the API used in the vertex fit. For tracks there + * should just be one return type (i.e. a variant of size 1), but there + * will be some runtime branching on the `type` to figure out how to get + * the relevant state. + * @todo Use `vertex_position` and its validity mask + * `vertex_position_mask` to determine if a sophisticated + * extrapolation is needed. + * @todo Implement this logic, don't just blindly call the + * `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]]. + */ + template <typename mask_v, typename float_v, typename proxy> + static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) { + using VertexTrack = + VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>; + return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; + } + }; + + /** Helper type for identifying and handling Composite and Resonance + * children. + * @todo Handle Resonance, adding a runtime check on PID and another type + * to go in the variant returned by prepareChild. + */ + template <> + struct ChildHelper<ChildType::Composite> { + + static constexpr auto child_type = ChildType::Composite; + + /** Given that `value` was `true`, and that `type` was returned by our + * `deriveType` function, return the narrowest-possible std::variant of + * types implementing the API used in the vertex fit. Because weakly + * decaying particles are handled differently in the fit, the return + * type here will have to be a variant of two different types. + * + * The original ParticleVertexFitter didn't use the estimated vertex + * position for `Resonance` and `Composite` children, so we don't + * either. + * @todo Add logic for the `ChildType::Resonance` value of `type` and + * add another member to the return type; remove [[maybe_unused]] + */ + template <typename mask_v, typename float_v, typename proxy> + static auto prepareChild( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) { + using VertexComposite = + VertexTraj<proxy, child_type /* flag this represents a weakly-decaying composite and not a track */>; + return VertexComposite{child, Sel::stateVectorFromComposite( child )}; + } + }; + + /** Helper type for the default `Neutral` case. + */ + template <> + struct ChildHelper<ChildType::Neutral> { + + static constexpr auto child_type = ChildType::Neutral; + + /** @todo Need to implement the helper class that handles generic + * children and return an instance of it. As a 1-member variant. + */ + template <typename mask_v, typename float_v, typename proxy> + static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) { + using VertexNeutral = + VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>; + return VertexNeutral{child, Sel::stateVectorFromNeutral( child )}; + } + }; + + template <class Proxy> + using child_helper_for_proxy_t = ChildHelper<child_type_for_proxy_v<Proxy>>; +} // namespace + namespace Sel::Fitters { /** @class Sel::Fitters::ParticleVertex @@ -44,91 +511,6 @@ namespace Sel::Fitters { } private: - template <std::size_t N, typename T> - using VecN = LHCb::LinAlg::Vec<T, N>; - - template <std::size_t N, typename T> - using SymNxN = LHCb::LinAlg::MatSym<T, N>; - - template <std::size_t N, std::size_t M, typename T> - using Matrix = LHCb::LinAlg::Mat<T, N, M>; - - template <std::size_t N, typename T> - using MatrixNxN = LHCb::LinAlg::Mat<T, N>; - - /* TODO look at merging this with the stuff in State4.h - */ - template <typename F> - struct MiniState { - SymNxN<5, F> cov{}; - auto& x() { return m_x; } - auto const& x() const { return m_x; } - auto& y() { return m_y; } - auto const& y() const { return m_y; } - auto& z() { return m_z; } - auto const& z() const { return m_z; } - auto& qOverP() { return m_qOverP; } - auto const& qOverP() const { return m_qOverP; } - auto& tx() { return m_tx; } - auto const& tx() const { return m_tx; } - auto& ty() { return m_ty; } - auto const& ty() const { return m_ty; } - - template <typename T> - MiniState( T const& other ) - : m_x{other.x()} - , m_y{other.y()} - , m_tx{other.tx()} - , m_ty{other.ty()} - , m_z{other.z()} - , m_qOverP{other.qOverP()} { - static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> ); - auto const& other_cov = other.covariance(); - LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { - LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } ); - } ); - } - - void linearTransportTo( F const& new_z ) { - const auto dz = new_z - m_z; - const auto dz2 = dz * dz; - m_x += dz * m_tx; - m_y += dz * m_ty; - cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 ); - cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) ); - cov( 2, 0 ) += dz * cov( 2, 2 ); - cov( 3, 0 ) += dz * cov( 3, 2 ); - cov( 4, 0 ) += dz * cov( 4, 2 ); - cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 ); - cov( 2, 1 ) += dz * cov( 3, 2 ); - cov( 3, 1 ) += dz * cov( 3, 3 ); - cov( 4, 1 ) += dz * cov( 4, 3 ); - m_z = new_z; - } - - SymNxN<2, F> covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); } - MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); } - - private: - F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{}; - }; - - // Different daughter types - // TODO look at using the meta enum thing - enum struct ChildType : unsigned char { - TrackWithVelo = 0, - TrackWithoutVelo, - Composite, - Resonance, - Neutral, - Other, - NumTypes - }; - - template <ChildType C> - static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance || - C == ChildType::Neutral || C == ChildType::Other ); - /** Indirectly stolen from TrackVertexUtils::poca */ template <typename float_v, typename mask_v> @@ -169,518 +551,21 @@ namespace Sel::Fitters { return {std::string{message}, "Sel::Fitters::ParticleVertex", StatusCode::FAILURE}; } - // Fast utility function for inverting a 2x2 matrix. This just comes - // from ROOT/Math/Dinv, but I took out the test on the determinant - // since we do not need it. - template <typename float_v> - static SymNxN<2, float_v> Invert2x2( SymNxN<2, float_v> const& rhs ) { - SymNxN<2, float_v> ret; - auto const det = rhs( 0, 0 ) * rhs( 1, 1 ) - rhs( 0, 1 ) * rhs( 0, 1 ); - auto const invdet = 1.f / det; - // if (det == T(0.)) { return false; } - ret( 0, 0 ) = rhs( 1, 1 ) * invdet; - ret( 1, 1 ) = rhs( 0, 0 ) * invdet; - ret( 0, 1 ) = -1.f * rhs( 0, 1 ) * invdet; - return ret; - } - - // Class for daughters with a straight-line trajectory - template <typename particle_t, ChildType child_type> - struct VertexTraj final { - using simd_t = typename particle_t::simd_t; - using int_v = typename simd_t::int_v; - using float_v = typename simd_t::float_v; - using state_t = std::conditional_t<has_full_state_v<child_type>, State4<float_v>, MiniState<float_v>>; - - private: - particle_t m_particle; // TODO maybe just cache some information from this? - state_t m_state; // TODO evaluate what is really needed - VecN<2, float_v> m_q; // predicted/fitted slope (tx,ty) - SymNxN<2, float_v> m_G; // weight matrix of (x,y) of state - Matrix<2, 3, float_v> m_A; // projection matrix for vertex position - - public: - VertexTraj( particle_t p, state_t state ) - : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {} - - void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX, - SymNxN<3, float_v>& halfD2ChisqDX2, float_v& chi2, int_v& ndof ) { - // move the state (not sure we should do it this way: maybe better to just cache the z.) - m_state.linearTransportTo( vertexpos( 2 ) ); - - // compute the weight matrix - // TODO Just use invChol? - m_G = Invert2x2( m_state.covXX() ); - - // compute residual - VecN<2, float_v> residual{vertexpos( 0 ) - m_state.x(), vertexpos( 1 ) - m_state.y()}; - - // fill the projection matrix: use the fitted momentum! - m_A( 0, 0 ) = m_A( 1, 1 ) = 1.f; - m_A( 0, 2 ) = -1.f * m_q( 0 ); - m_A( 1, 2 ) = -1.f * m_q( 1 ); - - // I tried to make this faster by writing it out, but H does - // not contain sufficiently many zeroes. Better to - // parallelise. - - // I am not sure that the compiler realizes that A^T*G is used - // more than once here, so I'll substitute the code from - // similarity. Note that similarity does not return an expression. - // halfD2ChisqDX2 += ROOT::Math::Similarity( ROOT::Math::Transpose(m_A), m_G ) ; - // halfDChisqDX += ( ROOT::Math::Transpose(m_A) * m_G) * res ; - auto const ATG = m_A.transpose() * m_G.cast_to_mat(); - // updating relevant fit results - halfD2ChisqDX2 = halfD2ChisqDX2 + ( ATG * m_A ).cast_to_sym(); - halfDChisqDX = halfDChisqDX + ATG * residual; - chi2 = chi2 + LHCb::LinAlg::similarity( residual, m_G ); - ndof = ndof + 2; - } - - void updateSlopes( VecN<3, float_v> const& vertexpos ) { - // first update the residual. (note the subtle difference with that in project!) - float_v const dz = vertexpos( 2 ) - m_state.z(); - VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_state.tx() * dz ), - vertexpos( 1 ) - ( m_state.y() + m_state.ty() * dz )}; - - // get the matrix that is the correlation of (x,y) and (tx,ty,qop) - // ROOT::Math::SMatrix<double,2,2> Vba = m_state.covariance().template - // Sub<ROOT::Math::SMatrix<double,2,2>>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit, - // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix<double,2,2> K = Vba*m_G ; compute the momentum vector - m_q = - VecN<2, float_v>{{m_state.tx(), m_state.ty()}} + m_state.covXT().transpose() * m_G.cast_to_mat() * residual; - } - - void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, - Matrix<4, 3, float_v>& gainmatrix ) const { - // Despatch is a bit crude, but it'll do...in the original - // ParticleVertexFitter this was done with a full specialisation, but - // now it would have to be a partial specialisation and it all gets a - // bit messy. - if constexpr ( child_type == ChildType::Composite ) { - addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); - } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) { - addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); - } else if constexpr ( child_type == ChildType::Neutral ) { - addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix ); - } else { - } - } - - /** Add the momenta of a neutral particle to the result of a combination - - Update the information of the momenta of the neutral objects based on the vertex of - the decaying particle. It is crucial that we add the neutral objects after - determining the decay vertex (i.e. at the end of the combination). - - @todo FIXME correctly process the covariance matrix - */ - void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&, - Matrix<4, 3, float_v>& ) const { - // TODO: make this a global function that modifies the neutral objects at the - // persistency stage - auto dx = m_state.x() - vertexpos( 0 ); - auto dy = m_state.y() - vertexpos( 1 ); - auto dz = m_state.z() - vertexpos( 2 ); - - auto d = sqrt( dx * dx + dy * dy + dz * dz ); - - auto p = m_particle.p(); - - p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()}; - } - - /** @todo FIXME implement logic for bremsstrahlung photons. - */ - void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, - Matrix<4, 3, float_v>& gain_matrix ) const { - // This specialisation is for TrackWithVelo and TrackWithoutVelo - // we first need to update q/p. since we also need the full gain matrix, we could as well redo that part. - float_v const dz = vertexpos( 2 ) - m_state.z(); - VecN<2, float_v> const residual{vertexpos( 0 ) - ( m_state.x() + m_q( 0 ) * dz ), - vertexpos( 1 ) - ( m_state.y() + m_q( 1 ) * dz )}; - - auto const Vba = m_state.cov.template sub<Matrix<3, 2, float_v>, 2, 0>(); - auto K = Vba * m_G.cast_to_mat(); - auto tx = m_state.tx() + ( K * residual )( 0 ); - auto ty = m_state.ty() + ( K * residual )( 1 ); - auto qOverP = m_state.qOverP() + ( K * residual )( 2 ); - // update 4-momentum - // TODO can we just store this value in the constructor? what else is needed? bremsstrahlung info? - auto const mass = m_particle.mass(); - auto const p = 1.f / abs( qOverP ); - auto const n2 = 1.f + tx * tx + ty * ty; - auto const n = sqrt( n2 ); - auto const pz = p / n; - auto const px = tx * pz; - auto const py = ty * pz; - auto const E = sqrt( p * p + mass * mass ); - p4 = p4 + VecN<4, float_v>{{px, py, pz, E}}; // TODO check whether this gives right results - // it was E px py pz before - - Matrix<4, 3, float_v> dp4dmom{}; // jacobi matrix of del(p4)/del(tx,ty,qOverP) - auto const n3 = n2 * n; - auto const invn3 = 1.f / n3; - dp4dmom( 0, 0 ) = p * ( 1.f + ty * ty ) * invn3; // dpx/dtx - dp4dmom( 0, 1 ) = p * tx * -ty * invn3; // dpx/dty - dp4dmom( 0, 2 ) = -px / qOverP; // dpx/dqop - - dp4dmom( 1, 0 ) = p * ty * -tx * invn3; // dpy/dtx - dp4dmom( 1, 1 ) = p * ( 1 + tx * tx ) * invn3; // dpy/dty - dp4dmom( 1, 2 ) = -py / qOverP; // dpy/dqop - - dp4dmom( 2, 0 ) = pz * -tx / n2; // dpz/dtx - dp4dmom( 2, 1 ) = pz * -ty / n2; // dpz/dtx - dp4dmom( 2, 2 ) = -pz / qOverP; // dpz/dqop - - dp4dmom( 3, 0 ) = 0.0; // dE/dtx - dp4dmom( 3, 1 ) = 0.0; // dE/dty - dp4dmom( 3, 2 ) = p / E * -p / qOverP; // dE/dqop - - auto FK = dp4dmom * K; - - p4cov = p4cov + LHCb::LinAlg::similarity( dp4dmom, m_state.cov.template sub<SymNxN<3, float_v>, 2, 2>() ) - - LHCb::LinAlg::similarity( FK, m_state.cov.template sub<SymNxN<2, float_v>, 0, 0>() ); - - gain_matrix = gain_matrix + FK * m_A; - // // For brem-recovered electrons, we need to do something - // // special. So, electrons are an ugly exception here. We could - // // also add to q/p instead, which is cheaper, but perhaps even - // // more ugly. - // if ( particle().particleID().abspid() == 11 ) { - // const double absqop = std::abs( m_state.qOverP() ); - // const double momentumFromParticle = particle().momentum().P(); - // // is 1% a reasonable threshold? - // if ( momentumFromParticle * absqop > 1.01 ) { - // const double momentumFromTrack = 1 / absqop; - // const double momentumError2FromTrack = m_state.covariance()( 4, 4 ) * std::pow( - // momentumFromTrack, 4 ); const double momentumError2FromParticle = - // particle().momCovMatrix()( 3, 3 ) * std::pow( particle().momentum().E() / - // momentumFromParticle, 2 ); - // const double bremEnergyCov = momentumError2FromParticle - momentumError2FromTrack; - // const double bremEnergy = momentumFromParticle - momentumFromTrack; - // // if the correction is unphysical, ignore it. - // if ( bremEnergyCov > 0 ) { - // const auto tx = q( 0 ); - // const auto ty = q( 1 ); - // auto t = std::sqrt( 1 + tx * tx + ty * ty ); - // // we could also 'scale' the momentum, but we anyway need the components for the - // jacobian Gaudi::LorentzVector p4brem{bremEnergy * tx / t, bremEnergy * ty / t, - // bremEnergy / t, bremEnergy}; p4 += p4brem; ROOT::Math::SMatrix<double, 4, 1> J; - // J( 0, 0 ) = tx / t; - // J( 1, 0 ) = ty / t; - // J( 2, 0 ) = 1 / t; - // J( 3, 0 ) = 1; - // p4cov += ROOT::Math::Similarity( J, Gaudi::SymMatrix1x1{bremEnergyCov} ); - // } - // } - // } - } - - void addToFourVector_Composite( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, - SymNxN<4, float_v>& p4cov, Matrix<4, 3, float_v>& gainmatrix ) const { - // first need to 'update' the momentum vector. but for that we - // first need to 'transport'. - using Sel::Utils::fourMomentum; - auto mom = fourMomentum( m_particle ); // 4-momentum - float_v const pz = Z( mom ); - float_v const tx = m_q( 0 ); - float_v const ty = m_q( 1 ); - - // first update the residual - using namespace LHCb::LinAlg; - using Sel::Utils::referencePoint; - auto particlepos = referencePoint( m_particle ); - float_v const dz = Z( vertexpos ) - particlepos.Z(); - VecN<2, float_v> res{X( vertexpos ) - ( particlepos.X() + tx * dz ), - Y( vertexpos ) - ( particlepos.Y() + ty * dz )}; - auto const& R = m_state.covXX(); - auto const& Rinv = m_G; - using Sel::Utils::momCovMatrix; - auto const momCov = momCovMatrix( m_particle ); - using Sel::Utils::momPosCovMatrix; - auto const momPosCov = momPosCovMatrix( m_particle ); - - // To do this right we need THREE projection matrices for the residual: - // r = A*x_vertex + B*mom_dau + C*x_dau - // Note that C=-A. - // Lets call the matrix F^T = (B^T C^T). then the uncertainty on the residual is - // R = F V77 F^T where V77 is the 7x7 covariance matrix of the daughter - // We will need R every iteration - // The correlation matrix between residual and the 7 parameters of the daughter is - // Vr_xmom = V77 F^T - // Consequently, the gain matrix for the momentum-decayvertex is - // K72 = V77 F^T R^-1 - // Now, we do not need the updated parameters for the decay - // vertex. If we just need the momentum part, then this reduces to - // K42 = ( V43 * C^T + V44 * B^T ) R^-1 - - // The easiest implementation is this ... - // Gaudi::Matrix3x2 CT{ ROOT::Math::Transpose(m_A) } ; - // ROOT::Math::SMatrix<double,4,2> BT ; - // BT(0,0) = +dz/pz; - // BT(2,0) = -dz/pz*tx ; - // BT(1,1) = +dz/pz ; - // BT(2,1) = -dz/pz*ty ; - // ROOT::Math::SMatrix<double,4,2> K42 = ( momPosCov * CT + momCov * BT ) * Rinv ; - // but since B and C are largely empty, we better write it out: - float_v const dzopz = dz / pz; - Matrix<4, 2, float_v> K42part; // uninitialised here - for ( int i = 0; i < 4; ++i ) { - K42part( i, 0 ) = - momPosCov( i, 0 ) - momPosCov( i, 2 ) * tx + momCov( i, 0 ) * dzopz - momCov( i, 2 ) * ( dzopz * tx ); - K42part( i, 1 ) = - momPosCov( i, 1 ) - momPosCov( i, 2 ) * ty + momCov( i, 1 ) * dzopz - momCov( i, 2 ) * ( dzopz * ty ); - } - Matrix<4, 2, float_v> const K42 = K42part * Rinv.cast_to_mat(); - - // update the fourvector - auto const deltap4 = K42 * res; - p4 = p4 + mom + deltap4; - // to understand this, compare to what is done for the track - p4cov = p4cov + momCov - similarity( K42, R ); - gainmatrix = gainmatrix + K42 * m_A; - } - }; - - /** Helper type for identifying and handling TrackWithVelo and - * TrackWithoutVelo children. - * @todo Handle TrackWithoutVelo. - */ - template <typename child_t> - struct TrackWithVeloChildHelper { - /** Determine if the given `child_t` type represents a [chunk of] tracks. - * @todo Make this more sophisticated and follow changes in the event - * model. - */ - constexpr static bool value = Sel::Utils::canBeExtrapolatedDownstream_v<child_t>; - - /** Given that `value` was `true`, derive (at runtime) the `ChildType` of - * `child_t`. This basically means "is it a long or downstream track". - * @todo Implement check and maybe return `ChildType::TrackWithoutVelo`. - * @todo Clearly state what is supported when `child` is a vector type. - * For example, if it is a proxy into a container that supports - * heterogeneous row types (e.g. some rows are long, some rows are - * downstream) then there is not a single correct return value we - * can give here. Is that a runtime error? Should the return type - * be changed so that can be supported? - */ - static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) { - return ChildType::TrackWithVelo; - } - - /** Given that `value` was `true`, and that `type` was returned by our - * `deriveType` function, return the narrowest-possible std::variant of - * types implementing the API used in the vertex fit. For tracks there - * should just be one return type (i.e. a variant of size 1), but there - * will be some runtime branching on the `type` to figure out how to get - * the relevant state. - * @todo Use `vertex_position` and its validity mask - * `vertex_position_mask` to determine if a sophisticated - * extrapolation is needed. - * @todo Implement this logic, don't just blindly call the - * `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]]. - */ - template <typename mask_v, typename float_v> - static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, - child_t const& child ) { - assert( type == ChildType::TrackWithVelo ); - using VertexTrack = - VertexTraj<child_t, - ChildType::TrackWithVelo /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; - } - }; - - /** Helper type for identifying and handling TrackWithVelo and - * TrackWithoutVelo children. - * @todo Handle TrackWithoutVelo. - */ - template <typename child_t> - struct TrackWithoutVeloChildHelper { - /** Determine if the given `child_t` type represents a [chunk of] tracks. - * @todo Make this more sophisticated and follow changes in the event - * model. - */ - constexpr static bool value = Sel::Utils::has_tracklike_API<child_t>; - - /** Given that `value` was `true`, derive (at runtime) the `ChildType` of - * `child_t`. This basically means "is it a long or downstream track". - * @todo Implement check and maybe return `ChildType::TrackWithoutVelo`. - * @todo Clearly state what is supported when `child` is a vector type. - * For example, if it is a proxy into a container that supports - * heterogeneous row types (e.g. some rows are long, some rows are - * downstream) then there is not a single correct return value we - * can give here. Is that a runtime error? Should the return type - * be changed so that can be supported? - */ - static ChildType deriveType( LHCb::IParticlePropertySvc const&, [[maybe_unused]] child_t const& child ) { - return ChildType::TrackWithoutVelo; - } - - /** Given that `value` was `true`, and that `type` was returned by our - * `deriveType` function, return the narrowest-possible std::variant of - * types implementing the API used in the vertex fit. For tracks there - * should just be one return type (i.e. a variant of size 1), but there - * will be some runtime branching on the `type` to figure out how to get - * the relevant state. - * @todo Use `vertex_position` and its validity mask - * `vertex_position_mask` to determine if a sophisticated - * extrapolation is needed. - * @todo Implement this logic, don't just blindly call the - * `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]]. - */ - template <typename mask_v, typename float_v> - static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, - child_t const& child ) { - assert( type == ChildType::TrackWithoutVelo ); - using VertexTrack = - VertexTraj<child_t, - ChildType::TrackWithoutVelo /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; - } - }; - - /** Helper type for identifying and handling Composite and Resonance - * children. - * @todo Handle Resonance, adding a runtime check on PID and another type - * to go in the variant returned by prepareChild. - */ - template <typename child_t> - struct CompositeChildHelper { - /** Determine if the given `child_t` represents a [chunk of] composite - * particles. - */ - constexpr static bool value = !Sel::Utils::canBeExtrapolatedDownstream_v<child_t>; - - /** Given that `value` was `true`, derive (at runtime) the `ChildType` of - * `child_t`. This basically means "is it a weakly decaying particle". - * @todo Implement check and maybe return `ChildType::Resonance`. - * @todo Can use `pp_svc` to get the "ctau" of the given `child`. Remove - * [[maybe_unused]] once we do that. - */ - static ChildType deriveType( [[maybe_unused]] LHCb::IParticlePropertySvc const& pp_svc, - [[maybe_unused]] child_t const& child ) { - return ChildType::Composite; - } - /** Given that `value` was `true`, and that `type` was returned by our - * `deriveType` function, return the narrowest-possible std::variant of - * types implementing the API used in the vertex fit. Because weakly - * decaying particles are handled differently in the fit, the return - * type here will have to be a variant of two different types. - * - * The original ParticleVertexFitter didn't use the estimated vertex - * position for `Resonance` and `Composite` children, so we don't - * either. - * @todo Add logic for the `ChildType::Resonance` value of `type` and - * add another member to the return type; remove [[maybe_unused]] - */ - template <typename mask_v, typename float_v> - static auto prepareChild( VecN<3, float_v> const&, mask_v const&, [[maybe_unused]] ChildType type, - child_t const& child ) { - assert( type == ChildType::Composite || type == ChildType::Resonance ); - using VertexComposite = - VertexTraj<child_t, - ChildType::Composite /* flag this represents a weakly-decaying composite and not a track */>; - return VertexComposite{child, stateVectorFromComposite( child )}; - } - }; - - /** Helper type for the default `Neutral` case. - */ - template <typename child_t> - struct NeutralChildHelper { - /** This helper just tries generic logic, it should be tried last and it - * accepts everything. If it doesn't work, a compile error is the - * correct outcome. - */ - constexpr static bool value = Sel::Utils::has_neutrals_API<child_t>; - /** Catch-all, don't expect other logic here. - */ - static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; } - /** @todo Need to implement the helper class that handles generic - * children and return an instance of it. As a 1-member variant. - */ - template <typename mask_v, typename float_v> - static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, - child_t const& child ) { - assert( type == ChildType::Neutral ); - using VertexNeutral = - VertexTraj<child_t, ChildType::Neutral /* flag this represents a track not a weakly-decaying composite */>; - return VertexNeutral{child, stateVectorFromNeutral( child )}; - } - }; - - /** Helper type for the default `Neutral` case. - */ - template <typename child_t> - struct LeftoverHelper { - /** This helper just tries generic logic, it should be tried last and it - * accepts everything. If it doesn't work, a compile error is the - * correct outcome. - */ - constexpr static bool value = true; - /** Catch-all, don't expect other logic here. - */ - static ChildType deriveType( LHCb::IParticlePropertySvc const&, child_t const& ) { return ChildType::Neutral; } - /** @todo Need to implement the helper class that handles generic - * children and return an instance of it. As a 1-member variant. - */ - template <typename mask_v, typename float_v> - static auto prepareChild( [[maybe_unused]] VecN<3, float_v> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, [[maybe_unused]] ChildType type, - child_t const& child ) { - assert( type == ChildType::Other ); - return bool{}; - } - }; - - /** Master list of helper types. The order matters; the **first** matching - * helper is preferred. - */ - template <typename child_t> - using all_helper_types = - boost::mp11::mp_list<TrackWithVeloChildHelper<child_t>, TrackWithoutVeloChildHelper<child_t>, - NeutralChildHelper<child_t>, CompositeChildHelper<child_t>, LeftoverHelper<child_t>>; - - /** Helper types that match (::value == true) for a given child type. - */ - template <typename child_t> - using matching_helper_types = boost::mp11::mp_copy_if<all_helper_types<child_t>, boost::mp11::mp_to_bool>; - - /** The first matching helper type for a given child type. - */ - template <typename child_t> - using first_matching_helper_type = boost::mp11::mp_front<matching_helper_types<child_t>>; - - /** Take a [SIMD vector worth of] particles and assign a single "type" - * value that will steer the fit algorithm below. It is an error - * (unfortunately only a runtime one...) to give children that are of - * heterogeneous `ChildType`, but it should be supported that there is - * a mix of actual PDG IDs (e.g. `child` could be a mix of D0 and D~0, but - * there should be a single return type (ChildType::Composite, in that - * case). - */ - template <typename child_t> - ChildType deriveType( child_t const& child ) const { - assert( m_pp_svc ); - using helper_t = first_matching_helper_type<child_t>; - return helper_t::deriveType( *m_pp_svc, child ); - } - /** Helper for getting the return type of the prepareChild method of a * helper type. */ - template <typename mask_v, typename float_v> - struct prepareChild_helper { - template <typename helper_t, typename variant_member_t> - using fn = decltype( helper_t::prepareChild( std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), - std::declval<ChildType>(), std::declval<variant_member_t>() ) ); + template <typename mask_v, typename float_v, typename VariantTypes> + struct return_types_helper; + + template <typename mask_v, typename float_v, template <class...> class Variant, typename... Proxy> + struct return_types_helper<mask_v, float_v, Variant<Proxy...>> { + using type = std::variant<decltype( child_helper_for_proxy_t<Proxy>::prepareChild( + std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>; }; + template <typename mask_v, typename float_v, typename VariantTypes> + using return_types_helper_t = typename return_types_helper<mask_v, float_v, VariantTypes>::type; + /** Choose the first match for the given child type in `all_helper_types` * and call its `prepareChild` method. All the extra complexity is to cope * with `child_t` being a variant type whose members might match different @@ -689,28 +574,21 @@ namespace Sel::Fitters { * @todo Include the `vertexinsidefoil` logic from ParticleVertexFitter.cpp */ template <typename mask_v, typename float_v, typename child_t> - auto prepareChild( VecN<3, float_v> const& vertex_position, mask_v const& vertex_position_mask, ChildType type, + auto prepareChild( VecN<3, float_v> const& vertex_position, mask_v const& vertex_position_mask, child_t const& possibly_variant_child ) const { // child_t might be variant<Ts...>, or it might be a value type itself. // if the former, this yields tuple<Ts...>, otherwise tuple<child_t> using variant_types = typename LHCb::is_variant<child_t>::tuple_type; - // get the [first] matching helper for each element in `variant_types` - using matched_helpers = boost::mp11::mp_transform<first_matching_helper_type, variant_types>; - // get the return types of the `prepareChild` methods of all of these - // helpers when they are invoked with the corresponding element of - // `variant_types`. e.g. this might be tuple<variant<A, B>, variant<C>> - using return_types = - boost::mp11::mp_transform_q<prepareChild_helper<mask_v, float_v>, matched_helpers, variant_types>; - // Turn this into a variant<...> type that can accommodate all relevant - // types, in the example just above `ret_t` would be variant<A, B, C> - using ret_t = boost::mp11::mp_rename<return_types, std::variant>; + + using ret_t = return_types_helper_t<mask_v, float_v, variant_types>; + static_assert( boost::mp11::mp_is_set<ret_t>::value, "Got duplicate entries in derived variant type. Looks like the code needs an mp_unique<...>." ); // Now we know what return type to impose, invoke/visit on `child` return LHCb::invoke_or_visit( [&]( auto const& child ) -> ret_t { - using helper_t = first_matching_helper_type<std::decay_t<decltype( child )>>; - return helper_t::prepareChild( vertex_position, vertex_position_mask, type, child ); + using helper_t = child_helper_for_proxy_t<std::decay_t<decltype( child )>>; + return helper_t::prepareChild( vertex_position, vertex_position_mask, child ); }, possibly_variant_child ); } @@ -748,7 +626,7 @@ namespace Sel::Fitters { LHCb::invoke_or_visit( [&]( auto const& child ) { using Sel::Utils::endVertexPos; - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) { + if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) { tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} ); } else { throw exception( "Resonance type didn't have an endVertex()" ); @@ -770,18 +648,17 @@ namespace Sel::Fitters { // a variant type. LHCb::invoke_or_visit( [&]( auto const& child ) { + using child_type = decltype( child ); + using Sel::Utils::trackState; using Sel::Utils::referencePoint; using Sel::Utils::threeMomentum; - if ( types[i] == ChildType::TrackWithVelo ) { - if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) { + + if constexpr ( Sel::Utils::isBasicParticle_v<child_type> ) { + if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v<child_type> ) velotrajs[nvelotrajs++] = {trackState( child )}; - } else { - throw exception( "TrackWithVelo type didn't have a state( " - "LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam)" ); - } - } else if ( types[i] == ChildType::Composite ) { - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) { + } else { + if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<child_type> ) { velotrajs[nvelotrajs++] = {referencePoint( child ), threeMomentum( child )}; } else { throw exception( "Composite type didn't one/both of referencePoint() and threeMomentum()" ); @@ -852,8 +729,9 @@ namespace Sel::Fitters { // the template pack `child_ts...`. If the input children are themselves // variants then `deriveType` is used as a visitor, increasing the amount // of runtime dispatch. - std::array const types{LHCb::invoke_or_visit( [this]( auto const& child ) { return this->deriveType( child ); }, - children.template get<IChild>() )...}; + std::array const types{LHCb::invoke_or_visit( + [this]( auto const& child ) { return child_type_for_proxy_v<std::decay_t<decltype( child )>>; }, + children.template get<IChild>() )...}; // Obtain an estimate of the vertex position for initialization auto [vertex_position, vertex_position_mask] = calculateInitialPosition<simd_t>( types, children ); @@ -872,7 +750,7 @@ namespace Sel::Fitters { // filled with variants; the different variant types encode the different // handling required for TrackWith[out]Velo/Composite/Resonance/Neutral/Other. std::tuple extrapolated_children{ - prepareChild( vertex_position, prelim_fit_mask, types[IChild], children.template get<IChild>() )...}; + prepareChild( vertex_position, prelim_fit_mask, children.template get<IChild>() )...}; // These are required for populating child relations at the end, but the // API for dealing with them should be agnostic to the child types @@ -883,7 +761,7 @@ namespace Sel::Fitters { auto num_descendants = ( LHCb::invoke_or_visit( []( [[maybe_unused]] auto const& child ) -> std::size_t { - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) + if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) return child.numDescendants(); else return 1; @@ -896,7 +774,7 @@ namespace Sel::Fitters { ( LHCb::invoke_or_visit( [¤t_unique_ids]( auto const& child ) -> void { - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v<decltype( child )> ) { + if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) { for ( auto i = 0u; i < child.numDescendants(); ++i ) current_unique_ids.emplace_back( child.descendant_unique_id( i ) ); } else { -- GitLab From 8e01508e4c0087477d9fac2efc889f2df53a952e Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Thu, 31 Mar 2022 16:35:25 +0100 Subject: [PATCH 23/26] Redefine states to allow neutrals --- Phys/ParticleMaker/src/NeutralMakers.cpp | 1 - Phys/SelTools/include/SelTools/Utilities.h | 109 ++++++- .../include/VertexFit/ParticleVertexFitter.h | 274 ++++++++---------- 3 files changed, 221 insertions(+), 163 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 62243130782..1fadde1eb4f 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -641,7 +641,6 @@ namespace LHCb::Phys::ParticleMakers { * * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices * - * */ class NeutralBasicsMaker : public Gaudi::Functional::Transformer<LHCb::Event::NeutralBasics( LHCb::UniqueIDGenerator const&, LHCb::Event::Calo::v2::Hypotheses const&, diff --git a/Phys/SelTools/include/SelTools/Utilities.h b/Phys/SelTools/include/SelTools/Utilities.h index 0dd1ebe108a..f2f9f847707 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -15,6 +15,97 @@ namespace Sel { + /** Common state type for non-composite objects with covariance matrix information + + TODO look at merging this with the stuff in State4.h + */ + template <typename FloatType> + struct MiniState { + LHCb::LinAlg::MatSym<FloatType, 5> cov{}; + auto& x() { return m_x; } + auto const& x() const { return m_x; } + auto& y() { return m_y; } + auto const& y() const { return m_y; } + auto& z() { return m_z; } + auto const& z() const { return m_z; } + auto& qOverP() { return m_qOverP; } + auto const& qOverP() const { return m_qOverP; } + auto& tx() { return m_tx; } + auto const& tx() const { return m_tx; } + auto& ty() { return m_ty; } + auto const& ty() const { return m_ty; } + + template <typename T> + MiniState( T const& other ) + : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} { + static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> ); + auto const& other_cov = other.covariance(); + LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { + LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } ); + } ); + } + + void linearTransportTo( FloatType const& new_z ) { + const auto dz = new_z - m_z; + const auto dz2 = dz * dz; + m_x += dz * m_tx; + m_y += dz * m_ty; + cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 ); + cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) ); + cov( 2, 0 ) += dz * cov( 2, 2 ); + cov( 3, 0 ) += dz * cov( 3, 2 ); + cov( 4, 0 ) += dz * cov( 4, 2 ); + cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 ); + cov( 2, 1 ) += dz * cov( 3, 2 ); + cov( 3, 1 ) += dz * cov( 3, 3 ); + cov( 4, 1 ) += dz * cov( 4, 3 ); + m_z = new_z; + } + + LHCb::LinAlg::MatSym<FloatType, 2> covXX() const { + return cov.template sub<LHCb::LinAlg::MatSym<FloatType, 2>, 0, 0>(); + } + LHCb::LinAlg::Mat<FloatType, 2> covXT() const { return cov.template sub<LHCb::LinAlg::Mat<FloatType, 2>, 2, 0>(); } + + private: + FloatType m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{}; + }; + + /** State without covariance matrix information + */ + template <typename FloatType> + class UnfittedState { + + public: + UnfittedState( FloatType x, FloatType y, FloatType z, FloatType tx, FloatType ty, FloatType qOverP ) + : m_x{std::move( x )} + , m_y{std::move( y )} + , m_z{std::move( z )} + , m_tx{std::move( tx )} + , m_ty{std::move( ty )} + , m_qOverP{std::move( qOverP )} {} + + template <class StateType> + UnfittedState( StateType const& other ) + : UnfittedState( other.x(), other.y(), other.z(), other.tx(), other.ty(), other.qOverP() ) {} + + auto& x() { return m_x; } + auto const& x() const { return m_x; } + auto& y() { return m_y; } + auto const& y() const { return m_y; } + auto& z() { return m_z; } + auto const& z() const { return m_z; } + auto& qOverP() { return m_qOverP; } + auto const& qOverP() const { return m_qOverP; } + auto& tx() { return m_tx; } + auto const& tx() const { return m_tx; } + auto& ty() { return m_ty; } + auto const& ty() const { return m_ty; } + + private: + FloatType m_x, m_y, m_z, m_tx, m_ty, m_qOverP; + }; + /** @fn stateVectorFromComposite * @brief Convert a particle (3-momentum + position) to 4-state (at some z). * @@ -95,8 +186,6 @@ namespace Sel { /** @fn stateVectorFromNeutral @brief Convert a particle (4-momentum + position) to 4-state (at some z). - - TODO: This must be properly filled */ template <typename particle_t> auto stateVectorFromNeutral( particle_t const& p ) { @@ -105,16 +194,12 @@ namespace Sel { auto const invpz = 1.f / p.pz(); - State4<float_v> s{}; - s.x() = p.hypothesis_position_x(); - s.y() = p.hypothesis_position_y(); - s.z() = p.hypothesis_position_z(); - s.tx() = p.px() * invpz; - s.ty() = p.py() * invpz; - - s.covXX() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>(); - s.covTT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::MatSym<float_v, 2>>(); - s.covXT() = LHCb::LinAlg::initialize_with_zeros<LHCb::LinAlg::Mat<float_v, 2, 2>>(); + UnfittedState<float_v> s{p.hypothesis_position_x(), + p.hypothesis_position_y(), + p.hypothesis_position_z(), + p.px() * invpz, + p.py() * invpz, + 1. / p.p()}; return s; } diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index 400448403c8..d2849f21967 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -33,9 +33,6 @@ namespace { template <std::size_t N, std::size_t M, typename T> using Matrix = LHCb::LinAlg::Mat<T, N, M>; - template <std::size_t N, typename T> - using MatrixNxN = LHCb::LinAlg::Mat<T, N>; - // Fast utility function for inverting a 2x2 matrix. This just comes // from ROOT/Math/Dinv, but I took out the test on the determinant // since we do not need it. @@ -63,79 +60,110 @@ namespace { NumTypes }; - template <ChildType C> - static constexpr bool has_full_state_v = ( C == ChildType::Composite || C == ChildType::Resonance || - C == ChildType::Neutral || C == ChildType::Other ); + /** Define the type of child based on the properties of a proxy - /* TODO look at merging this with the stuff in State4.h + Note that it is hardly possible to define functions and classes based on the + proxy types due to the presence of zip, cv-qualifiers and the high volatility + of the containers in their definition. We rely on static member values that + determine which operations can done and what information can be accessed. */ - template <typename F> - struct MiniState { - SymNxN<5, F> cov{}; - auto& x() { return m_x; } - auto const& x() const { return m_x; } - auto& y() { return m_y; } - auto const& y() const { return m_y; } - auto& z() { return m_z; } - auto const& z() const { return m_z; } - auto& qOverP() { return m_qOverP; } - auto const& qOverP() const { return m_qOverP; } - auto& tx() { return m_tx; } - auto const& tx() const { return m_tx; } - auto& ty() { return m_ty; } - auto const& ty() const { return m_ty; } + template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition> + struct child_type_for_conditions; - template <typename T> - MiniState( T const& other ) - : m_x{other.x()}, m_y{other.y()}, m_tx{other.tx()}, m_ty{other.ty()}, m_z{other.z()}, m_qOverP{other.qOverP()} { - static_assert( Sel::Utils::canBeExtrapolatedDownstream_v<T> ); - auto const& other_cov = other.covariance(); - LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { - LHCb::Utils::unwind<i, 5>( [this, i, &other_cov]( auto j ) { cov( i, j ) = other_cov( i, j ); } ); - } ); - } + template <> + struct child_type_for_conditions<true /* particle */, true /* track */, true /* extrapolation */> { + static constexpr auto value = ChildType::TrackWithVelo; + }; - void linearTransportTo( F const& new_z ) { - const auto dz = new_z - m_z; - const auto dz2 = dz * dz; - m_x += dz * m_tx; - m_y += dz * m_ty; - cov( 0, 0 ) += dz2 * cov( 2, 2 ) + 2 * dz * cov( 2, 0 ); - cov( 1, 0 ) += dz2 * cov( 3, 2 ) + dz * ( cov( 3, 0 ) + cov( 2, 1 ) ); - cov( 2, 0 ) += dz * cov( 2, 2 ); - cov( 3, 0 ) += dz * cov( 3, 2 ); - cov( 4, 0 ) += dz * cov( 4, 2 ); - cov( 1, 1 ) += dz2 * cov( 3, 3 ) + 2 * dz * cov( 3, 1 ); - cov( 2, 1 ) += dz * cov( 3, 2 ); - cov( 3, 1 ) += dz * cov( 3, 3 ); - cov( 4, 1 ) += dz * cov( 4, 3 ); - m_z = new_z; - } + template <> + struct child_type_for_conditions<true /* particle */, false /* track */, false /* extrapolation */> { + static constexpr auto value = ChildType::Neutral; + }; + + template <> + struct child_type_for_conditions<false /* particle */, false /* track */, false /* extrapolation */> { + static constexpr auto value = ChildType::Composite; + }; + + /// Determine the child type for a given proxy + template <class Proxy> + struct child_type_for_proxy + : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>, + Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {}; + + /// Child type for a given proxy + template <class Proxy> + static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value; + + template <typename Proxy> + static constexpr auto is_composite_v = child_type_for_proxy_v<Proxy> == ChildType::Composite; + + template <typename Proxy> + static constexpr bool has_covariance_v = ( child_type_for_proxy_v<Proxy> == ChildType::Composite || + child_type_for_proxy_v<Proxy> == ChildType::Resonance || + child_type_for_proxy_v<Proxy> == ChildType::TrackWithVelo || + child_type_for_proxy_v<Proxy> == ChildType::TrackWithoutVelo || + child_type_for_proxy_v<Proxy> == ChildType::Other ); - SymNxN<2, F> covXX() const { return cov.template sub<SymNxN<2, F>, 0, 0>(); } - MatrixNxN<2, F> covXT() const { return cov.template sub<MatrixNxN<2, F>, 2, 0>(); } + /** Basic proxy for a particle and a state (no covariance matrix) + */ + template <typename ProxyType> + struct UnfittedParticle final { + + public: + using proxy_type = ProxyType; + using simd_type = typename proxy_type::simd_t; + using float_v = typename simd_type::float_v; + using state_type = Sel::UnfittedState<float_v>; private: - F m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{}; + proxy_type m_particle; + state_type m_state; + + public: + UnfittedParticle( proxy_type p, state_type state ) : m_particle{std::move( p )}, m_state{std::move( state )} {} + + /** Add the momenta of a particle with no covariance information to the result of a combination + + Update the information of the momenta of the objects based on the vertex of + the decaying particle. It is crucial that we add the objects after + determining the decay vertex (i.e. at the end of the combination). + */ + void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&, + Matrix<4, 3, float_v>& ) const { + + auto dx = m_state.x() - vertexpos( 0 ); + auto dy = m_state.y() - vertexpos( 1 ); + auto dz = m_state.z() - vertexpos( 2 ); + + auto d = sqrt( dx * dx + dy * dy + dz * dz ); + + auto p = m_particle.p(); + + p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()}; + } }; // Class for daughters with a straight-line trajectory - template <typename particle_t, ChildType child_type> - struct VertexTraj final { - using simd_t = typename particle_t::simd_t; - using int_v = typename simd_t::int_v; - using float_v = typename simd_t::float_v; - using state_t = std::conditional_t<has_full_state_v<child_type>, Sel::State4<float_v>, MiniState<float_v>>; + template <typename ProxyType> + struct FittedParticle final { + + public: + using proxy_type = ProxyType; + using simd_type = typename proxy_type::simd_t; + using int_v = typename simd_type::int_v; + using float_v = typename simd_type::float_v; + using state_type = std::conditional_t<is_composite_v<proxy_type>, Sel::State4<float_v>, Sel::MiniState<float_v>>; private: - particle_t m_particle; // TODO maybe just cache some information from this? - state_t m_state; // TODO evaluate what is really needed + proxy_type m_particle; // TODO maybe just cache some information from this? + state_type m_state; // TODO evaluate what is really needed VecN<2, float_v> m_q; // predicted/fitted slope (tx,ty) SymNxN<2, float_v> m_G; // weight matrix of (x,y) of state Matrix<2, 3, float_v> m_A; // projection matrix for vertex position public: - VertexTraj( particle_t p, state_t state ) + FittedParticle( proxy_type p, state_type state ) : m_particle{std::move( p )}, m_state{std::move( state )}, m_q{m_state.tx(), m_state.ty()} {} void project( VecN<3, float_v> const& vertexpos, VecN<3, float_v>& halfDChisqDX, SymNxN<3, float_v>& halfD2ChisqDX2, @@ -187,6 +215,9 @@ namespace { void addToFourVector( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, Matrix<4, 3, float_v>& gainmatrix ) const { + + static constexpr auto child_type = child_type_for_proxy_v<proxy_type>; + // Despatch is a bit crude, but it'll do...in the original // ParticleVertexFitter this was done with a full specialisation, but // now it would have to be a partial specialisation and it all gets a @@ -195,35 +226,10 @@ namespace { addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); } else if constexpr ( child_type == ChildType::TrackWithVelo || child_type == ChildType::TrackWithoutVelo ) { addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); - } else if constexpr ( child_type == ChildType::Neutral ) { - addToFourVector_Neutral( vertexpos, p4, p4cov, gainmatrix ); } else { } } - /** Add the momenta of a neutral particle to the result of a combination - - Update the information of the momenta of the neutral objects based on the vertex of - the decaying particle. It is crucial that we add the neutral objects after - determining the decay vertex (i.e. at the end of the combination). - - @todo FIXME correctly process the covariance matrix - */ - void addToFourVector_Neutral( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>&, - Matrix<4, 3, float_v>& ) const { - // TODO: make this a global function that modifies the neutral objects at the - // persistency stage - auto dx = m_state.x() - vertexpos( 0 ); - auto dy = m_state.y() - vertexpos( 1 ); - auto dz = m_state.z() - vertexpos( 2 ); - - auto d = sqrt( dx * dx + dy * dy + dz * dz ); - - auto p = m_particle.p(); - - p4 = VecN<4, float_v>{p * dx / d, p * dy / d, p * dz / d, m_particle.energy()}; - } - /** @todo FIXME implement logic for bremsstrahlung photons. */ void addToFourVector_Track( VecN<3, float_v> const& vertexpos, VecN<4, float_v>& p4, SymNxN<4, float_v>& p4cov, @@ -377,43 +383,16 @@ namespace { } }; - template <bool ParticleCondition, bool TrackCondition, bool ExtrapolationCondition> - struct child_type_for_conditions; - - template <> - struct child_type_for_conditions<true, true, true> { - static constexpr auto value = ChildType::TrackWithVelo; - }; - - template <> - struct child_type_for_conditions<true, false, false> { - static constexpr auto value = ChildType::Neutral; - }; - - template <> - struct child_type_for_conditions<false, false, false> { - static constexpr auto value = ChildType::Composite; - }; - - template <class Proxy> - struct child_type_for_proxy - : child_type_for_conditions<Sel::Utils::isBasicParticle_v<Proxy>, Sel::Utils::hasTrack_v<Proxy>, - Sel::Utils::canBeExtrapolatedDownstream_v<Proxy>> {}; - - template <class Proxy> - static constexpr auto child_type_for_proxy_v = child_type_for_proxy<Proxy>::value; - + /// Define the constructor of child objects for each child type template <ChildType> - struct ChildHelper; + struct prepare_child_t; /** Helper type for identifying and handling TrackWithVelo and * TrackWithoutVelo children. * @todo Handle TrackWithoutVelo. */ template <> - struct ChildHelper<ChildType::TrackWithVelo> { - - static constexpr auto child_type = ChildType::TrackWithVelo; + struct prepare_child_t<ChildType::TrackWithVelo> { /** Given that `value` was `true`, and that `type` was returned by our * `deriveType` function, return the narrowest-possible std::variant of @@ -428,11 +407,9 @@ namespace { * `state(StateLocation::ClosestToBeam)` accessor. Remove [[maybe_unused]]. */ template <typename mask_v, typename float_v, typename proxy> - static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) { - using VertexTrack = - VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>; - return VertexTrack{child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam )}; + constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const { + return FittedParticle<proxy>( child, child.state( LHCb::Event::v3::Tracks::StateLocation::ClosestToBeam ) ); } }; @@ -442,9 +419,7 @@ namespace { * to go in the variant returned by prepareChild. */ template <> - struct ChildHelper<ChildType::Composite> { - - static constexpr auto child_type = ChildType::Composite; + struct prepare_child_t<ChildType::Composite> { /** Given that `value` was `true`, and that `type` was returned by our * `deriveType` function, return the narrowest-possible std::variant of @@ -459,34 +434,31 @@ namespace { * add another member to the return type; remove [[maybe_unused]] */ template <typename mask_v, typename float_v, typename proxy> - static auto prepareChild( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) { - using VertexComposite = - VertexTraj<proxy, child_type /* flag this represents a weakly-decaying composite and not a track */>; - return VertexComposite{child, Sel::stateVectorFromComposite( child )}; + constexpr auto operator()( LHCb::LinAlg::Vec<float_v, 3> const&, mask_v const&, proxy const& child ) const { + return FittedParticle<proxy>( child, Sel::stateVectorFromComposite( child ) ); } }; /** Helper type for the default `Neutral` case. */ template <> - struct ChildHelper<ChildType::Neutral> { - - static constexpr auto child_type = ChildType::Neutral; + struct prepare_child_t<ChildType::Neutral> { /** @todo Need to implement the helper class that handles generic * children and return an instance of it. As a 1-member variant. */ template <typename mask_v, typename float_v, typename proxy> - static auto prepareChild( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, - [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) { - using VertexNeutral = - VertexTraj<proxy, child_type /* flag this represents a track not a weakly-decaying composite */>; - return VertexNeutral{child, Sel::stateVectorFromNeutral( child )}; + constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec<float_v, 3> const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const { + return UnfittedParticle<proxy>{child, Sel::stateVectorFromNeutral( child )}; } }; - template <class Proxy> - using child_helper_for_proxy_t = ChildHelper<child_type_for_proxy_v<Proxy>>; + /// Return the object that will handle the information of a particle and its state + template <class Proxy, typename... Args> + auto prepare_child( Args&&... args ) { + return prepare_child_t<child_type_for_proxy_v<Proxy>>{}( std::forward<Args>( args )... ); + } } // namespace namespace Sel::Fitters { @@ -559,8 +531,8 @@ namespace Sel::Fitters { template <typename mask_v, typename float_v, template <class...> class Variant, typename... Proxy> struct return_types_helper<mask_v, float_v, Variant<Proxy...>> { - using type = std::variant<decltype( child_helper_for_proxy_t<Proxy>::prepareChild( - std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>; + using type = std::variant<decltype( + prepare_child<Proxy>( std::declval<VecN<3, float_v>>(), std::declval<mask_v>(), std::declval<Proxy>() ) )...>; }; template <typename mask_v, typename float_v, typename VariantTypes> @@ -580,15 +552,10 @@ namespace Sel::Fitters { // if the former, this yields tuple<Ts...>, otherwise tuple<child_t> using variant_types = typename LHCb::is_variant<child_t>::tuple_type; - using ret_t = return_types_helper_t<mask_v, float_v, variant_types>; - - static_assert( boost::mp11::mp_is_set<ret_t>::value, - "Got duplicate entries in derived variant type. Looks like the code needs an mp_unique<...>." ); // Now we know what return type to impose, invoke/visit on `child` return LHCb::invoke_or_visit( - [&]( auto const& child ) -> ret_t { - using helper_t = child_helper_for_proxy_t<std::decay_t<decltype( child )>>; - return helper_t::prepareChild( vertex_position, vertex_position_mask, child ); + [&]( auto const& child ) -> return_types_helper_t<mask_v, float_v, variant_types> { + return prepare_child<std::decay_t<decltype( child )>>( vertex_position, vertex_position_mask, child ); }, possibly_variant_child ); } @@ -800,9 +767,12 @@ namespace Sel::Fitters { auto halfDChisqDX = LHCb::LinAlg::initialize_with_zeros<VecN<3, float_v>>(); // add all particles that contribute to the vertex chi2 - ( std::visit( [&, &vertex_position = vertex_position]( - auto& child ) { child.project( vertex_position, halfDChisqDX, halfD2ChisqDX2, chi2, ndof ); }, - std::get<IChild>( extrapolated_children ) ), + ( std::visit( + [&, &vertex_position = vertex_position]( auto& child ) { + if constexpr ( has_covariance_v<typename std::decay_t<decltype( child )>::proxy_type> ) + child.project( vertex_position, halfDChisqDX, halfD2ChisqDX2, chi2, ndof ); + }, + std::get<IChild>( extrapolated_children ) ), ... ); // calculate the covariance and the change in the position @@ -822,8 +792,12 @@ namespace Sel::Fitters { // again"? Right now then if we hit max_iters we might update the // slopes but not the chi2 if ( !all( converged ) ) { - ( std::visit( [& vertex_position = vertex_position]( auto& child ) { child.updateSlopes( vertex_position ); }, - std::get<IChild>( extrapolated_children ) ), + ( std::visit( + [& vertex_position = vertex_position]( auto& child ) { + if constexpr ( has_covariance_v<typename std::decay_t<decltype( child )>::proxy_type> ) + child.updateSlopes( vertex_position ); + }, + std::get<IChild>( extrapolated_children ) ), ... ); } } -- GitLab From aac7f64b528a4f0496a0771b849161f3a820bd73 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 1 Apr 2022 11:29:17 +0100 Subject: [PATCH 24/26] Fixes since there is no implementation of Resonance objects with the new event model yet --- .../include/VertexFit/ParticleVertexFitter.h | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index d2849f21967..48aefc23619 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -564,6 +564,9 @@ namespace Sel::Fitters { * 1. if there are Resonance daughters, take the vertex position of the first. * 2. if not, use the states of the first two TrackWithVelo or Composite * 3. if not, try with TrackWithoutVelo and trajpoca + + TODO: Remove the case 1-bis once the differentiation between composites and + resonances is done */ template <typename simd_t, std::size_t N, typename Combination> std::tuple<VecN<3, typename simd_t::float_v>, typename simd_t::mask_v> @@ -604,6 +607,28 @@ namespace Sel::Fitters { } ); assert( tmp.has_value() ); return *tmp; // no validity check, but we have the assert just above + } else if ( count( ChildType::TrackWithVelo ) == 0 && count( ChildType::Composite ) == 1 ) { + // Case 1-bis + std::optional<std::tuple<VecN<3, float_v>, mask_v>> tmp{std::nullopt}; + LHCb::Utils::unwind<0, N>( [&]( auto i ) { + // Find the first `Composite` child + if ( !tmp.has_value() && types[i] == ChildType::Composite ) { + // Use invoke_or_visit in case the relevant entry in `children` is + // a variant type. + LHCb::invoke_or_visit( + [&]( auto const& child ) { + using Sel::Utils::endVertexPos; + if constexpr ( !Sel::Utils::isBasicParticle_v<decltype( child )> ) { + tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} ); + } else { + throw exception( "Resonance type didn't have an endVertex()" ); + } + }, + children.template get<i>() ); + } + } ); + assert( tmp.has_value() ); + return *tmp; // no validity check, but we have the assert just above } else if ( count( ChildType::TrackWithVelo ) + count( ChildType::Composite ) >= 2 ) { // Case 2 // Use the POCA of TrackWithVelo + closestToBeamState or Composite + position/momentum -- GitLab From cd1a067096d14ee5c9ae6c58e226ae976390e192 Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 1 Apr 2022 17:09:13 +0100 Subject: [PATCH 25/26] Make the pointing strategy an enum-based configurable --- Phys/ParticleMaker/src/NeutralMakers.cpp | 26 +++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 1fadde1eb4f..6f6997ebfc4 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -22,6 +22,7 @@ #include "CaloFutureUtils/CaloMomentum.h" #include "Kernel/IParticlePropertySvc.h" #include "Kernel/ParticleProperty.h" +#include "Kernel/meta_enum.h" // File containing definition AND implementation of neutral makers (photons, pi0s (...)) /** @class PhotonMaker @@ -637,6 +638,8 @@ namespace LHCb::Phys::ParticleMakers { debug() << "--------------------" << endmsg; } + meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV ); + /** @class NeutralBasicsMaker * * @brief Create a container of neutral objects from the output of the reconstruction sequence and a set of vertices @@ -660,7 +663,10 @@ namespace LHCb::Phys::ParticleMakers { LHCb::Event::NeutralBasics neutrals( unique_id_gen ); - if ( m_AveragePVs.value() ) { + switch ( m_PointingStrategy.value() ) { + + case PointingStrategy::AveragePVs: { + // Calculate the average the primary vertex positions neutrals.reserve( calo_hypos.size() ); @@ -672,21 +678,31 @@ namespace LHCb::Phys::ParticleMakers { } ); for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, avg ); - } else { + + break; + } + case PointingStrategy::OneCandidatePerPV: { + // Create one particle per primary vertex and ECAL hypothesis neutrals.reserve( calo_hypos.size() * primary_vertices.size() ); for ( auto const& ch : calo_hypos ) for ( auto const& vtx : primary_vertices ) neutrals.emplace_back( unique_id_gen, ch, vtx.position() ); + + break; + } + case PointingStrategy::Unknown: + // we will never end-up here since Gaudi::Property will trigger an error beforehand + throw std::runtime_error( "Pointing strategy set to \"Unknown\" in NeutralBasicsMaker" ); } return neutrals; } private: - Gaudi::Property<bool> m_AveragePVs{this, "AveragePVs", true, - "If set to \"true\"; uses the average of the primary vertices in order to " - "determine the pointing direction of the neutral particles"}; + Gaudi::Property<PointingStrategy> m_PointingStrategy{ + this, "PointingStrategy", PointingStrategy::AveragePVs, + "Strategy to adopt in order to determine the direction of the momenta"}; }; DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker ) -- GitLab From d267ae0728355131f194d3630d02237c95add02b Mon Sep 17 00:00:00 2001 From: CI Runner <ci-runner@gitlab.cern> Date: Fri, 1 Apr 2022 17:16:49 +0100 Subject: [PATCH 26/26] Add a new option to NeutralBasicsMaker --- Phys/ParticleMaker/src/NeutralMakers.cpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/Phys/ParticleMaker/src/NeutralMakers.cpp b/Phys/ParticleMaker/src/NeutralMakers.cpp index 6f6997ebfc4..032040c6ffb 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -638,7 +638,7 @@ namespace LHCb::Phys::ParticleMakers { debug() << "--------------------" << endmsg; } - meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV ); + meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV, OriginAsPV ); /** @class NeutralBasicsMaker * @@ -667,7 +667,7 @@ namespace LHCb::Phys::ParticleMakers { case PointingStrategy::AveragePVs: { - // Calculate the average the primary vertex positions + // calculate the average the primary vertex positions neutrals.reserve( calo_hypos.size() ); auto avg = std::accumulate( primary_vertices.cbegin(), primary_vertices.cend(), Gaudi::XYZPoint{0.f, 0.f, 0.f}, @@ -683,7 +683,7 @@ namespace LHCb::Phys::ParticleMakers { } case PointingStrategy::OneCandidatePerPV: { - // Create one particle per primary vertex and ECAL hypothesis + // create one particle per primary vertex and ECAL hypothesis neutrals.reserve( calo_hypos.size() * primary_vertices.size() ); for ( auto const& ch : calo_hypos ) @@ -691,6 +691,17 @@ namespace LHCb::Phys::ParticleMakers { break; } + case PointingStrategy::OriginAsPV: { + + // all particles have the (0, 0, 0) as the origin + neutrals.reserve( calo_hypos.size() ); + + Gaudi::XYZPoint origin{0.f, 0.f, 0.f}; + + for ( auto const& ch : calo_hypos ) neutrals.emplace_back( unique_id_gen, ch, origin ); + + break; + } case PointingStrategy::Unknown: // we will never end-up here since Gaudi::Property will trigger an error beforehand throw std::runtime_error( "Pointing strategy set to \"Unknown\" in NeutralBasicsMaker" ); -- GitLab