diff --git a/Phys/FunctorCore/CMakeLists.txt b/Phys/FunctorCore/CMakeLists.txt index 5895820f35b84593f24d0de39bddcf02a3a14062..c3b1131ae520fddaa3e5a7f9de8af8d74da470e6 100644 --- a/Phys/FunctorCore/CMakeLists.txt +++ b/Phys/FunctorCore/CMakeLists.txt @@ -240,7 +240,7 @@ string(STRIP ${CMAKE_CXX_COMPILER_CONTENT} CMAKE_CXX_COMPILER_CONTENT) # Specify the libraries a JIT compiled functor library will link against. The # list here is defined based upon the includes present in # `FunctorCore/include/JIT_includes.h` -set(JIT_LINK_LIBS "-lFunctorCoreLib -lTrackEvent -lPhysEvent -lMCEvent -lRecEvent -lHltEvent -lTMVA -lGaudiKernel") +set(JIT_LINK_LIBS "-lFunctorCoreLib -lTrackEvent -lPhysEvent -lMCEvent -lRecEvent -lHltEvent -lTMVA -lTrackKernel -lSelToolsLib -lGaudiKernel") file(GENERATE OUTPUT "functor_jitter_tmp" diff --git a/Phys/FunctorCore/include/Functors/Combination.h b/Phys/FunctorCore/include/Functors/Combination.h index b232a1a8d8181e8ac367dbb4ca7645ef4e37aad8..ac36ac48cab6b0f2108b44a108b6607e62ca2264 100644 --- a/Phys/FunctorCore/include/Functors/Combination.h +++ b/Phys/FunctorCore/include/Functors/Combination.h @@ -19,8 +19,12 @@ #include "GaudiKernel/GaudiException.h" #include "GaudiKernel/System.h" #include "Kernel/IDistanceCalculator.h" +#include "Kernel/ITrajPoca.h" #include "Kernel/STLExtensions.h" +#include "LHCbMath/GoldenSectionMinimizer.h" #include "SelTools/DistanceCalculator.h" +#include "TrackInterfaces/ITrackStateProvider.h" +#include "TrackKernel/TrackTraj.h" /** @file Combination.h * @brief Definitions of functors for tuple-likes of track-like objects. @@ -363,6 +367,116 @@ namespace Functors::detail { float m_thresh; }; + class INonlinearPOCAProvider { + + class Geometry { + const void* m_geom = nullptr; + const LHCb::DetDesc::ConditionContext& ( *m_cond )( void const*, EventContext const& ) = nullptr; + MsgStream& ( *m_get_warning )( void const* ) = nullptr; + LHCb::DetDesc::ConditionAccessor m_det; + + public: + template + Geometry( T* ptr, const std::string& geom_name ) + : m_geom{ ptr } + , m_cond{ []( const void* geom, const EventContext& ctx ) -> decltype( auto ) { + return static_cast( geom )->getConditionContext( ctx ); + } } + , m_get_warning{ []( const void* geom ) -> decltype( auto ) { + return static_cast( geom )->warning(); + } } + , m_det( ptr, geom_name ) {} + + MsgStream& warning() const { return ( *m_get_warning )( m_geom ); } + + auto* geometry( EventContext const& ctx ) const { return m_det.get( ( *m_cond )( m_geom, ctx ) ).geometry(); } + }; + + std::string m_tool_name, m_trajpoca_name, m_geom_name; + bool m_allow_default_tool_creation = false; + std::optional m_det; + std::optional> m_stateProvider; + std::optional> m_trajpoca; + + public: + INonlinearPOCAProvider( std::string tool_name, std::string trajpoca_name, bool allow_default_tool_creation = false, + std::string geom_name = LHCb::standard_geometry_top ) + : m_tool_name{ std::move( tool_name ) } + , m_trajpoca_name{ std::move( trajpoca_name ) } + , m_geom_name{ std::move( geom_name ) } + , m_allow_default_tool_creation{ allow_default_tool_creation } {} + + INonlinearPOCAProvider( INonlinearPOCAProvider&& rhs ) + : m_tool_name{ std::move( rhs.m_tool_name ) } + , m_trajpoca_name{ std::move( rhs.m_trajpoca_name ) } + , m_geom_name{ std::move( rhs.m_geom_name ) } + , m_allow_default_tool_creation{ rhs.m_allow_default_tool_creation } {} + + INonlinearPOCAProvider& operator=( INonlinearPOCAProvider&& rhs ) { + m_tool_name = std::move( rhs.m_tool_name ); + m_trajpoca_name = std::move( rhs.m_trajpoca_name ); + m_geom_name = std::move( rhs.m_geom_name ); + m_allow_default_tool_creation = rhs.m_allow_default_tool_creation; + return *this; + } + + void emplace( TopLevelInfo& top_level ) { + m_stateProvider.emplace( m_tool_name, m_allow_default_tool_creation ); + m_trajpoca.emplace( m_trajpoca_name, m_allow_default_tool_creation ); + + constexpr auto check_tool = []( const StatusCode code ) { + if ( code.isFailure() ) { + throw GaudiException{ "INonlinearPOCAProvider::bind", + "failed to retrieve tools. Please, ensure the tools are properly instantiated (for " + "example with ForceToolInstantiation algorithm) and check the spelling of the tool " + "names! This exception might be disabled by setting allow_default_tool_creation=True", + StatusCode::FAILURE }; + } + }; + + auto try_emplace = [&]( auto* det ) { + if ( det ) { + m_det.emplace( det, m_geom_name ); + check_tool( m_stateProvider->retrieve() ); + check_tool( m_trajpoca->retrieve() ); + } + return det; + }; + + if ( try_emplace( dynamic_cast*>( top_level.algorithm() ) ) ) return; + if ( try_emplace( dynamic_cast*>( top_level.algorithm() ) ) ) + return; + + throw GaudiException{ "INonlinearPOCAProvider::bind", + "attempting to bind an alorithm which is not a condition holder", StatusCode::FAILURE }; + } + auto prepare( EventContext const& ctx ) const { + assert( m_stateProvider.has_value() ); + assert( m_det.has_value() ); + if ( !m_stateProvider->get() ) { m_det.value().warning() << "retrieval of extrapolator too late!!!" << endmsg; } + return [trajpoca = m_trajpoca->get(), geom = m_det->geometry( ctx ), det = LHCb::get_pointer( m_det ), + stateProvider = m_stateProvider->get()]( auto const& trA, auto const& trB ) -> std::optional { + auto const* trackA = &Sel::Utils::deref_if_ptr( trA ); + auto const* trackB = &Sel::Utils::deref_if_ptr( trB ); + + if constexpr ( requires { + Sel::NonlinearDistanceCalculator::findPOCA( + *( stateProvider->trajectory( *( trackA ), *geom ) ), + *( stateProvider->trajectory( *( trackB ), *geom ) ), *trajpoca ); + } ) { + return Sel::NonlinearDistanceCalculator::findPOCA( *( stateProvider->trajectory( *( trackA ), *geom ) ), + *( stateProvider->trajectory( *( trackB ), *geom ) ), + *trajpoca ); + } else { + det->warning() << " requested to use INonlinearPOCAProvider, but that does not support the arguments -- " + << System::typeinfoName( typeid( trackA ) ) << " " << System::typeinfoName( typeid( trackB ) ) + << endmsg; + return std::nullopt; + } + }; + } + }; + } // namespace Functors::detail /** @namespace Functors::Combination @@ -440,4 +554,49 @@ namespace Functors::Combination { using LHCb::Event::charge; return charge( item ); } }; + + /** @brief Retrive the distance of closest approach between two tracks using the TrackStateProvider tool for + * extrapolation + */ + struct NonlinearDoca : public Function { + NonlinearDoca( std::string stateprovider_name = "TrackStateProvider", std::string trajpoca_name = "TrajPoca", + bool allow_default_tool_creation = false ) + : m_extrapolator{ std::move( stateprovider_name ), std::move( trajpoca_name ), allow_default_tool_creation } {} + + void bind( TopLevelInfo& top_level ) { m_extrapolator.emplace( top_level ); } + + auto prepare( EventContext const& evtCtx, TopLevelInfo const& ) const { + return [extrapolate = m_extrapolator.prepare( evtCtx )]( const auto& trA, + const auto& trB ) -> Functors::Optional { + if ( const auto res = extrapolate( trA, trB ); res.has_value() ) { return res->doca; } + return std::nullopt; + }; + } + + private: + detail::INonlinearPOCAProvider m_extrapolator; + }; + + /** @brief Retrive the point of closest approach between two tracks using the TrackStateProvider tool for + * extrapolation + */ + struct NonlinearPoca : public Function { + NonlinearPoca( std::string stateprovider_name = "TrackStateProvider", std::string trajpoca_name = "TrajPoca", + bool allow_default_tool_creation = false ) + : m_extrapolator{ std::move( stateprovider_name ), std::move( trajpoca_name ), allow_default_tool_creation } {} + + void bind( TopLevelInfo& top_level ) { m_extrapolator.emplace( top_level ); } + + auto prepare( EventContext const& evtCtx, TopLevelInfo const& ) const { + return [extrapolate = m_extrapolator.prepare( evtCtx )]( + const auto& trA, const auto& trB ) -> Functors::Optional { + if ( const auto res = extrapolate( trA, trB ); res.has_value() ) { return res->poca; } + return std::nullopt; + }; + }; + + private: + detail::INonlinearPOCAProvider m_extrapolator; + }; + } // namespace Functors::Combination diff --git a/Phys/FunctorCore/include/Functors/TrackLike.h b/Phys/FunctorCore/include/Functors/TrackLike.h index 91246851747ea3d0ebd6409630755a58cbe252df..e751dbfdbbc333e00517601d2754072274c8f1b0 100644 --- a/Phys/FunctorCore/include/Functors/TrackLike.h +++ b/Phys/FunctorCore/include/Functors/TrackLike.h @@ -29,10 +29,12 @@ #include "GaudiKernel/System.h" #include "GaudiKernel/Vector4DTypes.h" #include "Kernel/HitPattern.h" +#include "LHCbMath/GoldenSectionMinimizer.h" #include "MCInterfaces/IMCReconstructed.h" #include "SelKernel/Utilities.h" #include "SelKernel/VertexRelation.h" #include "TrackInterfaces/ITrackExtrapolator.h" +#include "TrackInterfaces/ITrackStateProvider.h" #include #include @@ -157,6 +159,119 @@ namespace Functors::Track::detail { } }; + /** + * @brief helper for ExtrapolateToPoca functor + */ + class ITrackExtrapolatorToPocaHolder { + + class Geometry { + const void* m_geom = nullptr; + const LHCb::DetDesc::ConditionContext& ( *m_cond )( void const*, EventContext const& ) = nullptr; + MsgStream& ( *m_get_warning )( void const* ) = nullptr; + LHCb::DetDesc::ConditionAccessor m_det; + + public: + template + Geometry( T* ptr, const std::string& geom_name ) + : m_geom{ ptr } + , m_cond{ []( const void* geom, const EventContext& ctx ) -> decltype( auto ) { + return static_cast( geom )->getConditionContext( ctx ); + } } + , m_get_warning{ []( const void* geom ) -> decltype( auto ) { + return static_cast( geom )->warning(); + } } + , m_det( ptr, geom_name ) {} + + MsgStream& warning() const { return ( *m_get_warning )( m_geom ); } + + auto* geometry( EventContext const& ctx ) const { return m_det.get( ( *m_cond )( m_geom, ctx ) ).geometry(); } + }; + + typedef LHCb::Math::GoldenSectionMinimizer GoldenSectionMinimizer; + + Gaudi::XYZPointF m_pos; + GoldenSectionMinimizer m_minimizer; + std::optional m_det; + std::string m_tool_name, m_geom_name; + bool m_allow_default_tool_creation = false; + std::optional> m_stateProvider; + + public: + ITrackExtrapolatorToPocaHolder( Gaudi::XYZPointF pos, std::string tool_name, float zmin, float zmax, + float tolerance, unsigned maxiter, bool allow_default_tool_creation = false, + std::string geom_name = LHCb::standard_geometry_top ) + : m_pos( pos ) + , m_minimizer{ zmin, zmax, tolerance, maxiter } + , m_tool_name{ std::move( tool_name ) } + , m_geom_name{ std::move( geom_name ) } + , m_allow_default_tool_creation{ allow_default_tool_creation } {} + + ITrackExtrapolatorToPocaHolder( ITrackExtrapolatorToPocaHolder&& rhs ) + : m_pos( rhs.m_pos ) + , m_minimizer{ std::move( rhs.m_minimizer ) } + , m_tool_name{ std::move( rhs.m_tool_name ) } + , m_geom_name{ std::move( rhs.m_geom_name ) } + , m_allow_default_tool_creation{ rhs.m_allow_default_tool_creation } {} + + ITrackExtrapolatorToPocaHolder& operator=( ITrackExtrapolatorToPocaHolder&& rhs ) { + m_pos = std::move( rhs.m_pos ); + m_minimizer = std::move( rhs.m_minimizer ); + m_tool_name = std::move( rhs.m_tool_name ); + m_geom_name = std::move( rhs.m_geom_name ); + m_allow_default_tool_creation = rhs.m_allow_default_tool_creation; + return *this; + } + + void emplace( TopLevelInfo& top_level ) { + m_stateProvider.emplace( m_tool_name, m_allow_default_tool_creation ); + + auto try_emplace = [&]( auto* det ) { + if ( det ) { + m_det.emplace( det, m_geom_name ); + if ( m_stateProvider->retrieve().isFailure() ) { + throw GaudiException{ "TrackExtrapolatorToPoca::bind", + "failed to retrieve tools. Please, ensure the tools are properly instantiated (for " + "example with ForceToolInstantiation algorithm) and check the spelling of the tool " + "names! This exception might be disabled by setting allow_default_tool_creation=True", + StatusCode::FAILURE }; + } + } + return det; + }; + + if ( try_emplace( dynamic_cast*>( top_level.algorithm() ) ) ) return; + if ( try_emplace( dynamic_cast*>( top_level.algorithm() ) ) ) + return; + + throw GaudiException{ "TrackExtrapolatorToPoca::bind", + "attempting to bind an alorithm which is not a condition holder", StatusCode::FAILURE }; + } + auto prepare( EventContext const& ctx ) const { + assert( m_stateProvider.has_value() ); + assert( m_det.has_value() ); + if ( !m_stateProvider->get() ) { m_det.value().warning() << "retrieval of extrapolator too late!!!" << endmsg; } + return [geom = m_det->geometry( ctx ), det = LHCb::get_pointer( m_det ), pos = m_pos, minimizer = m_minimizer, + stateProvider = m_stateProvider->get()]( auto const& tr ) -> Functors::Optional { + auto const* track = &Sel::Utils::deref_if_ptr( tr ); + if constexpr ( requires { + Sel::GSDistanceCalculator::findClosestStateToPoint( + pos, *stateProvider->trajectory( *track, *geom ), minimizer ); + } ) { + const auto res = Sel::GSDistanceCalculator::findClosestStateToPoint( + pos, *stateProvider->trajectory( *track, *geom ), minimizer ); + if ( res ) + return res.value(); + else + return std::nullopt; + } else { + det->warning() << " requested to use ITrackExtrapolator, but that does not support the argument -- " + << System::typeinfoName( typeid( decltype( track ) ) ) << " " << endmsg; + return std::nullopt; + } + }; + } + }; + } // namespace Functors::Track::detail namespace Functors::Track { @@ -175,6 +290,23 @@ namespace Functors::Track { detail::ITrackExtrapolatorHolder m_extrapolator; }; + /** @brief Retrieve state of track extrapolated to the POCA to a given position + */ + struct ExtrapolateToPoca : public Function { + ExtrapolateToPoca( float x, float y, float z, std::string tool_name = "TrackStateProvider", float zmin = 1.f, + float zmax = 10000.f, float tolerance = 0.01f, unsigned maxiter = 50, + bool allow_default_tool_creation = false ) + : m_extrapolator{ Gaudi::XYZPointF( x, y, z ), std::move( tool_name ), zmin, zmax, tolerance, maxiter, + allow_default_tool_creation } {} + + void bind( TopLevelInfo& top_level ) { m_extrapolator.emplace( top_level ); } + + auto prepare( EventContext const& evtCtx, TopLevelInfo const& ) const { return m_extrapolator.prepare( evtCtx ); } + + private: + detail::ITrackExtrapolatorToPocaHolder m_extrapolator; + }; + /** * @brief Retrieve const Container with pointers to all the states */ diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py index f00a1a415faabe679fe11af2a6187bb26916e7a2..c0d11e76cd1da8befae38ca00748702de6999869 100644 --- a/Phys/FunctorCore/python/Functors/__init__.py +++ b/Phys/FunctorCore/python/Functors/__init__.py @@ -413,6 +413,100 @@ EXTRAPOLATE_TRACK = Functor( Params=[("z", "z coordinate to which the track is extrapolated", float)], ) + +def EXTRAPOLATE_TRACK_TO_POCA( + x: float, + y: float, + z: float, + tool_name: str = "TrackStateProvider", + zmin: float = 1.0, + zmax: float = 10000.0, + tolerance: float = 0.01, + maxiter: int = 50, + allow_default_tool_creation: bool = False, +): + """Get the track state at the closest position to the given point (x,y,z) using the used-provided TrackStateProvider. + This algorithm uses golden-section minimisation algorithm. + """ + return EXTRAPOLATE_TRACK_TO_POCA._F( + x=x, + y=y, + z=z, + tool_name=tool_name, + zmin=zmin, + zmax=zmax, + tolerance=tolerance, + maxiter=maxiter, + allow_default_tool_creation=allow_default_tool_creation, + ) + + +EXTRAPOLATE_TRACK_TO_POCA._F = Functor( + "EXTRAPOLATE_TRACK_TO_POCA", + "Track::ExtrapolateToPoca", + "Get the track state at the closest position to the given point (x,y,z)", + Params=[ + ("x", "x coordinate of target point", float), + ("y", "y coordinate of target point", float), + ("z", "z coordinate of target point", float), + ("tool_name", "TrackStateProvider tool to be used", str), + ("zmin", "Minimum z of the search range", float), + ("zmax", "Maximum z of the search range", float), + ("tolerance", "Tolerance of the minimisation", float), + ("maxiter", "Maximum number of iterations", int), + ( + "allow_default_tool_creation", + "Allow creation of default tool if not found in ToolSvc", + bool, + ), + ], +) + + +def GSIP( + pvs, + tool_name="TrackStateProvider", + zmin=-1000.0, + zmax=10000.0, + tolerance=0.01, + maxiter=50, + allow_default_tool_creation=False, +): + """ + Calculates the impact parameter (IP) of a track with respect to a best primary vertex (see BPV definition). + The best primary vertex is choosen using a closest track state to point (0,0,0). + + Args: + pvs: The event primary vertices. + tool_name (str, optional): The name of the track state provider tool to use for extrapolation. + zmin (float, optional): Minimum z of the search range. + zmax (float, optional): Maximum z of the search range. + tolerance (float, optional): Tolerance of the minimisation. + maxiter (int, optional): Maximum number of iterations. + allow_default_tool_creation (bool, optional): Allow the creation of default tools if the provided ones are not found. Defaults to False. + + CAUTION: + If allow_default_tool_creation is set to False (default), please, make sure a tool with the provided name exist and is used by any other GaudiAlgorithm that is executed BEFORE the one using this functor. + + One should be very careful if allow_default_tool_creation is set to True, since if either + - a tool with provided name doesn't exist + - it wasn't included in the public_tools list of the top configure call + the DEFAULT version of the corresponding tool will be created instead! + """ + closest_state = EXTRAPOLATE_TRACK_TO_POCA( + x=0.0, + y=0.0, + z=0.0, + tool_name=tool_name, + zmin=zmin, + zmax=zmax, + tolerance=tolerance, + maxiter=maxiter, + allow_default_tool_creation=allow_default_tool_creation, + ) + return IP.bind(TOLINALG @ POSITION @ BPV(pvs) @ closest_state, closest_state) + + # PHI already has the meaning of phi(slopes) so for now we use PHI_COORDINATE PHI_COORDINATE = Functor( "PHI_COORDINATE", @@ -2043,6 +2137,105 @@ DOCA = __create_doca_functor( ) +def NONLINEAR_DOCA( + child1: int, + child2: int, + stateprovider_name: str = "TrackStateProvider", + trajpoca_name: str = "TrajPoca", + allow_default_tool_creation: bool = False, +): + """Compute the distance of closest approach between two track-like objects using the TrajPoca algorithm with a user-defined TrackStateProvider tool. + The starting point for the minimisation z_init is computed as y-z interception. In case, if the computed z_init is outside of the range [StateParameters::ZEndVelo, StateParameters::ZEndRich2], the StateParameters::ZMidMag is used instead. + + Args: + stateprovider_name (str, optional): TrackStateProvider tool to be used. Defaults to "TrackStateProvider". + trajpoca_name (str, optional): ITrajPoca tool to be used. Defaults to "TrajPoca". + allow_default_tool_creation (bool, optional): Allow the creation of default tools if the provided ones are not found. Defaults to False. + + CAUTION: + If allow_default_tool_creation is set to False (default), please, make sure a tool with the provided name exist and is used by any other GaudiAlgorithm that is executed BEFORE the one using this functor. + + One should be very careful if allow_default_tool_creation is set to True, since if either + - a tool with provided name doesn't exist + - it wasn't included in the public_tools list of the top configure call + the DEFAULT version of the corresponding tool will be created instead! + """ + return NONLINEAR_DOCA._F( + stateprovider_name=stateprovider_name, + trajpoca_name=trajpoca_name, + allow_default_tool_creation=allow_default_tool_creation, + ).bind( + TRACK @ CHILD(child1), + TRACK @ CHILD(child2), + ) + + +NONLINEAR_DOCA._F = Functor( + "_NONLINEAR_DOCA", + "Combination::NonlinearDoca", + """Compute the distance of closest approach between two track-like objects using the TrajPoca algorithm with a user-defined TrackStateProvider tool. + """, + Params=[ + ("stateprovider_name", "TrackStateProvider tool to be used", str), + ("trajpoca_name", "ITrajPoca tool to be used", str), + ( + "allow_default_tool_creation", + "Allow the creation of default tools if the provided ones are not found", + bool, + ), + ], +) + + +def NONLINEAR_POCA( + child1: int, + child2: int, + stateprovider_name: str = "TrackStateProvider", + trajpoca_name: str = "TrajPoca", + allow_default_tool_creation: bool = False, +): + """Compute the position of closest approach between two track-like objects using the TrajPoca with a user-defined TrackStateProvider tool. + The starting point for the minimisation z_init is computed as y-z interception. In case, if the computed z_init is outside of the range [StateParameters::ZEndVelo, StateParameters::ZEndRich2], the StateParameters::ZMidMag is used instead. + + Args: + stateprovider_name (str, optional): TrackStateProvider tool to be used. Defaults to "TrackStateProvider". + trajpoca_name (str, optional): ITrajPoca tool to be used. Defaults to "TrajPoca". + allow_default_tool_creation (bool, optional): Allow the creation of default tools if the provided ones are not found. Defaults to False. + + CAUTION: + If allow_default_tool_creation is set to False (default), please, make sure a tool with the provided name exist and is used by any other GaudiAlgorithm that is executed BEFORE the one using this functor. + + One should be very careful if allow_default_tool_creation is set to True, since if either + - a tool with provided name doesn't exist + - it wasn't included in the public_tools list of the top configure call + the DEFAULT version of the corresponding tool will be created instead! + """ + return NONLINEAR_POCA._F( + stateprovider_name=stateprovider_name, + trajpoca_name=trajpoca_name, + allow_default_tool_creation=allow_default_tool_creation, + ).bind( + TRACK @ CHILD(child1), + TRACK @ CHILD(child2), + ) + + +NONLINEAR_POCA._F = Functor( + "_NONLINEAR_POCA", + "Combination::NonlinearPoca", + """Compute the position of closest approach between two track-like objects using a user-defined TrackStateProvider tool.""", + Params=[ + ("stateprovider_name", "TrackStateProvider tool to be used", str), + ("trajpoca_name", "ITrajPoca tool to be used", str), + ( + "allow_default_tool_creation", + "Allow the creation of default tools if the provided ones are not found", + bool, + ), + ], +) + + SDOCACHI2 = __create_doca_functor( "SDOCACHI2", "Combination::SDistanceOfClosestApproachChi2", diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp index e91d2c6e38313a03a70418161f5126cb2cb787d1..e05b9bd4f5d70828679eee668a4dc2b3f95eeaa6 100644 --- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp +++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp @@ -22,6 +22,7 @@ namespace FPID = Functors::PID; namespace FC = Functors::Common; namespace FF = Functors::Functional; namespace FP = Functors::Particle; +namespace FCMB = Functors::Combination; using mask_arg_t = Functors::mask_arg_t; // @@ -436,6 +437,31 @@ using test_extrapolate_track_result_type2 = std::invoke_result_t; using test_extrapolate_track_result_type4 = std::invoke_result_t; +auto test_extrapolate_to_poca_track = FT::ExtrapolateToPoca( 0, 0, 0 ); +using test_extrapolate_to_poca_track_type = + decltype( test_extrapolate_to_poca_track.prepare( EventContext{}, Functors::TopLevelInfo{} ) ); +using test_extrapolate_to_poca_track_result_type = std::invoke_result_t; +using test_extrapolate_to_poca_track_result_type2 = + std::invoke_result_t; +using test_extrapolate_to_poca_track_result_type3 = + std::invoke_result_t; +using test_extrapolate_to_poca_track_result_type4 = + std::invoke_result_t; + +auto rkdoca = FCMB::NonlinearDoca(); +using rkdoca_type = decltype( rkdoca.prepare( EventContext{}, Functors::TopLevelInfo{} ) ); +using rkdoca_result_type = std::invoke_result_t; +using rkdoca_result_type2 = std::invoke_result_t; +using rkdoca_result_type3 = std::invoke_result_t; +using rkdoca_result_type4 = std::invoke_result_t; + +auto rkpoca = FCMB::NonlinearPoca(); +using rkpoca_type = decltype( rkpoca.prepare( EventContext{}, Functors::TopLevelInfo{} ) ); +using rkpoca_result_type = std::invoke_result_t; +using rkpoca_result_type2 = std::invoke_result_t; +using rkpoca_result_type3 = std::invoke_result_t; +using rkpoca_result_type4 = std::invoke_result_t; + BOOST_AUTO_TEST_CASE( instantiate_all ) { // empty on purpose -- just need to compile and instantiate... } diff --git a/Phys/SelTools/CMakeLists.txt b/Phys/SelTools/CMakeLists.txt index 9789f6c5ea957daade0b0aad785822a782e03c04..c127e59d24c1d88ecd37e47814df4fa720581961 100644 --- a/Phys/SelTools/CMakeLists.txt +++ b/Phys/SelTools/CMakeLists.txt @@ -16,6 +16,7 @@ Phys/SelTools gaudi_add_library(SelToolsLib SOURCES src/Utilities.cpp + src/DistanceCalculator.cpp LINK Gaudi::GaudiKernel LHCb::EventBase @@ -24,4 +25,5 @@ gaudi_add_library(SelToolsLib Rec::SelKernelLib BLAS::BLAS # needed because ROOT::TMVA does not link correctly ROOT::TMVA + Rec::TrackInterfacesLib ) diff --git a/Phys/SelTools/include/SelTools/DistanceCalculator.h b/Phys/SelTools/include/SelTools/DistanceCalculator.h index a1b0cdbce9c8edfc31a13e46477c748a6b551765..fe341e2db47a8a628d368393b7f5003d8980928c 100644 --- a/Phys/SelTools/include/SelTools/DistanceCalculator.h +++ b/Phys/SelTools/include/SelTools/DistanceCalculator.h @@ -17,6 +17,8 @@ #include "LHCbMath/MatrixUtils.h" #include "LHCbMath/SIMDWrapper.h" +#include "TrackInterfaces/ITrackStateProvider.h" + #include "Gaudi/Accumulators.h" #include "Gaudi/Algorithm.h" #include "GaudiKernel/GenericMatrixTypes.h" @@ -27,6 +29,18 @@ #include "GaudiKernel/Vector3DTypes.h" #include "GaudiKernel/Vector4DTypes.h" +#include "Kernel/ITrajPoca.h" +#include "TrackKernel/TrackTraj.h" + +namespace LHCb { + namespace Math { + + template + class GoldenSectionMinimizer; + + } +} // namespace LHCb + namespace { /** Extract a new_dim x new_dim symmetric matrix from another symmetric @@ -175,6 +189,24 @@ namespace Sel { } }; + struct POCA { + Gaudi::XYZPoint poca; + double doca; + }; + + namespace GSDistanceCalculator { + std::optional findClosestStateToPoint( const Gaudi::XYZPointF& pos, const LHCb::TrackTraj& track, + const LHCb::Math::GoldenSectionMinimizer& minimizer ); + }; + + namespace NonlinearDistanceCalculator { + + typedef std::array TrackStates; + + std::optional findPOCA( const LHCb::Trajectory& trajA, const LHCb::Trajectory& trajB, + const ITrajPoca& trajpoca ); + }; // namespace NonlinearDistanceCalculator + namespace helper { // FIXME: we should _not_ need this version... diff --git a/Phys/SelTools/src/DistanceCalculator.cpp b/Phys/SelTools/src/DistanceCalculator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..083df2b116846d61291aa98fb3211f70a841e30e --- /dev/null +++ b/Phys/SelTools/src/DistanceCalculator.cpp @@ -0,0 +1,59 @@ +/*****************************************************************************\ +* (c) Copyright 2025 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 "SelTools/DistanceCalculator.h" +#include "Event/StateParameters.h" +#include "LHCbMath/GoldenSectionMinimizer.h" + +namespace Sel { + + std::optional NonlinearDistanceCalculator::findPOCA( const LHCb::Trajectory& trajA, + const LHCb::Trajectory& trajB, + const ITrajPoca& trajpoca ) { + + Gaudi::XYZVector deltaX; + + const float tyA = trajA.direction( 0 ).Y() / trajA.direction( 0 ).Z(); + const float tyB = trajB.direction( 0 ).Y() / trajB.direction( 0 ).Z(); + const float yAt0A = trajA.position( 0 ).Y() - tyA * trajA.position( 0 ).Z(); + const float yAt0B = trajB.position( 0 ).Y() - tyB * trajB.position( 0 ).Z(); + const float dy = yAt0A - yAt0B; + const float dty = ( tyA - tyB ); + + // First and the last defined states + constexpr auto outside_detector = []( const float z ) { + constexpr auto minZ = -StateParameters::ZEndVelo; + constexpr auto maxZ = StateParameters::ZEndRich2; + return ( z < minZ ) || ( z > maxZ ); + }; + + const auto startZ = + ( LHCb::essentiallyZero( dty ) || outside_detector( -dy / dty ) ) ? StateParameters::ZMidMag : -dy / dty; + + double muA( startZ ), muB( startZ ); + if ( trajpoca.minimize( trajA, muA, trajB, muB, deltaX, 0.1 * Gaudi::Units::mm ).isFailure() ) return std::nullopt; + + const auto pos0 = trajA.position( muA ); + const auto pos1 = trajB.position( muB ); + const auto cntr = pos0 + 0.5 * ( pos1 - pos0 ); + return POCA{ .poca = cntr, .doca = deltaX.R() }; + } + + std::optional + GSDistanceCalculator::findClosestStateToPoint( const Gaudi::XYZPointF& pos, const LHCb::TrackTraj& track, + const LHCb::Math::GoldenSectionMinimizer& minimizer ) { + + const auto result = + minimizer( [pos, &track]( const float zs ) { return ( track.state( zs ).position() - pos ).Mag2(); } ); + if ( !result ) return std::nullopt; + + return track.state( result.value ); + } +} // namespace Sel diff --git a/Pr/PrAlgorithms/CMakeLists.txt b/Pr/PrAlgorithms/CMakeLists.txt index de3251475dfc136c0e203df18371e21ba8dc9b93..9bac6b5864909b57b6132d40608a1763a5bb1369 100644 --- a/Pr/PrAlgorithms/CMakeLists.txt +++ b/Pr/PrAlgorithms/CMakeLists.txt @@ -28,6 +28,7 @@ gaudi_add_module(PrAlgorithms src/PrStoreSciFiHits.cpp src/PrStoreUTHit.cpp src/HitsEmptyProducer.cpp + src/ForceToolInstantiation.cpp LINK Boost::container Boost::headers diff --git a/Pr/PrAlgorithms/src/ForceToolInstantiation.cpp b/Pr/PrAlgorithms/src/ForceToolInstantiation.cpp new file mode 100644 index 0000000000000000000000000000000000000000..73ecb104bc16d50a88853eaf20a97eb19a21b9ea --- /dev/null +++ b/Pr/PrAlgorithms/src/ForceToolInstantiation.cpp @@ -0,0 +1,37 @@ +/***********************************************************************************\ +* (c) Copyright 1998-2019 CERN for the benefit of the LHCb and ATLAS collaborations * +* * +* This software is distributed under the terms of the Apache version 2 licence, * +* copied verbatim in the file "LICENSE". * +* * +* 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 "LHCbAlgs/Consumer.h" + +/** + * @brief Forces instantiation of the provided tool + * + * This template class is designed to handle cases where certain algorithms (such as FunTuple) + * cannot directly create required tools due to framework limitations. + * By creating a separate algorithm that forces instantiation of the tool, the latter can be referred to + * by name. + * + * The current algorithm has to be added to the control flow for correct operation. + * + * @param Tool The ToolHandle of the required tool + */ + +class ForceToolInstantiation final : public LHCb::Algorithm::Consumer { +public: + ForceToolInstantiation( const std::string& name, ISvcLocator* pSvcLocator ) : Consumer( name, pSvcLocator, {} ){}; + + void operator()() const override{}; + +private: + PublicToolHandle m_tool{ this, "Tool", "" }; +}; + +DECLARE_COMPONENT_WITH_ID( ForceToolInstantiation, "ForceToolInstantiation" )