Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • rrabadan/LHCb
  • talin/LHCb
  • imjelde/LHCb
  • mstahl/LHCb
  • padeken/LHCb
  • mimazure/LHCb
  • roiser/LHCb
  • conrad/LHCb
  • kklimasz/LHCb
  • rcurrie/LHCb
  • wkrzemie/LHCb
  • fkeizer/LHCb
  • valassi/LHCb
  • hschrein/LHCb
  • anstahll/LHCb
  • jonrob/LHCb
  • graven/LHCb
  • clemenci/LHCb
  • chaen/LHCb
  • sstahl/LHCb
  • lhcb/LHCb
21 results
Show changes
Commits on Source (40)
Hide whitespace changes
Inline Side-by-side
Showing
with 678 additions and 64 deletions
Event/PhysEvent/include/Event/Particle.h
View file @ a89c7663
......@@ -13,6 +13,7 @@
// Include files
#include "Event/ProtoParticle.h"
#include "Event/VertexBase.h"
#include "GaudiKernel/GaudiException.h"
#include "GaudiKernel/GenericMatrixTypes.h"
#include "GaudiKernel/KeyedContainer.h"
#include "GaudiKernel/KeyedObject.h"
......@@ -27,6 +28,7 @@
#include "GaudiKernel/Vector4DTypes.h"
#include "GaudiKernel/VectorMap.h"
#include "Kernel/ParticleID.h"
#include "Kernel/Traits.h"
#include "LHCbMath/MatVec.h"
#include "LHCbMath/Vec3.h"
#include <algorithm>
......@@ -445,9 +447,10 @@ namespace LHCb {
friend auto pid( const Particle& p ) { return p.particleID().pid(); }
// flag which indicates client code should go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov
static constexpr bool hasVertex = true;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites), maybe
// (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::maybe;
private:
LHCb::ParticleID m_particleID; ///< PDG code
......@@ -481,6 +484,7 @@ namespace LHCb {
// -----------------------------------------------------------------------------
// Including forward declarations
#include "Event/ParticleCombination.h"
#include "Event/Vertex.h"
namespace LHCb {
......@@ -493,6 +497,30 @@ namespace LHCb {
return {pos.x(), pos.y(), pos.z()};
}
inline auto decayProducts( const LHCb::Particle& p ) {
auto const& d = p.daughtersVector();
switch ( d.size() ) {
case 2:
return ParticleCombination<LHCb::Particle>{d[0], d[1]};
case 3:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2]};
case 4:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2], d[3]};
case 5:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2], d[3], d[4]};
case 6:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2], d[3], d[4], d[5]};
case 7:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2], d[3], d[4], d[5], d[7]};
case 8:
return ParticleCombination<LHCb::Particle>{d[0], d[1], d[2], d[3], d[4], d[5], d[7], d[7]};
default:
throw GaudiException{"Combination of less than 2, or more than 8 -- this should never happen",
"decayProducts(LHCb::Particle const&)", StatusCode::FAILURE};
__builtin_unreachable();
}
}
inline Particle& Particle::operator=( const Particle& rhs ) {
if ( this != &rhs ) {
m_particleID = rhs.m_particleID;
......
Event/PhysEvent/include/Event/ParticleCombination.h 0 → 100644
View file @ a89c7663
/*****************************************************************************\
* (c) Copyright 2019-2021 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#pragma once
#include "Event/Particle.h"
#include "GaudiKernel/GaudiException.h"
#include "LHCbMath/Vec3.h"
#include "boost/container/static_vector.hpp"
#include <array>
#include <cmath>
#include <functional>
namespace LHCb {
namespace Combination::details {
using fp = float;
}
/** @brief A list of particle objects.
*
* A 'combination' is simply an ordered list of particle objects. This
* combination object holds such a list and has several member functions for
* computing kinematic 'properties' of the list.
*
* The 'four momentum' of a combination is defined as the sum of the four
* momentum of the particles held by it. Most other intrinsic properties can be
* derived from this, such as the 'invariant mass' of a combination.
*
* Some other methods exist to make this object compatible with the API expected
* by the ThOr functor framework, allowing an instance of ParticleCombination to
* be passed in as an argument to a functor.
*/
template <typename... InputTypes>
struct ParticleCombination;
template <typename Particle_t>
struct ParticleCombination<Particle_t> {
/// Pointers to each of the N child objects in this combination
/// 8 children should be more than enough...
static_assert( !std::is_pointer_v<Particle_t> );
boost::container::static_vector<Particle_t const*, 8> m_children{};
public:
template <typename... Args,
typename = std::enable_if_t<( sizeof...( Args ) > 1 ) &&
std::conjunction_v<std::is_convertible<Particle_t const*, Args>...>>>
ParticleCombination( Args&&... args ) : m_children{std::forward<Args>( args )...} {
assert( std::none_of( m_children.begin(), m_children.end(), []( auto* i ) { return i == nullptr; } ) );
}
template <typename Predicate>
bool pairwise_none_of( Predicate predicate ) const {
auto e = m_children.end();
for ( auto i1 = m_children.begin(); i1 != e; ++i1 ) {
for ( auto i2 = std::next( i1 ); i2 != e; ++i2 ) {
if ( predicate( **i1, **i2 ) ) return false;
}
}
return true;
}
template <typename Transform, typename Reduce = std::plus<>>
auto transform_reduce( Transform transform, Reduce reduce = {} ) const {
auto i1 = m_children.begin();
auto end = m_children.end();
if ( i1 == end ) {
throw GaudiException{"Empty collection of children -- this should never happen", "transform_reduce",
StatusCode::FAILURE};
}
using value_t = std::invoke_result_t<Transform, decltype( **i1 )>;
using reduced_t = std::invoke_result_t<Reduce, value_t, value_t>;
reduced_t current = transform( **i1 );
for ( ; i1 != m_children.end(); ++i1 ) { current = reduce( std::move( current ), transform( **i1 ) ); }
return current;
}
template <typename Transform, typename Reduce = std::plus<>>
auto pairwise_transform_reduce( Transform transform, Reduce reduce = {} ) const {
auto i1 = m_children.begin();
auto end = m_children.end();
if ( i1 == end ) {
throw GaudiException{"Empty collection -- this should never happen", "pairwise_transform_reduce",
StatusCode::FAILURE};
}
auto i2 = std::next( i1 );
if ( i2 == end ) {
throw GaudiException{"single particle collection -- this should never happen", "pairwise_transform_reduce",
StatusCode::FAILURE};
}
// Calculate the first value explicitly to avoid having to pass in an initial value, and thus having to specify
// its type
using value_t = std::invoke_result_t<Transform, decltype( **i1 ), decltype( **i1 )>;
using reduced_t = std::invoke_result_t<Reduce, value_t, value_t>;
reduced_t current = transform( **i1, **i2 );
do {
++i2;
for ( ; i2 != end; ++i2 ) { current = reduce( std::move( current ), transform( **i1, **i2 ) ); }
++i1;
i2 = i1;
} while ( i1 != end );
return current;
}
/** @brief Number of children.
*/
auto size() const { return m_children.size(); }
auto numChildren() const { return m_children.size(); }
Particle_t const& operator[]( size_t i ) const {
assert( i < m_children.size() );
assert( m_children[i] != nullptr );
return *m_children[i];
}
template <size_t N>
Particle_t const& get() const {
return ( *this )[N];
}
/** @brief the decayProducts of a combination _is_ the combination
*/
friend ParticleCombination const& decayProducts( ParticleCombination const& pc ) { return pc; }
auto begin() const { return m_children.begin(); }
auto end() const { return m_children.end(); }
/** @brief three momentum sum of the children.
*/
auto threeMomentum() const {
Combination::details::fp px = 0;
Combination::details::fp py = 0;
Combination::details::fp pz = 0;
for ( auto* i : m_children ) {
px += i->momentum().px();
py += i->momentum().py();
pz += i->momentum().pz();
}
return LHCb::LinAlg::Vec3{px, py, pz};
}
auto slopes() const {
auto mom = threeMomentum();
return mom / Z( mom );
}
/** @brief Squared-magnitude of the three momentum sum of the children.
*/
Combination::details::fp mom2() const { return threeMomentum().mag2(); }
template <std::size_t... idxs>
ParticleCombination subCombination() const {
assert( ( ( 0 <= idxs && idxs < m_children.size() ) && ... ) );
return ParticleCombination{m_children[idxs]...};
}
/** @brief Magnitude of the four momentum sum of the children.
*/
Combination::details::fp p() const { return std::sqrt( mom2() ); }
/** @brief Transverse momentum component of the four momentum sum of the children.
*/
Combination::details::fp pt() const {
Combination::details::fp px = 0.;
Combination::details::fp py = 0.;
for ( auto* i : m_children ) {
px += i->momentum().px();
py += i->momentum().py();
}
return std::sqrt( px * px + py * py );
}
/** @brief Invariant mass of the four momentum sum of the children.
*/
Combination::details::fp mass2() const {
// FIXME/TODO: use stable computation from https://gitlab.cern.ch/lhcb/Analysis/-/issues/23
using std::sqrt;
Combination::details::fp energy = 0.;
for ( auto* i : m_children ) { energy += sqrt( pow( i->p(), 2 ) + pow( i->measuredMass(), 2 ) ); }
return energy * energy - mom2();
}
Combination::details::fp mass() const {
using std::sqrt;
return sqrt( mass2() );
}
};
template <typename Particle_t, typename... InputTypes>
struct ParticleCombination<Particle_t, InputTypes...> : ParticleCombination<InputTypes...> {
static_assert( sizeof...( InputTypes ) > 0 );
static_assert( std::conjunction_v<std::is_same<Particle_t, InputTypes>...> );
using ParticleCombination<InputTypes...>::ParticleCombination;
};
} // namespace LHCb
Event/PhysEvent/include/Event/Particle_v2.h
View file @ a89c7663
......@@ -15,6 +15,7 @@
#include "Event/UniqueIDGenerator.h"
#include "Kernel/ParticleID.h"
#include "Kernel/STLExtensions.h"
#include "Kernel/Traits.h"
#include "Kernel/Variant.h"
#include "LHCbMath/MatVec.h"
#include "LHCbMath/SIMDWrapper.h"
......@@ -256,7 +257,11 @@ namespace LHCb::Event {
[&]( auto i ) { unwind<0, i + 1>( [&]( auto j ) { out( i, j ) = this->posCovElement( i, j ); } ); } );
return out;
}
static constexpr bool hasVertex = true;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = CanBeExtrapolatedDownstream::no;
private:
// Helper for returning bits of the momentum covariance matrix
......
Event/TrackEvent/include/Event/PrFittedForwardTracks.h
View file @ a89c7663
......@@ -16,6 +16,7 @@
#include "Event/UniqueIDGenerator.h"
#include "Kernel/EventLocalAllocator.h"
#include "Kernel/LHCbID.h"
#include "Kernel/Traits.h"
#include "LHCbMath/SIMDWrapper.h"
#include "LHCbMath/Vec3.h"
#include "TrackEnums.h"
......@@ -48,7 +49,11 @@ namespace LHCb::Pr::Fitted::Forward {
[[nodiscard, gnu::always_inline]] friend auto referencePoint( const FittedState& p ) { return p.position(); }
[[nodiscard, gnu::always_inline]] friend auto threeMomentum( const FittedState& p ) { return p.momentum(); }
static constexpr bool hasVertex = TrackProxy::hasVertex;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::yes;
private:
TrackProxy m_proxy;
......@@ -267,7 +272,8 @@ namespace LHCb::Pr::Fitted::Forward {
[[nodiscard, gnu::always_inline]] friend auto referencePoint( const FittedProxy& fp ) {
return fp.closestToBeamStatePos();
}
static constexpr bool hasVertex = false;
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::yes;
[[nodiscard, gnu::always_inline]] friend auto trackState( FittedProxy const& fp ) {
return fp.closestToBeamState();
}
......
Event/TrackEvent/include/Event/PrLongTracks.h
View file @ a89c7663
......@@ -17,6 +17,7 @@
#include "Event/PrVeloTracks.h"
#include "Kernel/EventLocalAllocator.h"
#include "Kernel/LHCbID.h"
#include "Kernel/Traits.h"
#include "Kernel/UTChannelID.h"
#include "Kernel/VPChannelID.h"
#include "LHCbMath/SIMDWrapper.h"
......@@ -175,7 +176,10 @@ namespace LHCb::Pr::Long {
[[nodiscard, gnu::always_inline]] friend auto trackState( const LongProxy& p ) {
return p.template get<Tag::States>( 0 );
}
static constexpr bool hasVertex = false;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::yes;
};
template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType>
......
Event/TrackEvent/include/Event/PrVeloTracks.h
View file @ a89c7663
......@@ -17,6 +17,7 @@
#include "Kernel/EventLocalAllocator.h"
#include "Kernel/HeaderMapping.h"
#include "Kernel/LHCbID.h"
#include "Kernel/Traits.h"
#include "Kernel/VPChannelID.h"
#include "LHCbMath/MatVec.h"
#include "LHCbMath/SIMDWrapper.h"
......@@ -138,7 +139,11 @@ namespace LHCb::Pr::Velo {
return ids;
}
static constexpr bool hasVertex = false;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::yes;
[[nodiscard, gnu::always_inline]] friend auto referencePoint( VeloProxy const& vp ) {
return vp.closestToBeamStatePos();
}
......
Event/TrackEvent/include/Event/Track_v1.h
View file @ a89c7663
......@@ -27,6 +27,7 @@
#include "Kernel/HeaderMapping.h"
#include "Kernel/LHCbID.h"
#include "Kernel/STLExtensions.h"
#include "Kernel/Traits.h"
#include "LHCbMath/Vec3.h"
#include "TrackEnums.h"
#include <algorithm>
......@@ -504,7 +505,11 @@ namespace LHCb::Event {
friend std::ostream& operator<<( std::ostream& str, const Track& obj ) { return obj.fillStream( str ); }
static constexpr bool hasVertex = false;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = CanBeExtrapolatedDownstream::yes;
[[nodiscard, gnu::always_inline]] friend auto trackState( Track const& t ) { return t.firstState(); }
[[nodiscard, gnu::always_inline]] friend auto threeMomentum( Track const& t ) {
auto m = t.momentum();
......
Event/TrackEvent/include/Event/Track_v3.h
View file @ a89c7663
......@@ -16,6 +16,7 @@
#include "Event/UniqueIDGenerator.h"
#include "Kernel/EventLocalAllocator.h"
#include "Kernel/LHCbID.h"
#include "Kernel/Traits.h"
#include "Kernel/meta_enum.h"
#include "LHCbMath/SIMDWrapper.h"
#include "LHCbMath/Vec3.h"
......@@ -148,7 +149,11 @@ namespace LHCb::Event::v3 {
[[nodiscard, gnu::always_inline]] auto z() const { return position().Z(); }
[[nodiscard, gnu::always_inline]] auto tx() const { return slopes().X(); }
[[nodiscard, gnu::always_inline]] auto ty() const { return slopes().Y(); }
static constexpr bool hasVertex = TrackProxy::hasVertex;
// flag which indicates client code can go for 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix` and not
// for a track-like stateCov -- possible values: yes (for track-like objects ) no (for neutrals, composites),
// maybe (particle, check at runtime, calls may return invalid results)
static constexpr auto canBeExtrapolatedDownstream = CanBeExtrapolatedDownstream::yes;
private:
TrackProxy m_proxy;
......@@ -321,7 +326,7 @@ namespace LHCb::Event::v3 {
[[nodiscard, gnu::always_inline]] friend auto referencePoint( FittedProxy const& fp ) {
return fp.position( fp.defaultState() );
}
static constexpr bool hasVertex = false;
static constexpr auto canBeExtrapolatedDownstream = Event::CanBeExtrapolatedDownstream::yes;
[[nodiscard, gnu::always_inline]] friend auto trackState( FittedProxy const& fp ) {
return fp.state( StateLocation::ClosestToBeam );
}
......
Kernel/LHCbKernel/include/Kernel/Traits.h 0 → 100644
View file @ a89c7663
/*****************************************************************************\
* (c) Copyright 2022 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#pragma once
// traits which indicate whether client code should use 'threeMomCovMatrix', `momPosCovMatrix` and `posCovMatrix`
// or whether to request a `trackState`
// possible values: yes (for track-like objects ) no (for neutrals, composites), maybe (particle, check at runtime, as
// calls may return invalid results -- not yet implemented)
//
namespace LHCb::Event::CanBeExtrapolatedDownstream {
struct yes_t : std::true_type {};
struct no_t : std::false_type {};
struct maybe_t : no_t {
}; // basically, if you need to know at compile time, the answer is 'no' as there is no guarantee -- but it could
// still eg. be a Particle v1 representing a track, which in turn could be extrapolated downstream..
inline constexpr auto yes = yes_t{};
inline constexpr auto no = no_t{};
inline constexpr auto maybe = maybe_t{};
} // namespace LHCb::Event::CanBeExtrapolatedDownstream
Kernel/LHCbKernel/include/Kernel/Variant.h
View file @ a89c7663
......@@ -10,16 +10,6 @@
\*****************************************************************************/
#pragma once
#include "Kernel/HeaderMapping.h"
#include <boost/version.hpp>
// ROOT/Cling has problems with std::variant in [at least] the version shipped
// in LCG_97a. Fortunately that version of LCG is new enough to contain a
// version of Boost that includes boost::variant2::variant, which ROOT/Cling
// seems to handle better.
#if BOOST_VERSION >= 107200
# define LHCB_VARIANT_USE_BOOST_VARIANT2
# include <boost/variant2/variant.hpp>
#endif
#include <tuple>
#include <variant>
......@@ -38,13 +28,6 @@ namespace LHCb {
using tuple_type = std::tuple<Ts...>;
};
#ifdef LHCB_VARIANT_USE_BOOST_VARIANT2
template <typename... Ts>
struct is_variant<boost::variant2::variant<Ts...>> : std::true_type {
using tuple_type = std::tuple<Ts...>;
};
#endif
template <typename T>
inline constexpr bool is_variant_v = is_variant<T>::value;
......@@ -52,12 +35,7 @@ namespace LHCb {
* on compile-time checks. See https://gitlab.cern.ch/lhcb/Rec/-/issues/124
*/
template <typename... Ts>
using variant =
#ifdef LHCB_VARIANT_USE_BOOST_VARIANT2
boost::variant2::variant<Ts...>;
#else
std::variant<Ts...>;
#endif
using variant = std::variant<Ts...>;
/** Check that the given list of types is not empty and contains only variant
* types.
......@@ -86,10 +64,3 @@ template <typename... Ts>
struct LHCb::header_map<std::variant<Ts...>> {
constexpr static auto value = ( LHCb::header_map_v<Ts> + ... ) + "<variant>";
};
#ifdef LHCB_VARIANT_USE_BOOST_VARIANT2
template <typename... Ts>
struct LHCb::header_map<boost::variant2::variant<Ts...>> {
constexpr static auto value = ( LHCb::header_map_v<Ts> + ... ) + "<boost/variant2/variant.hpp>";
};
#endif
Kernel/LHCbMath/CMakeLists.txt
View file @ a89c7663
......@@ -44,6 +44,7 @@ gaudi_add_library(LHCbMathLib
src/SIMDWrapper.cpp
src/Similarity.cpp
src/Spline.cpp
src/StateVertexUtils.cpp
src/ValueWithError.cpp
src/Vector3DWithError.cpp
src/WStatEntity.cpp
......
Kernel/LHCbMath/include/LHCbMath/SIMDWrapper.h
View file @ a89c7663
......@@ -247,6 +247,8 @@ namespace SIMDWrapper {
constexpr friend bool any( mask_v mask ) { return mask.data != 0; }
constexpr friend bool testbit( mask_v mask, int ) { return mask.data == 1; }
constexpr friend int ffs( mask_v mask ) { return mask.data; }
constexpr friend int popcount( mask_v mask ) { return mask.cast(); }
constexpr friend int select( mask_v mask, int a, int b ) { return mask.cast() ? a : b; }
friend std::ostream& operator<<( std::ostream& os, mask_v mask ) { return os << "scalar{" << mask.data << "}"; }
......@@ -254,10 +256,6 @@ namespace SIMDWrapper {
int data{};
};
constexpr int popcount( mask_v mask ) { return mask.cast(); }
constexpr int select( mask_v mask, int a, int b ) { return mask.cast() ? a : b; }
template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
constexpr T hmax( T data, mask_v mask ) {
return mask.cast() ? data : std::numeric_limits<T>::lowest();
......@@ -276,7 +274,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask ) { return scalar::popcount( mask ); }
static mask_v loop_mask( int offset, int size ) { return size > offset; }
static constexpr InstructionSet instructionSet() { return Scalar; }
};
......@@ -583,6 +580,8 @@ namespace SIMDWrapper {
friend bool testbit( mask_v mask, int bit ) { return ( movemask_u32( mask.native() ) & ( 1 << bit ) ) != 0; }
// TODO: c++20 countl_zero
friend int ffs( mask_v mask ) { return __builtin_ffs( movemask_u32( mask.native() ) ); }
// TODO: c++20 std::popcount
friend int popcount( mask_v mask ) { return __builtin_popcount( movemask_u32( mask.native() ) ); }
friend std::ostream& operator<<( std::ostream& os, mask_v x ) {
return print_bitset<4>( os, static_cast<unsigned long long>( movemask_u32( x.native() ) ), "neon" );
......@@ -594,9 +593,6 @@ namespace SIMDWrapper {
uint32x4_t data{};
};
// TODO: c++20 std::popcount
inline int popcount( mask_v mask ) { return __builtin_popcount( movemask_u32( mask.native() ) ); }
struct types {
constexpr static std::size_t size = 4;
using int_v = neon::int_v;
......@@ -606,7 +602,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask ) { return neon::popcount( mask ); }
static mask_v loop_mask( int i, int n ) { return vcltq_s32( int32x4_t{0, 1, 2, 3}, vdupq_n_s32( n - i ) ); }
static constexpr InstructionSet instructionSet() { return Neon; }
};
......@@ -1014,6 +1009,7 @@ namespace SIMDWrapper {
friend bool testbit( mask_v mask, int bit ) { return ( _mm_movemask_ps( mask.native() ) & ( 1 << bit ) ) != 0; }
// TODO: c++20 countl_zero
friend int ffs( mask_v mask ) { return __builtin_ffs( _mm_movemask_ps( mask.native() ) ); }
friend int popcount( mask_v mask ) { return _mm_popcnt_u32( _mm_movemask_ps( mask.native() ) ); }
friend std::ostream& operator<<( std::ostream& os, mask_v x ) {
return print_bitset<4>( os, static_cast<unsigned long long>( _mm_movemask_ps( x.native() ) ), "sse" );
......@@ -1025,8 +1021,6 @@ namespace SIMDWrapper {
__m128 data{};
};
inline int popcount( mask_v mask ) { return _mm_popcnt_u32( _mm_movemask_ps( mask.native() ) ); }
struct types {
constexpr static std::size_t size = 4;
using int_v = sse::int_v;
......@@ -1036,7 +1030,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask ) { return sse::popcount( mask ); }
static mask_v loop_mask( int i, int n ) {
return _mm_cmplt_ps( _mm_setr_ps( 0, 1, 2, 3 ), _mm_set1_ps( n - i ) );
}
......@@ -1506,6 +1499,7 @@ namespace SIMDWrapper {
}
// TODO: c++20 countl_zero
friend int ffs( mask_v mask ) { return __builtin_ffs( _mm256_movemask_ps( mask.native() ) ); }
friend int popcount( mask_v mask ) { return _mm_popcnt_u32( _mm256_movemask_ps( mask.native() ) ); }
__m256 native() const { return data; }
......@@ -1522,7 +1516,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask );
static mask_v loop_mask( int i, int n );
static constexpr InstructionSet instructionSet() { return AVX2; }
};
......@@ -1650,8 +1643,6 @@ namespace SIMDWrapper {
__m256 data{};
};
inline int popcount( mask_v mask ) { return _mm_popcnt_u32( _mm256_movemask_ps( mask.native() ) ); }
class int_v {
public:
constexpr static std::size_t size() { return avx2::types::size; }
......@@ -1840,8 +1831,6 @@ namespace SIMDWrapper {
inline int_v types::indices( int start ) { return indices() + start; }
inline int types::popcount( mask_v mask ) { return avx2::popcount( mask ); }
inline mask_v types::loop_mask( int i, int n ) {
return _mm256_cmp_ps( _mm256_setr_ps( 0, 1, 2, 3, 4, 5, 6, 7 ), _mm256_set1_ps( n - i ), _CMP_LT_OS );
}
......@@ -1905,6 +1894,7 @@ namespace SIMDWrapper {
friend bool testbit( mask_v mask, int bit ) { return ( mask.native() & ( 1 << bit ) ) != 0; }
// TODO: c++20 countl_zero
friend int ffs( mask_v mask ) { return __builtin_ffs( mask.native() ); }
friend int popcount( mask_v mask ) { return _mm_popcnt_u32( mask.native() ); }
friend std::ostream& operator<<( std::ostream& os, mask_v x ) {
return print_bitset<8>( os, x.native(), "avx256" );
......@@ -1916,8 +1906,6 @@ namespace SIMDWrapper {
__mmask8 data{};
};
inline int popcount( mask_v mask ) { return _mm_popcnt_u32( mask.native() ); }
class float_v;
struct types {
......@@ -1929,7 +1917,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask ) { return avx256::popcount( mask ); }
static mask_v loop_mask( int i, int n ) { return ( ( i + 8 ) > n ) ? ~( 0xFF << ( n - i ) ) : 0xFF; }
static constexpr InstructionSet instructionSet() { return AVX256; }
};
......@@ -2249,6 +2236,7 @@ namespace SIMDWrapper {
friend bool testbit( mask_v mask, int bit ) { return ( mask.native() & ( 1 << bit ) ) != 0; }
// TODO: c++20 countl_zero
friend int ffs( mask_v mask ) { return __builtin_ffs( mask.native() ); }
friend int popcount( mask_v mask ) { return _mm_popcnt_u32( mask.native() ); }
friend std::ostream& operator<<( std::ostream& os, mask_v x ) {
return print_bitset<16>( os, x.native(), "avx512" );
......@@ -2260,8 +2248,6 @@ namespace SIMDWrapper {
__mmask16 data{};
};
inline int popcount( mask_v mask ) { return _mm_popcnt_u32( mask.native() ); }
class float_v;
struct types {
......@@ -2273,7 +2259,6 @@ namespace SIMDWrapper {
static mask_v mask_false() { return mask_v{false}; }
static int_v indices();
static int_v indices( int start );
static int popcount( mask_v mask ) { return avx512::popcount( mask ); }
static mask_v loop_mask( int i, int n ) { return ( ( i + 16 ) > n ) ? ~( 0xFFFF << ( n - i ) ) : 0xFFFF; }
static constexpr InstructionSet instructionSet() { return AVX512; }
};
......
Kernel/LHCbMath/include/LHCbMath/StateVertexUtils.h 0 → 100644
View file @ a89c7663
/***************************************************************************** \
* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#pragma once
#include "GaudiKernel/Point3DTypes.h"
#include "GaudiKernel/Vector3DTypes.h"
#include "LHCbMath/MatVec.h"
namespace ROOT::Math {
template <typename T>
auto X( const T& t ) {
return t.x();
}
template <typename T>
auto Y( const T& t ) {
return t.y();
}
template <typename T>
auto Z( const T& t ) {
return t.z();
}
} // namespace ROOT::Math
namespace LHCb {
namespace StateVertexUtils {
enum ReturnStatus { Failure, Success };
///////////////////////////////////////////////////////////////////////////
/// Return the doca between two track states
///////////////////////////////////////////////////////////////////////////
template <typename StateVector>
auto doca( const StateVector& stateA, const StateVector& stateB );
///////////////////////////////////////////////////////////////////////////
/// Return the distance between a track state and a point
///////////////////////////////////////////////////////////////////////////
template <typename StateVector, typename XYZPoint>
auto doca( const StateVector& stateA, const XYZPoint& pos );
/////////////////////////////////////////////////////////////////////////
/// Compute the chi2 and decaylength of a 'particle' with respect
/// to a vertex. Return 1 if successful.
/// This should probably go into LHCb math.
/////////////////////////////////////////////////////////////////////////
ReturnStatus computeChiSquare( const Gaudi::XYZPoint& pos, const Gaudi::XYZVector& mom,
const Gaudi::SymMatrix6x6& cov6, const Gaudi::XYZPoint& motherpos,
const Gaudi::SymMatrix3x3& mothercov, double& chi2, double& decaylength,
double& decaylengtherr );
/////////////////////////////////////////////////////////////////////////
/// Compute the chi2, decaylength and decaylength uncertainty of a 'particle' with respect
/// to a vertex.
/////////////////////////////////////////////////////////////////////////
template <typename Vec_t, typename Sym6x6_t, typename Sym3x3_t, typename Vec_t1>
auto computeChiSquare( Vec_t const& pos, Vec_t const& mom, Sym6x6_t const& cov6, Vec_t1 const& motherpos,
Sym3x3_t const& mothercov ) {
// pos: decay vertex of particle
// vec: direction or momentum of particle (does not need to be normalized)
// cov6: corresponding covariance matrix
// This calculation is basically a 1-iteration beamspot fit. The
// constraint is
//
// r = x - lambda p/|p| - xbs
//
// where x and p are the position of the decay vertex of the
// candidate and its momentum, lambda is the decaylength and xbs
// the position of the beamspot. The covariance in the constraint
// is
//
// V = Vbs + Vxx - a * Vxp - a Vxp^T + a^2 * Vpp
//
// where a=lambda/|p|^2. It needs an initial estimate for the
// flightlength, for which we simply take the projection of deltaX
// on the direction. We now minimize the chisquare contribution
//
// chi^2 = r^T V^{-1} r
//
// for lambda.
//
// @todo add an implementation for ipchi2-of-composite that just returns
// it without calculating the decay length [error] too. OL tried
// this in https://gitlab.cern.ch/snippets/894, but without a
// better study of which implementation copes better when the
// matrix inversion fails (which it inevitably does sometimes)
// lets not switch to using it.
// const auto dx = pos - motherpos;
const auto dx = Vec_t{pos.X() - motherpos.X(), pos.Y() - motherpos.Y(), pos.Z() - motherpos.Z()};
const auto p3mag = mom.mag();
const auto dir = mom * ( 1.f / p3mag ); // rsqrt
const auto a = dot( dir, dx ) / p3mag;
LHCb::LinAlg::resize_t<Sym6x6_t, 3> W{};
W.fill( [&]( auto i, auto j ) {
return mothercov( i, j ) + cov6( i, j ) + a * a * cov6( i + 3, j + 3 ) -
a * ( cov6( i + 3, j ) + cov6( j + 3, i ) );
} );
// TODO check if easier to just give 3x3 cov matrices instead of creating 6x6 cov matrix
W = W.invChol();
const auto halfdChi2dLam2 = similarity( dir, W );
const auto decaylength = dir.dot( W * dx ) / halfdChi2dLam2;
const auto decaylengtherr = sqrt( 1 / halfdChi2dLam2 ); // TODO: put rsqrt here
const auto res = dx - dir * decaylength;
const auto chi2 = similarity( res, W );
return std::tuple{chi2, decaylength, decaylengtherr};
}
/////////////////////////////////////////////////////////////////////////
/// Compute the point of the doca of two track states. Return 1 if successful.
/////////////////////////////////////////////////////////////////////////
template <typename StateVector, typename XYZPoint>
ReturnStatus poca( const StateVector& stateA, const StateVector& stateB, XYZPoint& vertex );
///////////////////////////////////////////////////////////////////////////
/// Return the doca between two track states
///////////////////////////////////////////////////////////////////////////
template <typename StateVector>
auto doca( const StateVector& stateA, const StateVector& stateB ) {
// first compute the cross product of the directions.
const auto txA = stateA.tx();
const auto tyA = stateA.ty();
const auto txB = stateB.tx();
const auto tyB = stateB.ty();
auto nx = tyA - tyB; // y1 * z2 - y2 * z1
auto ny = txB - txA; // - x1 * z2 + x2 * z1
auto nz = txA * tyB - tyA * txB; // x1 * y2 - x2 * y1
auto n = std::sqrt( nx * nx + ny * ny + nz * nz );
// compute the doca
auto dx = stateA.x() - stateB.x();
auto dy = stateA.y() - stateB.y();
auto dz = stateA.z() - stateB.z();
auto ndoca = dx * nx + dy * ny + dz * nz;
return ndoca / n;
}
///////////////////////////////////////////////////////////////////////////
/// Return the distance between a track state and a point
///////////////////////////////////////////////////////////////////////////
template <typename StateVector, typename XYZPoint>
auto doca( const StateVector& state, const XYZPoint& pos ) {
auto tx = state.tx();
auto ty = state.ty();
auto dz = pos.z() - state.z();
auto dx = state.x() + dz * tx - pos.x();
auto dy = state.y() + dz * ty - pos.y();
return std::sqrt( ( dx * dx + dy * dy ) / ( 1.0 + tx * tx + ty * ty ) );
}
/////////////////////////////////////////////////////////////////////////
/// Compute the point of the doca of two track states. Return 1 if successful.
/////////////////////////////////////////////////////////////////////////
template <typename StateVector, typename XYZPoint>
ReturnStatus poca( const StateVector& stateA, const StateVector& stateB, XYZPoint& vertex ) {
// find the z-positions of the doca. then take the average of the two points.
// see also Gaudi::Math::closestPointParams.
const auto zA = stateA.z();
const auto xA = stateA.x();
const auto yA = stateA.y();
const auto txA = stateA.tx();
const auto tyA = stateA.ty();
const auto zB = stateB.z();
const auto xB = stateB.x();
const auto yB = stateB.y();
const auto txB = stateB.tx();
const auto tyB = stateB.ty();
// define d^2 = ( xA + muA*txA - xB - muB*txB)^2 + ( yA + muA*tyA - xB - muB*tyB)^2 + (zA + muA - zB - muB)^2
// compute (half) 2nd derivative to muA and muB in point muA=muB=0
auto secondAA = txA * txA + tyA * tyA + 1;
auto secondBB = txB * txB + tyB * tyB + 1;
auto secondAB = -txA * txB - tyA * tyB - 1;
// compute inverse matrix, but stop if determinant not positive
auto det = secondAA * secondBB - secondAB * secondAB;
ReturnStatus rc = Success;
if ( std::abs( det ) > 0 ) {
auto secondinvAA = secondBB / det;
auto secondinvBB = secondAA / det;
auto secondinvAB = -secondAB / det;
// compute (half) first derivative
auto firstA = txA * ( xA - xB ) + tyA * ( yA - yB ) + ( zA - zB );
auto firstB = -txB * ( xA - xB ) - tyB * ( yA - yB ) - ( zA - zB );
// compute muA and muB with delta-mu = - seconderivative^-1 * firstderivative
auto muA = -( secondinvAA * firstA + secondinvAB * firstB );
auto muB = -( secondinvBB * firstB + secondinvAB * firstA );
// return the average point
vertex.SetX( 0.5 * ( xA + muA * txA + xB + muB * txB ) );
vertex.SetY( 0.5 * ( yA + muA * tyA + yB + muB * tyB ) );
vertex.SetZ( 0.5 * ( zA + muA + zB + muB ) );
} else {
rc = Failure; // parallel tracks
}
return rc;
}
///////////////////////////////////////////////////////////////////////////
/// Return the chi2 of a track state with respect to a
/// vertex. This is also known as the 'IPCHI2'.
///////////////////////////////////////////////////////////////////////////
template <typename STATE, typename POS, typename POSCOV>
auto vertexChi2( const STATE& state, const POS& vertexpos, const POSCOV& vertexcov ) {
auto tx = state.tx();
auto ty = state.ty();
auto dz = Z( vertexpos ) - state.z();
auto dx = state.x() + dz * tx - X( vertexpos );
auto dy = state.y() + dz * ty - Y( vertexpos );
// compute the covariance matrix. first only the trivial parts:
const auto& statecov = state.covariance();
auto cov00 = vertexcov( 0, 0 ) + statecov( 0, 0 );
auto cov10 = vertexcov( 1, 0 ) + statecov( 1, 0 );
auto cov11 = vertexcov( 1, 1 ) + statecov( 1, 1 );
// add the contribution from the extrapolation
cov00 += dz * dz * statecov( 2, 2 ) + 2 * dz * statecov( 2, 0 );
cov10 += dz * dz * statecov( 3, 2 ) + dz * ( statecov( 3, 0 ) + statecov( 2, 1 ) );
cov11 += dz * dz * statecov( 3, 3 ) + 2 * dz * statecov( 3, 1 );
// add the contribution from pv Z
cov00 += tx * tx * vertexcov( 2, 2 ) - 2 * tx * vertexcov( 2, 0 );
cov10 += tx * ty * vertexcov( 2, 2 ) - ty * vertexcov( 2, 0 ) - tx * vertexcov( 2, 1 );
cov11 += ty * ty * vertexcov( 2, 2 ) - 2 * ty * vertexcov( 2, 1 );
// invert the covariance matrix on the fly
return ( dx * dx * cov11 - 2 * dx * dy * cov10 + dy * dy * cov00 ) / ( cov00 * cov11 - cov10 * cov10 );
}
///////////////////////////////////////////////////////////////////////////
/// Return the chi2 of the vertex of two track states
///////////////////////////////////////////////////////////////////////////
template <typename STATEA, typename STATEB>
auto vertexChi2( const STATEA& stateA, const STATEB& stateB ) {
// first compute the cross product of the directions. we'll need this in any case
const auto txA = stateA.tx();
const auto tyA = stateA.ty();
const auto txB = stateB.tx();
const auto tyB = stateB.ty();
auto nx = tyA - tyB; // y1 * z2 - y2 * z1
auto ny = txB - txA; // - x1 * z2 + x2 * z1
auto nz = txA * tyB - tyA * txB; // x1 * y2 - x2 * y1
// auto n2 = nx*nx + ny*ny + nz*nz ;
// compute doca. we don't divide by the normalization to save time. we call it 'ndoca'
auto dx = stateA.x() - stateB.x();
auto dy = stateA.y() - stateB.y();
auto dz = stateA.z() - stateB.z();
auto ndoca = dx * nx + dy * ny + dz * nz;
// figure out what floating point type to use for the covariance matrix manipulations
using float_t = std::decay_t<decltype( stateA.covariance()( 0, 0 ) )>;
using Vector4 = ROOT::Math::SVector<float_t, 4>;
using SymMatrix4x4 = ROOT::Math::SMatrix<float_t, 4, 4, ROOT::Math::MatRepSym<float_t, 4>>;
// the hard part: compute the jacobians :-)
Vector4 jacA, jacB;
jacA( 0 ) = nx;
jacA( 1 ) = ny;
jacA( 2 ) = -dy + dz * tyB;
jacA( 3 ) = dx - dz * txB;
jacB( 0 ) = -nx;
jacB( 1 ) = -ny;
jacB( 2 ) = dy - dz * tyA;
jacB( 3 ) = -dx + dz * txA;
// compute the variance on ndoca
decltype( auto ) covA = stateA.covariance();
decltype( auto ) covB = stateB.covariance();
using ROOT::Math::Similarity;
auto const varndoca = Similarity( jacA, covA.template Sub<SymMatrix4x4>( 0, 0 ) ) +
Similarity( jacB, covB.template Sub<SymMatrix4x4>( 0, 0 ) );
// return the chi2
return ndoca * ndoca / varndoca;
}
} // namespace StateVertexUtils
} // namespace LHCb
Kernel/LHCbMath/src/StateVertexUtils.cpp 0 → 100644
View file @ a89c7663
/*****************************************************************************\
* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
* *
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#include "LHCbMath/StateVertexUtils.h"
#include "GaudiKernel/GenericVectorTypes.h"
namespace LHCb {
namespace StateVertexUtils {
/////////////////////////////////////////////////////////////////////////
/// Compute the chi2 and decaylength of a 'particle' with respect
/// to a vertex.
/////////////////////////////////////////////////////////////////////////
Gaudi::Vector3 transform( const Gaudi::XYZVector& vec ) { return Gaudi::Vector3( vec.X(), vec.Y(), vec.Z() ); }
ReturnStatus computeChiSquare( const Gaudi::XYZPoint& pos, const Gaudi::XYZVector& vec,
const Gaudi::SymMatrix6x6& cov6, const Gaudi::XYZPoint& motherpos,
const Gaudi::SymMatrix3x3& mothercov, double& chi2, double& decaylength,
double& decaylengtherr ) {
// pos: decay vertex of particle
// vec: direction or momentum of particle (does not need to be normalized)
// cov6: corresponding covariance matrix
// This calculation is basically a 1-iteration beamspot fit. The
// constraint is
//
// r = x - lambda p/|p| - xbs
//
// where x and p are the position of the decay vertex of the
// candidate and its momentum, lambda is the decaylength and xbs
// the position of the beamspot. The covariance in the constraint
// is
//
// V = Vbs + Vxx - a * Vxp - a Vxp^T + a^2 * Vpp
//
// where a=lambda/|p|^2. It needs an initial estimate for the
// flightlength, for which we simply take the projection of deltaX
// on the direction. We now minimize the chisquare contribution
//
// chi^2 = r^T V^{-1} r
//
// for lambda.
Gaudi::Vector3 dx = transform( pos - motherpos );
double p3mag = vec.R();
Gaudi::Vector3 dir = transform( vec.Unit() );
double a = ROOT::Math::Dot( dir, dx ) / p3mag;
Gaudi::SymMatrix3x3 W;
for ( size_t row = 0; row < 3; ++row )
for ( size_t col = 0; col <= row; ++col )
W( row, col ) = mothercov( row, col ) + cov6( row, col ) + a * a * cov6( row + 3, col + 3 ) -
a * ( cov6( row + 3, col ) + cov6( col + 3, row ) );
int OK = W.Invert();
double halfdChi2dLam2 = ROOT::Math::Similarity( W, dir );
decaylength = ROOT::Math::Dot( dir, W * dx ) / halfdChi2dLam2;
decaylengtherr = std::sqrt( 1 / halfdChi2dLam2 );
Gaudi::Vector3 res = dx - decaylength * dir;
chi2 = ROOT::Math::Similarity( W, res );
return OK ? Success : Failure;
}
} // namespace StateVertexUtils
} // namespace LHCb