diff --git a/Muon/MuonID/python/MuonID/ConfiguredMuonIDs.py b/Muon/MuonID/python/MuonID/ConfiguredMuonIDs.py index 9026e77500e74531fd39c12530b57c267f1b9aaf..f689b6ca535358affa48b72b125f1ae6bec0e4b2 100644 --- a/Muon/MuonID/python/MuonID/ConfiguredMuonIDs.py +++ b/Muon/MuonID/python/MuonID/ConfiguredMuonIDs.py @@ -51,7 +51,7 @@ class ConfiguredMuonIDs(object): ## check if modules exist and load them locals = {} try: - exec ("from MuonID import " + mod[0] + " as info", {}, locals) + exec("from MuonID import " + mod[0] + " as info", {}, locals) except: if version != "def": log.warning( @@ -62,9 +62,9 @@ class ConfiguredMuonIDs(object): "ConfiguredMuonIDs: Default seems not available for DATA=%s. Loading older default %s" % (data, mod[2])) try: - exec ("from MuonID import " + mod[1] + " as info", {}, locals) + exec("from MuonID import " + mod[1] + " as info", {}, locals) except: - exec ("from MuonID import " + mod[2] + " as info", {}, locals) + exec("from MuonID import " + mod[2] + " as info", {}, locals) info = locals["info"] diff --git a/Muon/MuonID/python/MuonID/loadModule.py b/Muon/MuonID/python/MuonID/loadModule.py index a85b7cfb6f21abd12dd24c40ec8fa92ef3f19bdf..43c0fbf42f29f489095abe9dc33ebf96051b1c15 100644 --- a/Muon/MuonID/python/MuonID/loadModule.py +++ b/Muon/MuonID/python/MuonID/loadModule.py @@ -25,5 +25,5 @@ def loadFromModule(file): myf = str(file) locals = {} - exec ("from MuonID import " + myf, {}, locals) + exec("from MuonID import " + myf, {}, locals) return getattr(locals[myf], myf) diff --git a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h index 17af73e8ca90d0cd6807d43f18ecc704b7ab184f..252db5d0eff341ccf8dcee4d2c7844bb8e7a7fb7 100644 --- a/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h +++ b/Phys/ParticleCombiners/include/CombKernel/ThOrCombiner.h @@ -169,11 +169,10 @@ namespace ThOr::detail::Combiner { template auto have_overlap( T0 const& p0, T1 const& p1 ) { - if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v && Sel::Utils::canBeExtrapolatedDownstream_v ) { + if constexpr ( Sel::Utils::isBasicParticle_v && Sel::Utils::isBasicParticle_v ) { // track-like/track-like return p0.unique_id() == p1.unique_id(); - } else if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v && - !Sel::Utils::canBeExtrapolatedDownstream_v ) { + } else if constexpr ( Sel::Utils::isBasicParticle_v && !Sel::Utils::isBasicParticle_v ) { // 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 && - Sel::Utils::canBeExtrapolatedDownstream_v ) { + } else if constexpr ( !Sel::Utils::isBasicParticle_v && Sel::Utils::isBasicParticle_v ) { // composite/track-like return have_overlap( p1, p0 ); } else { diff --git a/Phys/ParticleCombiners/src/Filter.cpp b/Phys/ParticleCombiners/src/Filter.cpp index d87b7eee794f86737c5f250c961a338c5dbbcccf..db12706628b5c106464e888fb3d121e447454656 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, "ChargedBasicsFilter" ) +DECLARE_COMPONENT_WITH_ID( Pr::Filter, "NeutralBasicsFilter" ) + DECLARE_COMPONENT_WITH_ID( Pr::Filter, "ParticleRangeFilter" ) diff --git a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp index b90e62625d7ec637faf684497e096c43d9fd3448..57808bdc7507f2d75859aaa964b26d529e5e23b5 100644 --- a/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp +++ b/Phys/ParticleCombiners/src/ThOrCombiner_Composites.cpp @@ -19,6 +19,10 @@ using ThOrCombiner__Composites3ChargedBasics = ThOr::CombinerBest; +using ThOrCombiner__CompositesNeutralBasics = ThOr::CombinerBest; + 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 43d39bdea1ac8a234b5f659d406c11d36374a3f5..032040c6ffb7ec8db691e904782c55b828e967a5 100644 --- a/Phys/ParticleMaker/src/NeutralMakers.cpp +++ b/Phys/ParticleMaker/src/NeutralMakers.cpp @@ -14,12 +14,15 @@ #include "GaudiAlg/Transformer.h" #include "Event/Particle.h" +#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" #include "Kernel/ParticleProperty.h" +#include "Kernel/meta_enum.h" // File containing definition AND implementation of neutral makers (photons, pi0s (...)) /** @class PhotonMaker @@ -634,4 +637,85 @@ namespace LHCb::Phys::ParticleMakers { debug() << " Skipped : " << m_SkipPi0sCounter << endmsg; debug() << "--------------------" << endmsg; } + + meta_enum_class( PointingStrategy, int, Unknown = -1, AveragePVs, OneCandidatePerPV, OriginAsPV ); + + /** @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 { + + 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", ""}} ) {} + + 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 { + + LHCb::Event::NeutralBasics neutrals( unique_id_gen ); + + switch ( m_PointingStrategy.value() ) { + + case PointingStrategy::AveragePVs: { + + // 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}, + [&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 ); + + 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::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" ); + } + + return neutrals; + } + + private: + Gaudi::Property m_PointingStrategy{ + this, "PointingStrategy", PointingStrategy::AveragePVs, + "Strategy to adopt in order to determine the direction of the momenta"}; + }; + + DECLARE_COMPONENT( LHCb::Phys::ParticleMakers::NeutralBasicsMaker ) + } // namespace LHCb::Phys::ParticleMakers diff --git a/Phys/SelAlgorithms/src/Monitor.cpp b/Phys/SelAlgorithms/src/Monitor.cpp index 8154407d84103f2a0fa411f6821a766588666e40..59230867ba7a07db8f812ec6f8953dce6f6ec370 100644 --- a/Phys/SelAlgorithms/src/Monitor.cpp +++ b/Phys/SelAlgorithms/src/Monitor.cpp @@ -31,6 +31,8 @@ namespace SelAlgorithms { DECLARE_COMPONENT_WITH_ID( Monitor, "Monitor__ChargedBasics" ) + DECLARE_COMPONENT_WITH_ID( Monitor, "Monitor__NeutralBasics" ) + DECLARE_COMPONENT_WITH_ID( Monitor, "Monitor__Composites" ) DECLARE_COMPONENT_WITH_ID( Monitor, "Monitor__ParticleRange" ) diff --git a/Phys/SelKernel/include/SelKernel/Utilities.h b/Phys/SelKernel/include/SelKernel/Utilities.h index f1c30e4f6ce6e313354e3bce267b2c1ca0bb52ea..7422842d3a91fbfdb4e81101c51c5c9c4d195fc5 100644 --- a/Phys/SelKernel/include/SelKernel/Utilities.h +++ b/Phys/SelKernel/include/SelKernel/Utilities.h @@ -294,6 +294,12 @@ namespace Sel::Utils { template constexpr auto canBeExtrapolatedDownstream_v = std::remove_pointer_t>::canBeExtrapolatedDownstream; + template + constexpr auto isBasicParticle_v = std::remove_pointer_t>::isBasicParticle; + + template + constexpr auto hasTrack_v = std::remove_pointer_t>::hasTrack; + template 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 567e36191707be4c99a0c3cbeb2a0521d8f91004..f2f9f847707886bd0f493a55d38cf61e1629a2ff 100644 --- a/Phys/SelTools/include/SelTools/Utilities.h +++ b/Phys/SelTools/include/SelTools/Utilities.h @@ -9,10 +9,103 @@ * or submit itself to any jurisdiction. * \*****************************************************************************/ #pragma once +#include "Event/Particle_v2.h" #include "LHCbMath/MatVec.h" #include "SelKernel/State4.h" 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 + struct MiniState { + LHCb::LinAlg::MatSym 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 + 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 ); + auto const& other_cov = other.covariance(); + LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { + LHCb::Utils::unwind( [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 covXX() const { + return cov.template sub, 0, 0>(); + } + LHCb::LinAlg::Mat covXT() const { return cov.template sub, 2, 0>(); } + + private: + FloatType m_x{}, m_y{}, m_tx{}, m_ty{}, m_z{}, m_qOverP{}; + }; + + /** State without covariance matrix information + */ + template + 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 + 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). * @@ -90,4 +183,24 @@ namespace Sel { } return s; } + + /** @fn stateVectorFromNeutral + @brief Convert a particle (4-momentum + position) to 4-state (at some z). + */ + template + auto stateVectorFromNeutral( particle_t const& p ) { + + using float_v = std::decay_t; + + auto const invpz = 1.f / p.pz(); + + UnfittedState s{p.hypothesis_position_x(), + p.hypothesis_position_y(), + p.hypothesis_position_z(), + p.px() * invpz, + p.py() * invpz, + 1. / p.p()}; + + return s; + } } // namespace Sel diff --git a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h index d47095ffb215cb8fe311cc9e45e4cab7e3bd47ee..48aefc23619082038158d75974c10bef8b20a0ed 100644 --- a/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h +++ b/Phys/VertexFit/include/VertexFit/ParticleVertexFitter.h @@ -22,6 +22,445 @@ #include "Gaudi/Algorithm.h" +namespace { + + template + using VecN = LHCb::LinAlg::Vec; + + template + using SymNxN = LHCb::LinAlg::MatSym; + + template + using Matrix = LHCb::LinAlg::Mat; + + // 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 + 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 + }; + + /** Define the type of child based on the properties of a proxy + + 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 + struct child_type_for_conditions; + + template <> + struct child_type_for_conditions { + static constexpr auto value = ChildType::TrackWithVelo; + }; + + template <> + struct child_type_for_conditions { + static constexpr auto value = ChildType::Neutral; + }; + + template <> + struct child_type_for_conditions { + static constexpr auto value = ChildType::Composite; + }; + + /// Determine the child type for a given proxy + template + struct child_type_for_proxy + : child_type_for_conditions, Sel::Utils::hasTrack_v, + Sel::Utils::canBeExtrapolatedDownstream_v> {}; + + /// Child type for a given proxy + template + static constexpr auto child_type_for_proxy_v = child_type_for_proxy::value; + + template + static constexpr auto is_composite_v = child_type_for_proxy_v == ChildType::Composite; + + template + static constexpr bool has_covariance_v = ( child_type_for_proxy_v == ChildType::Composite || + child_type_for_proxy_v == ChildType::Resonance || + child_type_for_proxy_v == ChildType::TrackWithVelo || + child_type_for_proxy_v == ChildType::TrackWithoutVelo || + child_type_for_proxy_v == ChildType::Other ); + + /** Basic proxy for a particle and a state (no covariance matrix) + */ + template + 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; + + private: + 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 + 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, Sel::State4, Sel::MiniState>; + + private: + 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: + 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, + 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 Vba = m_state.covariance().template + // Sub>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit, + // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix 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 { + + static constexpr auto child_type = child_type_for_proxy_v; + + // 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 { + } + } + + /** @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, 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, 2, 2>() ) - + LHCb::LinAlg::similarity( FK, m_state.cov.template sub, 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 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 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 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; + } + }; + + /// Define the constructor of child objects for each child type + template + struct prepare_child_t; + + /** Helper type for identifying and handling TrackWithVelo and + * TrackWithoutVelo children. + * @todo Handle TrackWithoutVelo. + */ + template <> + struct prepare_child_t { + + /** 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 + constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const { + return FittedParticle( 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 prepare_child_t { + + /** 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 + constexpr auto operator()( LHCb::LinAlg::Vec const&, mask_v const&, proxy const& child ) const { + return FittedParticle( child, Sel::stateVectorFromComposite( child ) ); + } + }; + + /** Helper type for the default `Neutral` case. + */ + template <> + struct prepare_child_t { + + /** @todo Need to implement the helper class that handles generic + * children and return an instance of it. As a 1-member variant. + */ + template + constexpr auto operator()( [[maybe_unused]] LHCb::LinAlg::Vec const& vertex_position, + [[maybe_unused]] mask_v const& vertex_position_mask, proxy const& child ) const { + return UnfittedParticle{child, Sel::stateVectorFromNeutral( child )}; + } + }; + + /// Return the object that will handle the information of a particle and its state + template + auto prepare_child( Args&&... args ) { + return prepare_child_t>{}( std::forward( args )... ); + } +} // namespace + namespace Sel::Fitters { /** @class Sel::Fitters::ParticleVertex @@ -44,86 +483,6 @@ namespace Sel::Fitters { } private: - template - using VecN = LHCb::LinAlg::Vec; - - template - using SymNxN = LHCb::LinAlg::MatSym; - - template - using Matrix = LHCb::LinAlg::Mat; - - template - using MatrixNxN = LHCb::LinAlg::Mat; - - /* TODO look at merging this with the stuff in State4.h - */ - template - 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 - 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 ); - auto const& other_cov = other.covariance(); - LHCb::Utils::unwind<0, 5>( [this, &other_cov]( auto i ) { - LHCb::Utils::unwind( [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, 0, 0>(); } - MatrixNxN<2, F> covXT() const { return cov.template sub, 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, - Other, - NumTypes - }; - /** Indirectly stolen from TrackVertexUtils::poca */ template @@ -164,413 +523,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 - 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 - 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, MiniState>; - - 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 Vba = m_state.covariance().template - // Sub>(2,0); compute the corresponding gain matrix (this is WBG in the BFR fit, - // but here it is Vba * Vaa^-1) ROOT::Math::SMatrix 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 ( represents_composite ) { - addToFourVector_Composite( vertexpos, p4, p4cov, gainmatrix ); - } else { - addToFourVector_Track( vertexpos, p4, p4cov, gainmatrix ); - } - } - - /** @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, 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, 2, 2>() ) - - LHCb::LinAlg::similarity( FK, m_state.cov.template sub, 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 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 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 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 - struct TrackChildHelper { - /** 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; - - /** 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 - 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 ); - using VertexTrack = - VertexTraj; - 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 CompositeChildHelper { - /** Determine if the given `child_t` represents a [chunk of] composite - * particles. - */ - constexpr static bool value = !Sel::Utils::canBeExtrapolatedDownstream_v; - - /** 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 - 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; - return VertexComposite{child, stateVectorFromComposite( child )}; - } - }; - - /** Helper type for the default `Other` case. - */ - template - 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::Other; } - /** @todo Need to implement the helper class that handles generic - * children and return an instance of it. As a 1-member variant. - */ - template - 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 ) { - assert( type == ChildType::Other ); - return bool{}; - } - }; - - /** Master list of helper types. The order matters; the **first** matching - * helper is preferred. - */ - template - using all_helper_types = - boost::mp11::mp_list, CompositeChildHelper, LeftoverHelper>; - - /** Helper types that match (::value == true) for a given child type. - */ - template - using matching_helper_types = boost::mp11::mp_copy_if, boost::mp11::mp_to_bool>; - - /** The first matching helper type for a given child type. - */ - template - using first_matching_helper_type = boost::mp11::mp_front>; - - /** 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 - ChildType deriveType( child_t const& child ) const { - assert( m_pp_svc ); - using helper_t = first_matching_helper_type; - return helper_t::deriveType( *m_pp_svc, child ); - } - /** Helper for getting the return type of the prepareChild method of a * helper type. */ - template - struct prepareChild_helper { - template - using fn = decltype( helper_t::prepareChild( std::declval>(), std::declval(), - std::declval(), std::declval() ) ); + template + struct return_types_helper; + + template class Variant, typename... Proxy> + struct return_types_helper> { + using type = std::variant( std::declval>(), std::declval(), std::declval() ) )...>; }; + template + using return_types_helper_t = typename return_types_helper::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 @@ -579,28 +546,16 @@ namespace Sel::Fitters { * @todo Include the `vertexinsidefoil` logic from ParticleVertexFitter.cpp */ template - 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, or it might be a value type itself. // if the former, this yields tuple, otherwise tuple using variant_types = typename LHCb::is_variant::tuple_type; - // get the [first] matching helper for each element in `variant_types` - using matched_helpers = boost::mp11::mp_transform; - // 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> - using return_types = - boost::mp11::mp_transform_q, 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 - using ret_t = boost::mp11::mp_rename; - static_assert( boost::mp11::mp_is_set::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>; - return helper_t::prepareChild( vertex_position, vertex_position_mask, type, child ); + [&]( auto const& child ) -> return_types_helper_t { + return prepare_child>( vertex_position, vertex_position_mask, child ); }, possibly_variant_child ); } @@ -609,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 std::tuple, typename simd_t::mask_v> @@ -638,7 +596,29 @@ namespace Sel::Fitters { LHCb::invoke_or_visit( [&]( auto const& child ) { using Sel::Utils::endVertexPos; - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v ) { + if constexpr ( !Sel::Utils::isBasicParticle_v ) { + tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} ); + } else { + throw exception( "Resonance type didn't have an endVertex()" ); + } + }, + children.template get() ); + } + } ); + 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, 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 ) { tmp.emplace( LHCb::LinAlg::convert( endVertexPos( child ) ), mask_v{true} ); } else { throw exception( "Resonance type didn't have an endVertex()" ); @@ -660,18 +640,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 ) { + + if constexpr ( Sel::Utils::isBasicParticle_v ) { + if constexpr ( Sel::Utils::canBeExtrapolatedDownstream_v ) velotrajs[nvelotrajs++] = {trackState( child )}; - } else { - throw exception( - "TrackWithVelo type didn't have a state( child_t::StateLocation::ClosestToBeam)" ); - } - } else if ( types[i] == ChildType::Composite ) { - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v ) { + } else { + if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v ) { velotrajs[nvelotrajs++] = {referencePoint( child ), threeMomentum( child )}; } else { throw exception( "Composite type didn't one/both of referencePoint() and threeMomentum()" ); @@ -700,7 +679,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() ); } } @@ -736,13 +715,15 @@ 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 // of runtime dispatch. - std::array const types{LHCb::invoke_or_visit( [this]( auto const& child ) { return this->deriveType( child ); }, - children.template get() )...}; + std::array const types{LHCb::invoke_or_visit( + [this]( auto const& child ) { return child_type_for_proxy_v>; }, + children.template get() )...}; // Obtain an estimate of the vertex position for initialization auto [vertex_position, vertex_position_mask] = calculateInitialPosition( types, children ); @@ -759,9 +740,9 @@ 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() )...}; + prepareChild( vertex_position, prelim_fit_mask, children.template get() )...}; // These are required for populating child relations at the end, but the // API for dealing with them should be agnostic to the child types @@ -772,7 +753,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 ) + if constexpr ( !Sel::Utils::isBasicParticle_v ) return child.numDescendants(); else return 1; @@ -785,7 +766,7 @@ namespace Sel::Fitters { ( LHCb::invoke_or_visit( [¤t_unique_ids]( auto const& child ) -> void { - if constexpr ( !Sel::Utils::canBeExtrapolatedDownstream_v ) { + if constexpr ( !Sel::Utils::isBasicParticle_v ) { for ( auto i = 0u; i < child.numDescendants(); ++i ) current_unique_ids.emplace_back( child.descendant_unique_id( i ) ); } else { @@ -811,9 +792,12 @@ namespace Sel::Fitters { auto halfDChisqDX = LHCb::LinAlg::initialize_with_zeros>(); // 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( extrapolated_children ) ), + ( std::visit( + [&, &vertex_position = vertex_position]( auto& child ) { + if constexpr ( has_covariance_v::proxy_type> ) + child.project( vertex_position, halfDChisqDX, halfD2ChisqDX2, chi2, ndof ); + }, + std::get( extrapolated_children ) ), ... ); // calculate the covariance and the change in the position @@ -833,8 +817,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( extrapolated_children ) ), + ( std::visit( + [& vertex_position = vertex_position]( auto& child ) { + if constexpr ( has_covariance_v::proxy_type> ) + child.updateSlopes( vertex_position ); + }, + std::get( extrapolated_children ) ), ... ); } }