From 2059d1589e123ce57174cc5491dfe4da5e7f8255 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Wed, 27 Apr 2022 23:08:58 +0200 Subject: [PATCH 01/14] first working hit search algo --- Pr/PrPixel/CMakeLists.txt | 2 + Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 192 ++++++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp diff --git a/Pr/PrPixel/CMakeLists.txt b/Pr/PrPixel/CMakeLists.txt index 7f56490f906..ce0a95d469b 100644 --- a/Pr/PrPixel/CMakeLists.txt +++ b/Pr/PrPixel/CMakeLists.txt @@ -20,6 +20,7 @@ gaudi_add_module(PrPixel src/VeloClusterTrackingSIMD.cpp src/VeloKalman.cpp src/PrVPHitsToVPLightClusters.cpp + src/VeloHeavyFlavourTracking.cpp LINK Gaudi::GaudiAlgLib LHCb::DAQEventLib @@ -32,4 +33,5 @@ gaudi_add_module(PrPixel LHCb::VPKernel LHCb::VPDAQLib Rec::PrKernel + Rec::SelKernelLib ) diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp new file mode 100644 index 00000000000..eb5f19e3886 --- /dev/null +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -0,0 +1,192 @@ +/*****************************************************************************\ +* (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 "DetDesc/GenericConditionAccessorHolder.h" +#include "DetDesc/IConditionDerivationMgr.h" +#include "Event/Particle.h" +#include "Event/PrHits.h" +#include "Event/PrimaryVertices.h" +#include "Event/RawEvent.h" +#include "Event/RelationTable.h" +#include "GaudiAlg/Transformer.h" +#include "Kernel/LHCbID.h" +#include "LHCbMath/MatVec.h" +#include "SelKernel/Utilities.h" +#include "SelKernel/VertexRelation.h" +#include "VPDet/DeVP.h" +#include "VPDet/VPDetPaths.h" +#include <vector> + +namespace LHCb::Pr::Velo { + + namespace { + + namespace VeloInfo { + int constexpr NSensors = 208; + } // namespace VeloInfo + + namespace TrackInfo { + struct nInSensor : LHCb::Event::int_field {}; + struct hasHitInFirstSensor : LHCb::Event::int_field {}; + struct mCorr : LHCb::Event::float_field {}; + struct cosAngle : LHCb::Event::float_field {}; + } // namespace TrackInfo + + // in/outputs + using Composites = LHCb::Particles; + using VeloHits = LHCb::Pr::VP::Hits; + using PVs = LHCb::Event::PV::PrimaryVertexContainer; + using Output = LHCb::Event::RelationTable1D<Composites, TrackInfo::nInSensor, TrackInfo::hasHitInFirstSensor, + TrackInfo::mCorr, TrackInfo::cosAngle>; + + // miscellaneous + using HitPos = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>; + using SensorID = LHCb::Detector::VPChannelID::SensorID; + + // to store relevant sensor information + struct SensorInfo { + float z; + SensorID id; + friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) { + return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id}; + } + }; + using SensorInfos = std::vector<SensorInfo>; + + // to store relevant hit information + struct HitInfo { + float d; + HitPos position; + friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; } + }; + using HitInfos = std::vector<HitInfo>; + + // functions + template <typename Position3, typename Vector3> + HitInfos getHitsForFirstCrossedSensor( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV, + Vector3 const& tvec, float max_distance = 1. * Gaudi::Units::mm ) { + HitInfos poshits; + if ( !sensors.size() ) return poshits; + // get first sensor (in z) and its reference point + auto firstsensor = sensors.front(); + auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec; + // look if this sensor has hits and if they are in search window + for ( auto hit : hits.scalar() ) { + auto id = LHCb::Detector::VPChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); + if ( id.sensor() == firstsensor.id ) { + auto pos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); + auto diff = ( pos.vec3() - ref ).mag(); + if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3()} ); } + } + } + return poshits; + } + + } // namespace + + /** + * @class VeloHeavyFlavourTracking + * Search for hits of charged b/c hadrons based on primary and secondary vertices + * + * @author Maarten van Veghel (Nikhef, RUG) + */ + + class VeloHeavyFlavourTracking + : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const& ), + LHCb::DetDesc::usesConditions<DeVP>> { + public: + // standard constructor + VeloHeavyFlavourTracking( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + {KeyValue{"Composites", ""}, KeyValue{"PVs", ""}, KeyValue{"Hits", ""}, + KeyValue{"DeVP", LHCb::Det::VP::det_path}}, + {KeyValue{"OutputLocation", ""}} ) {} + + // main execution + Output operator()( Composites const& composites, PVs const& pvs, VeloHits const& hits, + DeVP const& velo ) const override { + // get z positions of sensors (should be in initialize TODO) + SensorInfos info_sensors( VeloInfo::NSensors ); + auto getSensorZs = [&info_sensors]( auto const& sensor ) { + info_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} ); + }; + velo.runOnAllSensors( getSensorZs ); + std::sort( info_sensors.begin(), info_sensors.end() ); + + // output + Output output_table( &composites ); + output_table.reserve( composites.size() ); + + for ( auto const composite : composites ) { + // get SV and PV information + using Sel::Utils::endVertexPos; + auto const& endVertex = endVertexPos( *composite ); + auto const& PV = endVertexPos( Sel::getBestPV( *composite, pvs ) ); + auto const diff = endVertex - PV; + auto const tvec = diff / diff.Z(); + + // search window crosses sensors? + int nInSensor = 0; + SensorInfos crossedSensors; + auto pos = PV; + for ( auto sinfo : info_sensors ) { + if ( sinfo.z < pos.Z() ) continue; + if ( sinfo.z > endVertex.Z() ) break; + // if in range look if search window intersects with sensor + pos = pos + ( sinfo.z - pos.Z() ) * tvec; + if ( velo.sensor( sinfo.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) { + nInSensor++; + crossedSensors.push_back( {pos.Z().cast(), sinfo.id} ); + } + } + + // look for hits in relevant sensors + // TODO maybe check all sensors in this plane, in case the search window is at the edge? + auto hits_first_sensor = getHitsForFirstCrossedSensor( crossedSensors, hits, PV, tvec ); + std::sort( hits_first_sensor.begin(), hits_first_sensor.end() ); + bool hasHitInFirstSensor = hits_first_sensor.size() > 0; + + // determine mother particle direction based on closest hit in first crossed sensor + // or using SV-PV if not available + auto direction = + hasHitInFirstSensor ? ( hits_first_sensor.front().position - PV ).normalize() : diff.normalize(); + + // obtain corrected mass and opening angle + using Sel::Utils::mass2; + using Sel::Utils::threeMomentum; + using std::sqrt; + auto const mom = threeMomentum( *composite ); + auto const perp = mom - direction * dot( mom, direction ); + auto const pt = perp.mag(); + auto const mcorr = pt + sqrt( mass2( *composite ) + pt * pt ); + auto const cosangle = dot( mom, direction ) / mom.mag(); + + // statistics + m_nInSensor += nInSensor; + if ( nInSensor > 0 ) m_nHitCandidates += hits_first_sensor.size(); + + // set output + output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor, + hasHitInFirstSensor, mcorr, cosangle ); + } + + return output_table; + } + + private: + // statistics monitoring + mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "number of sensors crossed"}; + mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "number of hit candidates per crossed sensor"}; + }; + +} // namespace LHCb::Pr::Velo + +DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::VeloHeavyFlavourTracking, "VeloHeavyFlavourTracking" ) -- GitLab From db358c253d48de604a367e6f9896668f9d259ec8 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Thu, 28 Apr 2022 16:46:35 +0200 Subject: [PATCH 02/14] searches whole module, not just sensor --- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 104 ++++++++++++-------- 1 file changed, 63 insertions(+), 41 deletions(-) diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index eb5f19e3886..81751211bed 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -48,8 +48,9 @@ namespace LHCb::Pr::Velo { TrackInfo::mCorr, TrackInfo::cosAngle>; // miscellaneous - using HitPos = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>; - using SensorID = LHCb::Detector::VPChannelID::SensorID; + using HitPos = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>; + using SensorID = LHCb::Detector::VPChannelID::SensorID; + using ChannelID = LHCb::Detector::VPChannelID; // to store relevant sensor information struct SensorInfo { @@ -65,31 +66,68 @@ namespace LHCb::Pr::Velo { struct HitInfo { float d; HitPos position; + ChannelID id; friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; } }; using HitInfos = std::vector<HitInfo>; // functions template <typename Position3, typename Vector3> - HitInfos getHitsForFirstCrossedSensor( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV, - Vector3 const& tvec, float max_distance = 1. * Gaudi::Units::mm ) { + SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& endVertex, + Position3 const& PV, Vector3 const& tvec ) { + SensorInfos crossedSensors; + auto pos = PV; + for ( auto sensor : sensors ) { + if ( sensor.z < pos.Z() ) continue; + if ( sensor.z > endVertex.Z() ) break; + // if in range look if search window intersects with sensor + pos = pos + ( sensor.z - pos.Z() ) * tvec; + if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) + crossedSensors.push_back( {pos.Z().cast(), sensor.id} ); + } + return crossedSensors; + } + + template <typename Position3, typename Vector3> + HitInfos getHitsForFirstCrossedModule( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV, + Vector3 const& tvec, float max_distance ) { HitInfos poshits; if ( !sensors.size() ) return poshits; - // get first sensor (in z) and its reference point + // get first sensor (in z) auto firstsensor = sensors.front(); - auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec; - // look if this sensor has hits and if they are in search window + // its associated module (to not be affected by borders between sensors) + auto dummy_id = ChannelID(); + dummy_id.setSensor( firstsensor.id ); + auto crossed_module = dummy_id.module(); + // first sensor reference point + auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec; + // look if this sensor (and associated station) has hits and if they are in search window for ( auto hit : hits.scalar() ) { - auto id = LHCb::Detector::VPChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); - if ( id.sensor() == firstsensor.id ) { + auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); + if ( id.module() == crossed_module ) { auto pos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); auto diff = ( pos.vec3() - ref ).mag(); - if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3()} ); } + if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3(), id} ); } } } return poshits; } + template <typename Composite_t, typename Vector3_t> + auto getMissingMomentumVars( Composite_t const& composite, Vector3_t const& direction ) { + using Sel::Utils::mass2; + using Sel::Utils::threeMomentum; + using std::sqrt; + auto const mom = threeMomentum( composite ); + // corrected mass + auto const perp = mom - direction * dot( mom, direction ); + auto const pt = perp.mag(); + auto const mcorr = pt + sqrt( mass2( composite ) + pt * pt ); + // operning angle + auto const cosangle = dot( mom, direction ) / mom.mag(); + return std::make_tuple( mcorr, cosangle ); + } + } // namespace /** @@ -134,57 +172,41 @@ namespace LHCb::Pr::Velo { auto const tvec = diff / diff.Z(); // search window crosses sensors? - int nInSensor = 0; - SensorInfos crossedSensors; - auto pos = PV; - for ( auto sinfo : info_sensors ) { - if ( sinfo.z < pos.Z() ) continue; - if ( sinfo.z > endVertex.Z() ) break; - // if in range look if search window intersects with sensor - pos = pos + ( sinfo.z - pos.Z() ) * tvec; - if ( velo.sensor( sinfo.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) { - nInSensor++; - crossedSensors.push_back( {pos.Z().cast(), sinfo.id} ); - } - } + auto crossedSensors = getCrossedSensors( velo, info_sensors, endVertex, PV, tvec ); + int nInSensor = crossedSensors.size(); // look for hits in relevant sensors - // TODO maybe check all sensors in this plane, in case the search window is at the edge? - auto hits_first_sensor = getHitsForFirstCrossedSensor( crossedSensors, hits, PV, tvec ); - std::sort( hits_first_sensor.begin(), hits_first_sensor.end() ); - bool hasHitInFirstSensor = hits_first_sensor.size() > 0; + auto first_hits = getHitsForFirstCrossedModule( crossedSensors, hits, PV, tvec, m_max_distance ); + std::sort( first_hits.begin(), first_hits.end() ); + bool hasHitInFirstModule = first_hits.size() > 0; // determine mother particle direction based on closest hit in first crossed sensor // or using SV-PV if not available - auto direction = - hasHitInFirstSensor ? ( hits_first_sensor.front().position - PV ).normalize() : diff.normalize(); + auto direction = hasHitInFirstModule ? ( first_hits.front().position - PV ).normalize() : diff.normalize(); // obtain corrected mass and opening angle - using Sel::Utils::mass2; - using Sel::Utils::threeMomentum; - using std::sqrt; - auto const mom = threeMomentum( *composite ); - auto const perp = mom - direction * dot( mom, direction ); - auto const pt = perp.mag(); - auto const mcorr = pt + sqrt( mass2( *composite ) + pt * pt ); - auto const cosangle = dot( mom, direction ) / mom.mag(); + auto [mcorr, cosangle] = getMissingMomentumVars( *composite, direction ); // statistics m_nInSensor += nInSensor; - if ( nInSensor > 0 ) m_nHitCandidates += hits_first_sensor.size(); + if ( nInSensor > 0 ) m_nHitCandidates += first_hits.size(); // set output output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor, - hasHitInFirstSensor, mcorr, cosangle ); + hasHitInFirstModule, mcorr, cosangle ); } return output_table; } private: + // properties + Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm, + "maximum hit distance from PV-SV search window"}; + // statistics monitoring - mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "number of sensors crossed"}; - mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "number of hit candidates per crossed sensor"}; + mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"}; + mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates per sensor"}; }; } // namespace LHCb::Pr::Velo -- GitLab From d840280b0848d54960600ea5f8c83719d399c38a Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Tue, 25 Oct 2022 14:49:52 +0200 Subject: [PATCH 03/14] large update; tracking (info) more general, saved in soa --- .../PrKernel/PrVeloHeavyFlavourUtils.h | 104 +++++++ Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 294 ++++++++++-------- 2 files changed, 271 insertions(+), 127 deletions(-) create mode 100644 Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h new file mode 100644 index 00000000000..5be01d6168b --- /dev/null +++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h @@ -0,0 +1,104 @@ +/*****************************************************************************\ +* (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 "Event/Particle.h" +#include "Event/Particle_v2.h" +#include "Event/PrTracksTag.h" +#include "Event/PrimaryVertices.h" +#include "Event/RelationTable.h" +#include "VPDet/DeVP.h" + +namespace LHCb::Pr::Velo { + + namespace HeavyFlavourTracking { + + // main heavy flavour track storage + namespace Tag { + // hit info + struct LHCbID : Event::lhcbid_field {}; + struct HitPosition : Event::vec3_field {}; + struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {}; + + // crossed sensor info + struct nSensors : Event::int_field {}; + + // geometry + struct OriginState : Event::state_field {}; + + // track quality info + struct Chi2 : Event::float_field {}; + + // main storage: relating track to associated composite and PV in full SoA event model + template <typename Composites, typename PVs> + using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>; + + // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects + struct CompositeKey : Event::int_field {}; + struct PrimaryVertexKey : Event::int_field {}; + + template <typename T> + using Track_t = Event::SOACollection<T, CompositeKey, PrimaryVertexKey, Hits, nSensors, OriginState, Chi2>; + } // namespace Tag + + namespace Proxy { + + // base class for proxy of for a track 'object' in the track container + template <typename ProxyBase> + struct TrackProxy : ProxyBase { + using ProxyBase::Proxy; + // hit info + [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); } + [[nodiscard]] auto hitID( std::size_t i ) const { + return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>(); + } + [[nodiscard]] auto hitPosition( std::size_t i ) const { + return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>(); + } + // crossed sensor info + [[nodiscard]] auto nSensors() const { return this->template field<Tag::nSensors>(); } + // geometry + [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); } + // track quality info + [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); } + }; + + } // namespace Proxy + + namespace v3 { + using Composites = LHCb::Event::Composites; + using PVs = LHCb::Event::PV::PrimaryVertexContainer; // FIXME + + template <typename Composites, typename PVs> + struct HeavyFlavourTracks : Tag::RelationTable2D_t<Composites, PVs> { + using Tag::RelationTable2D_t<Composites, PVs>::RelationTable2D; + + template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> + using proxy_type = Proxy::TrackProxy< + typename Tag::RelationTable2D_t<Composites, PVs>::proxy_type<simd, behaviour, ContainerType>>; + }; + } // namespace v3 + + namespace v1 { + using Composites = LHCb::Particles; + using PVs = LHCb::Event::PV::PrimaryVertexContainer; + + struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> { + using base_t = typename Tag::Track_t<HeavyFlavourTracks>; + using base_t::base_t; + + template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> + using proxy_type = Proxy::TrackProxy<typename Event::Proxy<simd, behaviour, ContainerType>>; + }; + } // namespace v1 + + } // namespace HeavyFlavourTracking + +} // namespace LHCb::Pr::Velo diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index 81751211bed..b66c44d3dab 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -11,63 +11,78 @@ #include "DetDesc/GenericConditionAccessorHolder.h" #include "DetDesc/IConditionDerivationMgr.h" -#include "Event/Particle.h" #include "Event/PrHits.h" -#include "Event/PrimaryVertices.h" -#include "Event/RawEvent.h" -#include "Event/RelationTable.h" #include "GaudiAlg/Transformer.h" #include "Kernel/LHCbID.h" #include "LHCbMath/MatVec.h" +#include "PrKernel/PrVeloHeavyFlavourUtils.h" #include "SelKernel/Utilities.h" #include "SelKernel/VertexRelation.h" -#include "VPDet/DeVP.h" #include "VPDet/VPDetPaths.h" -#include <vector> +#include "VeloKalmanHelpers.h" namespace LHCb::Pr::Velo { namespace { - namespace VeloInfo { - int constexpr NSensors = 208; - } // namespace VeloInfo - - namespace TrackInfo { - struct nInSensor : LHCb::Event::int_field {}; - struct hasHitInFirstSensor : LHCb::Event::int_field {}; - struct mCorr : LHCb::Event::float_field {}; - struct cosAngle : LHCb::Event::float_field {}; - } // namespace TrackInfo - - // in/outputs - using Composites = LHCb::Particles; - using VeloHits = LHCb::Pr::VP::Hits; - using PVs = LHCb::Event::PV::PrimaryVertexContainer; - using Output = LHCb::Event::RelationTable1D<Composites, TrackInfo::nInSensor, TrackInfo::hasHitInFirstSensor, - TrackInfo::mCorr, TrackInfo::cosAngle>; + // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA) + using namespace LHCb::Pr::Velo::HeavyFlavourTracking::v1; + using VeloHits = LHCb::Pr::VP::Hits; + using Output = HeavyFlavourTracks; // miscellaneous - using HitPos = LHCb::LinAlg::Vec3<SIMDWrapper::scalar::float_v>; - using SensorID = LHCb::Detector::VPChannelID::SensorID; - using ChannelID = LHCb::Detector::VPChannelID; + namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag; + using FType = SIMDWrapper::scalar::float_v; + using Position = LHCb::LinAlg::Vec3<FType>; + using SensorID = LHCb::Detector::VPChannelID::SensorID; + using ChannelID = LHCb::Detector::VPChannelID; // to store relevant sensor information struct SensorInfo { - float z; - SensorID id; + FType z; + SensorID id; + // sorting friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) { return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id}; } + // module information + auto module() const { + auto cid = ChannelID(); + cid.setSensor( this->id ); + return cid.module(); + } }; using SensorInfos = std::vector<SensorInfo>; + class Sensors { + private: + SensorInfos m_sensors; + + public: + // constructors + Sensors() = default; + // conditions related constructor + Sensors( DeVP const& velo ) { + m_sensors = SensorInfos{}; + m_sensors.reserve( velo.numberSensors() ); + auto getSensorZs = [&]( auto const& sensor ) { + m_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} ); + }; + velo.runOnAllSensors( getSensorZs ); + std::sort( m_sensors.begin(), m_sensors.end() ); + }; + // main access + SensorInfos const& operator()() const { return m_sensors; } + }; + // to store relevant hit information struct HitInfo { - float d; - HitPos position; + FType d; + Position position; ChannelID id; - friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.d < rhs.d; } + friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { + return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()}; + } }; using HitInfos = std::vector<HitInfo>; @@ -89,126 +104,151 @@ namespace LHCb::Pr::Velo { } template <typename Position3, typename Vector3> - HitInfos getHitsForFirstCrossedModule( SensorInfos const& sensors, VeloHits const& hits, Position3 const& PV, - Vector3 const& tvec, float max_distance ) { + HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start, + Vector3 const& tvec, FType max_distance ) { HitInfos poshits; if ( !sensors.size() ) return poshits; - // get first sensor (in z) - auto firstsensor = sensors.front(); - // its associated module (to not be affected by borders between sensors) - auto dummy_id = ChannelID(); - dummy_id.setSensor( firstsensor.id ); - auto crossed_module = dummy_id.module(); - // first sensor reference point - auto ref = PV + ( firstsensor.z - PV.Z() ) * tvec; - // look if this sensor (and associated station) has hits and if they are in search window - for ( auto hit : hits.scalar() ) { - auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); - if ( id.module() == crossed_module ) { - auto pos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); - auto diff = ( pos.vec3() - ref ).mag(); - if ( diff < max_distance ) { poshits.push_back( {diff.cast(), pos.vec3(), id} ); } + poshits.reserve( sensors.size() ); + // starting point search + auto pos = start; + // check sensors + for ( auto const sensor : sensors ) { + // its associated module (to not be affected by borders between sensors) + auto crossed_module = sensor.module(); + // go to sensor location + pos = pos + ( sensor.z - pos.z.cast() ) * tvec; + // look if this sensor (and associated station) has hits and if they are in search window + for ( auto hit : hits.scalar() ) { + auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); + if ( id.module() == crossed_module ) { + auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); + auto diff = ( hitpos.vec3() - pos ).mag(); + if ( diff < max_distance ) poshits.push_back( {diff.cast(), hitpos.vec3(), id} ); + } } } return poshits; } - template <typename Composite_t, typename Vector3_t> - auto getMissingMomentumVars( Composite_t const& composite, Vector3_t const& direction ) { - using Sel::Utils::mass2; - using Sel::Utils::threeMomentum; - using std::sqrt; - auto const mom = threeMomentum( composite ); - // corrected mass - auto const perp = mom - direction * dot( mom, direction ); - auto const pt = perp.mag(); - auto const mcorr = pt + sqrt( mass2( composite ) + pt * pt ); - // operning angle - auto const cosangle = dot( mom, direction ) / mom.mag(); - return std::make_tuple( mcorr, cosangle ); - } - } // namespace /** - * @class VeloHeavyFlavourTracking + * @class HeavyFlavourTrackFinder * Search for hits of charged b/c hadrons based on primary and secondary vertices * * @author Maarten van Veghel (Nikhef, RUG) */ - class VeloHeavyFlavourTracking - : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const& ), - LHCb::DetDesc::usesConditions<DeVP>> { + class HeavyFlavourTrackFinder + : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const&, + Sensors const& ), + LHCb::DetDesc::usesConditions<DeVP, Sensors>> { public: // standard constructor - VeloHeavyFlavourTracking( const std::string& name, ISvcLocator* pSvcLocator ) - : Transformer( name, pSvcLocator, - {KeyValue{"Composites", ""}, KeyValue{"PVs", ""}, KeyValue{"Hits", ""}, - KeyValue{"DeVP", LHCb::Det::VP::det_path}}, - {KeyValue{"OutputLocation", ""}} ) {} - - // main execution - Output operator()( Composites const& composites, PVs const& pvs, VeloHits const& hits, - DeVP const& velo ) const override { - // get z positions of sensors (should be in initialize TODO) - SensorInfos info_sensors( VeloInfo::NSensors ); - auto getSensorZs = [&info_sensors]( auto const& sensor ) { - info_sensors.push_back( {(float)sensor.z(), sensor.sensorNumber()} ); - }; - velo.runOnAllSensors( getSensorZs ); - std::sort( info_sensors.begin(), info_sensors.end() ); - - // output - Output output_table( &composites ); - output_table.reserve( composites.size() ); - - for ( auto const composite : composites ) { - // get SV and PV information - using Sel::Utils::endVertexPos; - auto const& endVertex = endVertexPos( *composite ); - auto const& PV = endVertexPos( Sel::getBestPV( *composite, pvs ) ); - auto const diff = endVertex - PV; - auto const tvec = diff / diff.Z(); - - // search window crosses sensors? - auto crossedSensors = getCrossedSensors( velo, info_sensors, endVertex, PV, tvec ); - int nInSensor = crossedSensors.size(); - - // look for hits in relevant sensors - auto first_hits = getHitsForFirstCrossedModule( crossedSensors, hits, PV, tvec, m_max_distance ); - std::sort( first_hits.begin(), first_hits.end() ); - bool hasHitInFirstModule = first_hits.size() > 0; - - // determine mother particle direction based on closest hit in first crossed sensor - // or using SV-PV if not available - auto direction = hasHitInFirstModule ? ( first_hits.front().position - PV ).normalize() : diff.normalize(); - - // obtain corrected mass and opening angle - auto [mcorr, cosangle] = getMissingMomentumVars( *composite, direction ); - - // statistics - m_nInSensor += nInSensor; - if ( nInSensor > 0 ) m_nHitCandidates += first_hits.size(); - - // set output - output_table.emplace_back<SIMDWrapper::InstructionSet::Scalar>().set( composite->key(), nInSensor, - hasHitInFirstModule, mcorr, cosangle ); - } + HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator ); - return output_table; - } + // initialize + StatusCode initialize() override; + + // main function/operator + Output operator()( Composites const&, PVs const&, VeloHits const&, DeVP const&, Sensors const& ) const override; private: // properties Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm, - "maximum hit distance from PV-SV search window"}; + "maximum hit distance from PV-SV line (search window)"}; // statistics monitoring mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"}; - mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates per sensor"}; - }; + mutable Gaudi::Accumulators::StatCounter<> m_nHitCandidates{this, "Nb hit candidates"}; + }; // namespace LHCb::Pr::Velo + + DECLARE_COMPONENT_WITH_ID( HeavyFlavourTrackFinder, "VeloHeavyFlavourTrackFinder" ) + + // ============================= IMPLEMENTATION =================================== + + // constructor + HeavyFlavourTrackFinder::HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator ) + : Transformer( name, pSvcLocator, + {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ), + KeyValue( "DeVP", LHCb::Det::VP::det_path ), +#ifdef USE_DD4HEP + KeyValue( "Sensors", {"/world:AlgorithmSpecific-" + name + "-sensorinfos"} )}, +#else + KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )}, +#endif + + {KeyValue{"OutputLocation", ""}} ) { + } + + // initialization + StatusCode LHCb::Pr::Velo::HeavyFlavourTrackFinder::initialize() { + auto sc = Transformer::initialize().andThen( [&] { + addConditionDerivation<Sensors( DeVP const& )>( {LHCb::Det::VP::det_path}, inputLocation<Sensors>() ); + } ); + return sc; + } + + // main execution + Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( Composites const& composites, PVs const& pvs, + VeloHits const& velohits, DeVP const& velo, + Sensors const& sensors ) const { + // output + Output output{}; // composites.zipIdentifier()}; + output.reserve( composites.size() ); + + for ( auto const composite : composites ) { + // set output + auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>(); + track.field<TrackTag::CompositeKey>().set( composite->key() ); + + // primary vertex + using Sel::Utils::endVertexPos; + auto const& PV = Sel::getBestPV( *composite, pvs ); + auto const& PVpos = endVertexPos( PV ); + track.field<TrackTag::CompositeKey>().set( PV.key() ); + + // secondary vertex + auto const& endVertex = endVertexPos( *composite ); + auto const diff = endVertex - PVpos; + auto const tvec = diff / diff.Z(); + + // crossed sensors in PV-SV line + auto crossedSensors = getCrossedSensors( velo, sensors(), endVertex, PVpos, tvec ); + track.field<TrackTag::nSensors>().set( crossedSensors.size() ); + + // look for hits in relevant sensors and sort in z and distance to search window + auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, PVpos, tvec, m_max_distance.value() ); + std::sort( found_hits.begin(), found_hits.end() ); + + // store sorted hits + auto saved_hits = track.field<TrackTag::Hits>(); + auto nhits = found_hits.size(); + saved_hits.resize( nhits ); + for ( long unsigned int i = 0; i < nhits; i++ ) { + auto found_hit = found_hits[i]; + saved_hits[i].field<TrackTag::LHCbID>().set( LHCbID( found_hit.id ) ); + saved_hits[i].field<TrackTag::HitPosition>().set( found_hit.position ); + } -} // namespace LHCb::Pr::Velo + // determine mother particle direction based on closest hit + // in first crossed sensor or using SV-PV if not available + auto direction = found_hits.size() ? found_hits.front().position - PVpos : diff; + auto origin_state = track.field<TrackTag::OriginState>(); + origin_state.setPosition( PVpos ); + origin_state.setDirection( direction.x / direction.z, direction.y / direction.z ); + origin_state.setQOverP( 1. / composite->p() ); + + // fitting + // TODO use stuff from VeloKalmanHelpers + track.field<TrackTag::Chi2>().set( -1.f ); + + // statistics + m_nInSensor += crossedSensors.size(); + m_nHitCandidates += nhits; + } -DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::VeloHeavyFlavourTracking, "VeloHeavyFlavourTracking" ) + return output; + } + +} // namespace LHCb::Pr::Velo -- GitLab From 562404700ab34cda9329ded7af9f080f620bc7af Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Tue, 25 Oct 2022 22:37:24 +0200 Subject: [PATCH 04/14] hit duplicate removal --- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 63 ++++++++++++--------- 1 file changed, 36 insertions(+), 27 deletions(-) diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index b66c44d3dab..b96369a569c 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -77,26 +77,29 @@ namespace LHCb::Pr::Velo { // to store relevant hit information struct HitInfo { - FType d; - Position position; - ChannelID id; + FType d; + Position position; + ChannelID id; + // sorting in z first, then distance to search window friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()}; } + // equal just by unique channel id + friend bool operator==( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.id == rhs.id; } }; using HitInfos = std::vector<HitInfo>; // functions template <typename Position3, typename Vector3> - SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& endVertex, - Position3 const& PV, Vector3 const& tvec ) { + SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start, + Position3 const& end, Vector3 const& slopes ) { SensorInfos crossedSensors; - auto pos = PV; + auto pos = start; for ( auto sensor : sensors ) { if ( sensor.z < pos.Z() ) continue; - if ( sensor.z > endVertex.Z() ) break; + if ( sensor.z > end.Z() ) break; // if in range look if search window intersects with sensor - pos = pos + ( sensor.z - pos.Z() ) * tvec; + pos = pos + ( sensor.z - pos.Z() ) * slopes; if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) crossedSensors.push_back( {pos.Z().cast(), sensor.id} ); } @@ -105,10 +108,10 @@ namespace LHCb::Pr::Velo { template <typename Position3, typename Vector3> HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start, - Vector3 const& tvec, FType max_distance ) { + Vector3 const& slopes, FType max_distance ) { HitInfos poshits; if ( !sensors.size() ) return poshits; - poshits.reserve( sensors.size() ); + poshits.reserve( 2 * sensors.size() ); // starting point search auto pos = start; // check sensors @@ -116,17 +119,24 @@ namespace LHCb::Pr::Velo { // its associated module (to not be affected by borders between sensors) auto crossed_module = sensor.module(); // go to sensor location - pos = pos + ( sensor.z - pos.z.cast() ) * tvec; - // look if this sensor (and associated station) has hits and if they are in search window + pos = pos + ( sensor.z - pos.z.cast() ) * slopes; + // look if this sensor (and associated module) has hits and if they are in search window for ( auto hit : hits.scalar() ) { auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); if ( id.module() == crossed_module ) { auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); - auto diff = ( hitpos.vec3() - pos ).mag(); - if ( diff < max_distance ) poshits.push_back( {diff.cast(), hitpos.vec3(), id} ); + // make sure they are in the same plane up to some tolerance + // (not all sensors in a module are exactly) + if ( abs( hitpos.vec3().z.cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) { + auto d = ( hitpos.vec3() - pos ).mag(); + if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos.vec3(), id} ); + } } } } + // sort and remove duplicates (can arise with only module/plane (not sensor) checking) + std::sort( poshits.begin(), poshits.end() ); + poshits.erase( std::unique( poshits.begin(), poshits.end() ), poshits.end() ); return poshits; } @@ -194,7 +204,7 @@ namespace LHCb::Pr::Velo { VeloHits const& velohits, DeVP const& velo, Sensors const& sensors ) const { // output - Output output{}; // composites.zipIdentifier()}; + Output output{}; output.reserve( composites.size() ); for ( auto const composite : composites ) { @@ -204,22 +214,21 @@ namespace LHCb::Pr::Velo { // primary vertex using Sel::Utils::endVertexPos; - auto const& PV = Sel::getBestPV( *composite, pvs ); - auto const& PVpos = endVertexPos( PV ); + auto const& PV = Sel::getBestPV( *composite, pvs ); + auto const& pv = endVertexPos( PV ); track.field<TrackTag::CompositeKey>().set( PV.key() ); // secondary vertex - auto const& endVertex = endVertexPos( *composite ); - auto const diff = endVertex - PVpos; - auto const tvec = diff / diff.Z(); + auto const& end = endVertexPos( *composite ); + auto const fd = end - pv; + auto const slopes = fd / fd.Z(); // crossed sensors in PV-SV line - auto crossedSensors = getCrossedSensors( velo, sensors(), endVertex, PVpos, tvec ); + auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes ); track.field<TrackTag::nSensors>().set( crossedSensors.size() ); - // look for hits in relevant sensors and sort in z and distance to search window - auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, PVpos, tvec, m_max_distance.value() ); - std::sort( found_hits.begin(), found_hits.end() ); + // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window) + auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); // store sorted hits auto saved_hits = track.field<TrackTag::Hits>(); @@ -233,14 +242,14 @@ namespace LHCb::Pr::Velo { // determine mother particle direction based on closest hit // in first crossed sensor or using SV-PV if not available - auto direction = found_hits.size() ? found_hits.front().position - PVpos : diff; + auto direction = found_hits.size() ? found_hits.front().position - pv : fd; auto origin_state = track.field<TrackTag::OriginState>(); - origin_state.setPosition( PVpos ); + origin_state.setPosition( pv ); origin_state.setDirection( direction.x / direction.z, direction.y / direction.z ); origin_state.setQOverP( 1. / composite->p() ); // fitting - // TODO use stuff from VeloKalmanHelpers + // TODO (using VeloKalmanHelpers) track.field<TrackTag::Chi2>().set( -1.f ); // statistics -- GitLab From e17338cbe7adde486d55a84b57cac4572dc0ddfc Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Tue, 1 Nov 2022 18:54:41 +0100 Subject: [PATCH 05/14] added mc checker of hf tracking --- .../PrKernel/PrVeloHeavyFlavourUtils.h | 83 +++--- Pr/PrMCTools/CMakeLists.txt | 2 + .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 249 ++++++++++++++++++ Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 4 +- 4 files changed, 291 insertions(+), 47 deletions(-) create mode 100644 Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h index 5be01d6168b..941c308233b 100644 --- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h +++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h @@ -36,66 +36,59 @@ namespace LHCb::Pr::Velo { // track quality info struct Chi2 : Event::float_field {}; - // main storage: relating track to associated composite and PV in full SoA event model - template <typename Composites, typename PVs> - using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>; - // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects - struct CompositeKey : Event::int_field {}; - struct PrimaryVertexKey : Event::int_field {}; + struct CompositeIndex : Event::int_field {}; + struct PrimaryVertexIndex : Event::int_field {}; template <typename T> - using Track_t = Event::SOACollection<T, CompositeKey, PrimaryVertexKey, Hits, nSensors, OriginState, Chi2>; - } // namespace Tag - - namespace Proxy { - - // base class for proxy of for a track 'object' in the track container - template <typename ProxyBase> - struct TrackProxy : ProxyBase { - using ProxyBase::Proxy; - // hit info - [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); } - [[nodiscard]] auto hitID( std::size_t i ) const { - return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>(); - } - [[nodiscard]] auto hitPosition( std::size_t i ) const { - return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>(); - } - // crossed sensor info - [[nodiscard]] auto nSensors() const { return this->template field<Tag::nSensors>(); } - // geometry - [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); } - // track quality info - [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); } - }; - - } // namespace Proxy - - namespace v3 { - using Composites = LHCb::Event::Composites; - using PVs = LHCb::Event::PV::PrimaryVertexContainer; // FIXME + using Track_t = Event::SOACollection<T, CompositeIndex, PrimaryVertexIndex, Hits, nSensors, OriginState, Chi2>; + /* + // main storage: relating track to associated composite and PV in full SoA event model template <typename Composites, typename PVs> - struct HeavyFlavourTracks : Tag::RelationTable2D_t<Composites, PVs> { - using Tag::RelationTable2D_t<Composites, PVs>::RelationTable2D; - - template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> - using proxy_type = Proxy::TrackProxy< - typename Tag::RelationTable2D_t<Composites, PVs>::proxy_type<simd, behaviour, ContainerType>>; - }; - } // namespace v3 + using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>; + */ + } // namespace Tag namespace v1 { using Composites = LHCb::Particles; + using Composite = LHCb::Particle; using PVs = LHCb::Event::PV::PrimaryVertexContainer; struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> { using base_t = typename Tag::Track_t<HeavyFlavourTracks>; using base_t::base_t; + template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType> + struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> { + using base_t = typename Event::Proxy<simd, behaviour, ContainerType>; + using base_t::Proxy; + + // composite / pv info + [[nodiscard]] auto compositeIndex() const { return this->template get<Tag::CompositeIndex>(); } + [[nodiscard]] auto primaryVertexIndex() const { return this->template get<Tag::PrimaryVertexIndex>(); } + + // hit info + [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); } + [[nodiscard]] auto hitID( std::size_t i ) const { + return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>(); + } + [[nodiscard]] auto hitPosition( std::size_t i ) const { + return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>(); + } + + // crossed sensor info + [[nodiscard]] auto nSensors() const { return this->template get<Tag::nSensors>(); } + + // geometry + [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); } + + // track quality info + [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); } + }; + template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> - using proxy_type = Proxy::TrackProxy<typename Event::Proxy<simd, behaviour, ContainerType>>; + using proxy_type = TrackProxy<simd, behaviour, ContainerType>; }; } // namespace v1 diff --git a/Pr/PrMCTools/CMakeLists.txt b/Pr/PrMCTools/CMakeLists.txt index 62b61171ae9..4cae29bd6a7 100644 --- a/Pr/PrMCTools/CMakeLists.txt +++ b/Pr/PrMCTools/CMakeLists.txt @@ -38,6 +38,7 @@ gaudi_add_module(PrMCTools src/PrVPHitsChecker.cpp src/PrUTHitsChecker.cpp src/PrFTHitsChecker.cpp + src/PrVeloHeavyFlavourTrackingChecker.cpp src/VPClusterEfficiency.cpp src/PrGhostDumper.cpp src/PrDownstreamDumper.cpp @@ -70,6 +71,7 @@ gaudi_add_module(PrMCTools Rec::PrFitParamsLib Rec::PrKernel Rec::TrackInterfacesLib + Rec::SelKernelLib ROOT::Core ROOT::RIO ROOT::Tree diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp new file mode 100644 index 00000000000..503e2afedff --- /dev/null +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -0,0 +1,249 @@ +/*****************************************************************************\ +* (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 "DetDesc/GenericConditionAccessorHolder.h" +#include "Event/LinksByKey.h" +#include "Event/MCHit.h" +#include "Event/MCParticle.h" +#include "Event/RawEvent.h" +#include "Event/VPDigit.h" +#include "Event/VPFullCluster.h" +#include "GaudiAlg/Consumer.h" +#include "GaudiAlg/GaudiTupleAlg.h" +#include "Linker/LinkedFrom.h" +#include "Linker/LinkedTo.h" +#include "PrKernel/PrVeloHeavyFlavourUtils.h" + +/** @class PrVeloHeavyFlavourTrackingChecker PrVeloHeavyFlavourTrackingChecker.cpp + * + * For efficiency checker of VeloHeavyFlavourTrackFinder on MC + * + */ + +namespace { + // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA) + using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks; + using Composites = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites; + using Composite = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite; + using MCParticles = LHCb::MCParticles; + using MCHits = LHCb::MCHits; + using Links = LHCb::LinksByKey; + + // helper classes + struct DecayMCParticles { + LHCb::MCParticle const* B = NULL; + LHCb::MCParticle const* tau = NULL; + std::array<LHCb::MCParticle const*, 3> pions = {NULL, NULL, NULL}; + }; + + // mcparticle filter + std::optional<DecayMCParticles> has_b2taunu_tau23pions( LHCb::MCParticle const* mcp ) { + if ( mcp->particleID().hasBottom() ) { + auto vtxs = mcp->endVertices(); + if ( vtxs.size() ) { + auto sv = vtxs.at( vtxs.size() - 1 ); + auto daus = sv->products(); + if ( daus.size() >= 2 ) { + bool has_taup = false; + bool has_ntau = false; + size_t n_photons = 0; + for ( LHCb::MCParticle const* dau : daus ) { + auto dID = dau->particleID().abspid(); + if ( dID == 15 ) has_taup = true; + if ( dID == 16 ) has_ntau = true; + if ( dID == 22 ) n_photons += 1; + } + if ( has_taup && has_ntau && ( n_photons == ( daus.size() - 2 ) ) ) { + auto parts = DecayMCParticles(); + parts.B = mcp; + for ( LHCb::MCParticle const* dau : daus ) { + if ( dau->particleID().abspid() == 15 ) parts.tau = dau; + } + auto tau_vtxs = parts.tau->endVertices(); + if ( tau_vtxs.size() ) { + auto tv = tau_vtxs.at( tau_vtxs.size() - 1 ); + int i_pion = 0; + for ( LHCb::MCParticle const* tdau : tv->products() ) { + if ( tdau->particleID().abspid() == 211 ) { + parts.pions[i_pion] = tdau; + i_pion += 1; + } + } + } + return parts; + } // getting all relevant MCParticles + } + } + } + return std::nullopt; + } + + // match MCParticle to Composites + template <typename LinkMap> + std::optional<Composite const*> matchToComposites( DecayMCParticles const& mcps, Composites const& composites, + LinkMap const& map, MCParticles const& mcparts ) { + for ( Composite const* comp : composites ) { + size_t nmatches = 0; + for ( Composite const* dau : comp->daughtersVector() ) { + if ( !dau->proto() ) continue; + for ( LHCb::MCParticle const* pion : mcps.pions ) { + map.applyToLinks( dau->proto()->track()->key(), [&nmatches, &mcparts, + &pion]( unsigned int, unsigned int mcPartKey, float ) { + if ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) nmatches += 1; + } ); + } + } + // assuming there is only one matching composite + if ( nmatches == mcps.pions.size() ) return comp; + } + return std::nullopt; + } + + // match MCParticles to MCHits + auto getMCHits( DecayMCParticles const& mcps, MCHits const& mchits ) { + std::vector<LHCb::MCHit const*> hits; + // which ones come from B or tau? + for ( LHCb::MCHit const* mchit : mchits ) { + auto hitmcp = mchit->mcParticle(); + if ( hitmcp == mcps.B || hitmcp == mcps.tau ) hits.push_back( mchit ); + } + // order them in z + auto ordering = []( LHCb::MCHit const* i1, LHCb::MCHit const* i2 ) { return i1->entry().z() < i2->entry().z(); }; + std::sort( hits.begin(), hits.end(), ordering ); + return hits; + } + + // get track + std::optional<int> getHeavyFlavourTrack( HeavyFlavourTracks const& tracks, int key ) { + for ( auto track : tracks.scalar() ) { + if ( track.compositeIndex() == key ) return track.indices().cast(); + } + return std::nullopt; + } + +} // namespace + +class PrVeloHeavyFlavourTrackingChecker + : public Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&, + MCHits const&, Links const&, Links const& ), + LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>> { +public: + using base_type = Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&, + MCHits const&, Links const&, Links const& ), + LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>; + using KeyValue = base_type::KeyValue; + + // standard constructor + PrVeloHeavyFlavourTrackingChecker( std::string const& name, ISvcLocator* pSvc ) + : base_type{name, + pSvc, + {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"Composites", ""}, + KeyValue{"HeavyFlavourTracks", ""}, KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""}, + KeyValue{"TrackLinks", ""}}} {} + + // main execution + void operator()( MCParticles const&, Composites const&, HeavyFlavourTracks const&, MCHits const&, Links const&, + Links const& ) const override; + +private: + // monitoring + mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""}; +}; + +DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavourTrackingChecker" ) + +//// implementation ///// + +// main execution +void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, Composites const& composites, + HeavyFlavourTracks const& tracks, MCHits const& mchits, + Links const& vplinks, Links const& tracklinks ) const { + // setup ntuple + auto tuple = this->nTuple( "MCParticleTuple" ); + + // table linking a LHCb::MCHit* to std::vector<LHCb::Detector::VPChannelID>> and vice versa + std::map<LHCb::MCHit const*, std::vector<unsigned int>> channelIDForMCHit; + vplinks.applyToAllLinks( [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float ) { + channelIDForMCHit[mchits[mcHitKey]].emplace_back( channelID ); + } ); + + // fill for all tracks + for ( auto const mcp : mcparts ) { + if ( !mcp ) continue; + // check if it is the right MCParticle + auto check_mcps = has_b2taunu_tau23pions( mcp ); + if ( !check_mcps.has_value() ) continue; + auto mcps = check_mcps.value(); + // get MCHits for the candidate + auto bhits = getMCHits( mcps, mchits ); + auto first_mchit = bhits.size() ? bhits.front()->midPoint() : Gaudi::XYZPoint( 0., 0., 0. ); + // first link MCParticle to daughters (pions) of composite + auto composite = matchToComposites( mcps, composites, tracklinks, mcparts ); + auto has_comp = composite.has_value(); + // get heavy flavour track + auto itrack = has_comp ? getHeavyFlavourTrack( tracks, composite.value()->key() ) : std::nullopt; + auto has_track = has_comp ? itrack.has_value() : false; + auto has_hits = has_track ? (int)tracks.scalar()[itrack.value()].nHits() > 0 : false; + //// fill ntuple with info + std::string base = "MCP_"; + // mc truth + auto sc = tuple->column( base + "TRUEID", mcp->particleID().pid() ); + sc &= tuple->column( base + "P", mcp->momentum().P() ); + // mc hits + sc &= tuple->column( base + "has_VP_MCHits", (bool)bhits.size() ); + sc &= tuple->column( base + "VP_first_MCHit_x", first_mchit.x() ); + sc &= tuple->column( base + "VP_first_MCHit_y", first_mchit.y() ); + sc &= tuple->column( base + "VP_first_MCHit_z", first_mchit.z() ); + // composite info + sc &= tuple->column( base + "has_Composite", has_comp ); + std::string compbase = base + "Composite_"; + auto endvtx = has_comp ? composite.value()->endVertex()->position() : Gaudi::XYZPoint( 0., 0., 0. ); + sc &= tuple->column( compbase + "EndVertex_x", endvtx.x() ); + sc &= tuple->column( compbase + "EndVertex_y", endvtx.y() ); + sc &= tuple->column( compbase + "EndVertex_z", endvtx.z() ); + auto compmom = + has_comp ? composite.value()->momentum() : Gaudi::LorentzVector( 0., 0., 0., -1 * Gaudi::Units::GeV ); + sc &= tuple->column( compbase + "P", compmom.R() ); + sc &= tuple->column( compbase + "PT", compmom.Pt() ); + sc &= tuple->column( compbase + "eta", compmom.Eta() ); + sc &= tuple->column( compbase + "phi", compmom.Phi() ); + // track info + std::string hfbase = base + "HeavyFlavourTracking_"; + sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks.scalar()[itrack.value()].nHits() : -1 ); + // track mc matching info + bool fh_from_mcps = false; + int fh_trueid = 0; + if ( has_hits ) { + float max_weight = 0.f; + LHCb::MCHit const* mchit = nullptr; + unsigned int hid = tracks.scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID(); + vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) { + if ( weight > max_weight ) { + max_weight = weight; + mchit = static_cast<LHCb::MCHit const*>( mchits.containedObject( mcHitKey ) ); + } + } ); + if ( mchit ) { + fh_trueid = mchit->mcParticle()->particleID().pid(); + fh_from_mcps = ( mchit->mcParticle() == mcps.B ) || ( mchit->mcParticle() == mcps.tau ); + } + } + sc &= tuple->column( hfbase + "Track_first_hit_x", + has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. ); + sc &= tuple->column( hfbase + "Track_first_hit_y", + has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. ); + sc &= tuple->column( hfbase + "Track_first_hit_z", + has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. ); + sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps ); + sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid ); + // write result + sc.andThen( [&] { return tuple->write(); } ).ignore(); + } +} diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index b96369a569c..5fab6af5c4b 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -210,13 +210,13 @@ namespace LHCb::Pr::Velo { for ( auto const composite : composites ) { // set output auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>(); - track.field<TrackTag::CompositeKey>().set( composite->key() ); + track.field<TrackTag::CompositeIndex>().set( composite->key() ); // primary vertex using Sel::Utils::endVertexPos; auto const& PV = Sel::getBestPV( *composite, pvs ); auto const& pv = endVertexPos( PV ); - track.field<TrackTag::CompositeKey>().set( PV.key() ); + track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() ); // secondary vertex auto const& end = endVertexPos( *composite ); -- GitLab From 8fee8c4381110ab43bae8fe81cc1e7adfdccfee2 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Tue, 8 Nov 2022 14:16:42 +0100 Subject: [PATCH 06/14] hf track mc checker running on all mcparticles --- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 227 +++++++++++++++--- 1 file changed, 190 insertions(+), 37 deletions(-) diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index 503e2afedff..a7caf773073 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -13,13 +13,17 @@ #include "Event/MCHit.h" #include "Event/MCParticle.h" #include "Event/RawEvent.h" +#include "Event/Track_v3.h" #include "Event/VPDigit.h" #include "Event/VPFullCluster.h" +#include "Functors/Common.h" #include "GaudiAlg/Consumer.h" #include "GaudiAlg/GaudiTupleAlg.h" #include "Linker/LinkedFrom.h" #include "Linker/LinkedTo.h" #include "PrKernel/PrVeloHeavyFlavourUtils.h" +#include "SelKernel/Utilities.h" +#include "SelKernel/VertexRelation.h" /** @class PrVeloHeavyFlavourTrackingChecker PrVeloHeavyFlavourTrackingChecker.cpp * @@ -32,15 +36,24 @@ namespace { using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks; using Composites = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites; using Composite = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite; + using PVs = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs; using MCParticles = LHCb::MCParticles; using MCHits = LHCb::MCHits; using Links = LHCb::LinksByKey; // helper classes struct DecayMCParticles { - LHCb::MCParticle const* B = NULL; - LHCb::MCParticle const* tau = NULL; - std::array<LHCb::MCParticle const*, 3> pions = {NULL, NULL, NULL}; + LHCb::MCParticle const* B = NULL; + LHCb::MCParticle const* tau = NULL; + std::array<LHCb::MCParticle const*, 3> pions = {NULL, NULL, NULL}; + std::array<Composite const*, 3> recpions = {NULL, NULL, NULL}; + }; + + struct WeightedHit { + unsigned int id; + float weight; + // reverse sorting (high to low) + friend bool operator<( WeightedHit& lhs, WeightedHit& rhs ) { return lhs.weight > rhs.weight; } }; // mcparticle filter @@ -85,18 +98,50 @@ namespace { return std::nullopt; } + // check if MCParticle has match to track + template <typename LinkMap> + bool hasTrack( LHCb::MCParticle const* mcp, LinkMap const& map, float min_weight = 0.7 ) { + bool has_track = false; + map.applyToAllLinks( [&has_track, &mcp, &min_weight]( unsigned int, unsigned int mcPartKey, float weight ) { + if ( ( (unsigned int)mcp->key() == mcPartKey ) && weight > min_weight ) has_track = true; + } ); + return has_track; + } + + // match MCParticle daughters (pions from tau) to Tracks (general), 'reconstructed' variable + template <typename LinkMap> + bool matchToTracks( DecayMCParticles const& mcps, LinkMap const& map, MCParticles const& mcparts, + float min_weight = 0.7 ) { + size_t nmatches = 0; + for ( LHCb::MCParticle const* pion : mcps.pions ) { + map.applyToAllLinks( + [&pion, &mcps, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) { + if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) && + ( weight > min_weight ) ) + nmatches += 1; + } ); + } + return nmatches == mcps.pions.size(); + } + // match MCParticle to Composites template <typename LinkMap> - std::optional<Composite const*> matchToComposites( DecayMCParticles const& mcps, Composites const& composites, - LinkMap const& map, MCParticles const& mcparts ) { + std::optional<Composite const*> matchToComposites( DecayMCParticles& mcps, Composites const& composites, + LinkMap const& map, MCParticles const& mcparts, + float min_weight = 0.7 ) { for ( Composite const* comp : composites ) { size_t nmatches = 0; for ( Composite const* dau : comp->daughtersVector() ) { if ( !dau->proto() ) continue; - for ( LHCb::MCParticle const* pion : mcps.pions ) { - map.applyToLinks( dau->proto()->track()->key(), [&nmatches, &mcparts, - &pion]( unsigned int, unsigned int mcPartKey, float ) { - if ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) nmatches += 1; + for ( int i = 0; i < (int)mcps.pions.size(); i++ ) { + LHCb::MCParticle const* pion = mcps.pions[i]; + map.applyToLinks( dau->proto()->track()->key(), [&i, &dau, &mcps, &nmatches, &mcparts, &pion, &min_weight]( + unsigned int, unsigned int mcPartKey, float weight ) { + if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) && + ( weight >= min_weight ) ) { + mcps.recpions[i] = dau; + nmatches += 1; + } } ); } } @@ -116,7 +161,7 @@ namespace { } // order them in z auto ordering = []( LHCb::MCHit const* i1, LHCb::MCHit const* i2 ) { return i1->entry().z() < i2->entry().z(); }; - std::sort( hits.begin(), hits.end(), ordering ); + if ( hits.size() > 1 ) std::sort( hits.begin(), hits.end(), ordering ); return hits; } @@ -128,31 +173,41 @@ namespace { return std::nullopt; } + // composite observables + float getMcorr( Composite const* composite, Gaudi::XYZVector const& direction ) { + auto mom = composite->momentum().Vect(); + auto mass2 = composite->momentum().mass2(); + auto pt = ( mom - direction * mom.Dot( direction ) ).r(); + return pt + std::sqrt( mass2 + pt * pt ); + } + } // namespace class PrVeloHeavyFlavourTrackingChecker - : public Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&, - MCHits const&, Links const&, Links const& ), + : public Gaudi::Functional::Consumer<void( MCParticles const&, PVs const&, MCHits const&, Links const&, + Links const& ), LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>> { public: - using base_type = Gaudi::Functional::Consumer<void( MCParticles const&, Composites const&, HeavyFlavourTracks const&, - MCHits const&, Links const&, Links const& ), - LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>; - using KeyValue = base_type::KeyValue; + using base_type = + Gaudi::Functional::Consumer<void( MCParticles const&, PVs const&, MCHits const&, Links const&, Links const& ), + LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg>>; + using KeyValue = base_type::KeyValue; // standard constructor PrVeloHeavyFlavourTrackingChecker( std::string const& name, ISvcLocator* pSvc ) : base_type{name, pSvc, - {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"Composites", ""}, - KeyValue{"HeavyFlavourTracks", ""}, KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""}, - KeyValue{"TrackLinks", ""}}} {} + {KeyValue{"MCParticles", LHCb::MCParticleLocation::Default}, KeyValue{"PVs", ""}, + KeyValue{"MCHits", "/Event/MC/VP/Hits"}, KeyValue{"VPLinks", ""}, KeyValue{"TrackLinks", ""}}} {} // main execution - void operator()( MCParticles const&, Composites const&, HeavyFlavourTracks const&, MCHits const&, Links const&, - Links const& ) const override; + void operator()( MCParticles const&, PVs const&, MCHits const&, Links const&, Links const& ) const override; private: + // data that might not be there + DataObjectReadHandle<Composites> m_composites{this, "Composites", ""}; + DataObjectReadHandle<HeavyFlavourTracks> m_tracks{this, "HeavyFlavourTracks", ""}; + // monitoring mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""}; }; @@ -162,17 +217,17 @@ DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavou //// implementation ///// // main execution -void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, Composites const& composites, - HeavyFlavourTracks const& tracks, MCHits const& mchits, +void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, PVs const& pvs, MCHits const& mchits, Links const& vplinks, Links const& tracklinks ) const { // setup ntuple auto tuple = this->nTuple( "MCParticleTuple" ); // table linking a LHCb::MCHit* to std::vector<LHCb::Detector::VPChannelID>> and vice versa - std::map<LHCb::MCHit const*, std::vector<unsigned int>> channelIDForMCHit; - vplinks.applyToAllLinks( [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float ) { - channelIDForMCHit[mchits[mcHitKey]].emplace_back( channelID ); - } ); + std::map<LHCb::MCHit const*, std::vector<WeightedHit>> channelIDForMCHit; + vplinks.applyToAllLinks( + [&channelIDForMCHit, &mchits]( unsigned int channelID, unsigned int mcHitKey, float weight ) { + channelIDForMCHit[mchits[mcHitKey]].emplace_back( WeightedHit{channelID, weight} ); + } ); // fill for all tracks for ( auto const mcp : mcparts ) { @@ -184,46 +239,117 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, // get MCHits for the candidate auto bhits = getMCHits( mcps, mchits ); auto first_mchit = bhits.size() ? bhits.front()->midPoint() : Gaudi::XYZPoint( 0., 0., 0. ); - // first link MCParticle to daughters (pions) of composite - auto composite = matchToComposites( mcps, composites, tracklinks, mcparts ); - auto has_comp = composite.has_value(); + // look at long tracks and see if pions match to all (as in 'reconstructed') + auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 ); + // link MCParticle to daughters (pions) of composite + auto composite = + m_composites.exist() ? matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ) : std::nullopt; + auto has_comp = composite.has_value(); // get heavy flavour track - auto itrack = has_comp ? getHeavyFlavourTrack( tracks, composite.value()->key() ) : std::nullopt; + auto const* tracks = m_tracks.getIfExists(); + auto itrack = ( has_comp && tracks ) ? getHeavyFlavourTrack( *tracks, composite.value()->key() ) : std::nullopt; auto has_track = has_comp ? itrack.has_value() : false; - auto has_hits = has_track ? (int)tracks.scalar()[itrack.value()].nHits() > 0 : false; + auto has_hits = has_track ? (int)tracks->scalar()[itrack.value()].nHits() > 0 : false; + int pv_key = has_track ? (int)tracks->scalar()[itrack.value()].primaryVertexIndex().cast() : -1; //// fill ntuple with info std::string base = "MCP_"; // mc truth auto sc = tuple->column( base + "TRUEID", mcp->particleID().pid() ); sc &= tuple->column( base + "P", mcp->momentum().P() ); + sc &= tuple->column( base + "PT", mcp->momentum().Pt() ); + sc &= tuple->column( base + "eta", mcp->momentum().Eta() ); + sc &= tuple->column( base + "phi", mcp->momentum().Phi() ); + sc &= tuple->column( base + "tx", mcp->momentum().Px() / mcp->momentum().Pz() ); + sc &= tuple->column( base + "ty", mcp->momentum().Py() / mcp->momentum().Pz() ); + auto ovtx = mcps.B->originVertex()->position(); + sc &= tuple->column( base + "OriginVertex_x", ovtx.x() ); + sc &= tuple->column( base + "OriginVertex_y", ovtx.y() ); + sc &= tuple->column( base + "OriginVertex_z", ovtx.z() ); + auto bvtx = mcps.B->endVertices().size() ? mcps.B->endVertices().at( mcps.B->endVertices().size() - 1 )->position() + : Gaudi::XYZPoint( 0., 0., -100 * Gaudi::Units::m ); + sc &= tuple->column( base + "EndVertex_x", bvtx.x() ); + sc &= tuple->column( base + "EndVertex_y", bvtx.y() ); + sc &= tuple->column( base + "EndVertex_z", bvtx.z() ); + auto tvtx = mcps.tau->endVertices().size() + ? mcps.tau->endVertices().at( mcps.tau->endVertices().size() - 1 )->position() + : Gaudi::XYZPoint( 0., 0., -110 * Gaudi::Units::m ); + sc &= tuple->column( base + "Tau_EndVertex_x", tvtx.x() ); + sc &= tuple->column( base + "Tau_EndVertex_y", tvtx.y() ); + sc &= tuple->column( base + "Tau_EndVertex_z", tvtx.z() ); + auto pionsmom = mcps.pions[0]->momentum() + mcps.pions[1]->momentum() + mcps.pions[2]->momentum(); + sc &= tuple->column( base + "Pions_P", pionsmom.P() ); + sc &= tuple->column( base + "Pions_PT", pionsmom.Pt() ); + sc &= tuple->column( base + "Pions_Mass", pionsmom.mass() ); + for ( int i = 0; i < (int)mcps.pions.size(); i++ ) { + std::string pbase = "Pion_" + std::to_string( i ) + "_"; + auto pion = mcps.pions[i]; + sc &= tuple->column( base + pbase + "TRUEID", pion->particleID().pid() ); + sc &= tuple->column( base + pbase + "P", pion->momentum().P() ); + sc &= tuple->column( base + pbase + "PT", pion->momentum().Pt() ); + sc &= tuple->column( base + pbase + "has_Track", hasTrack( pion, tracklinks ) ); + } // mc hits sc &= tuple->column( base + "has_VP_MCHits", (bool)bhits.size() ); sc &= tuple->column( base + "VP_first_MCHit_x", first_mchit.x() ); sc &= tuple->column( base + "VP_first_MCHit_y", first_mchit.y() ); sc &= tuple->column( base + "VP_first_MCHit_z", first_mchit.z() ); + // associated VP hits + std::vector<WeightedHit> hits_for_mchit = + bhits.size() ? channelIDForMCHit[bhits.front()] : std::vector<WeightedHit>(); + if ( hits_for_mchit.size() > 1 ) std::sort( hits_for_mchit.begin(), hits_for_mchit.end() ); + sc &= tuple->column( base + "VP_first_MCHit_has_VPHit", hits_for_mchit.size() > 0 ); + sc &= + tuple->column( base + "VP_first_MCHit_VPHit_ChannelID", hits_for_mchit.size() ? hits_for_mchit.front().id : 0 ); // composite info + sc &= tuple->column( base + "reconstructed", reconstructed ); sc &= tuple->column( base + "has_Composite", has_comp ); std::string compbase = base + "Composite_"; auto endvtx = has_comp ? composite.value()->endVertex()->position() : Gaudi::XYZPoint( 0., 0., 0. ); + float vtxchi2 = has_comp ? composite.value()->endVertex()->chi2() : -1.f; sc &= tuple->column( compbase + "EndVertex_x", endvtx.x() ); sc &= tuple->column( compbase + "EndVertex_y", endvtx.y() ); sc &= tuple->column( compbase + "EndVertex_z", endvtx.z() ); + sc &= tuple->column( compbase + "EndVertex_chi2", vtxchi2 ); auto compmom = has_comp ? composite.value()->momentum() : Gaudi::LorentzVector( 0., 0., 0., -1 * Gaudi::Units::GeV ); sc &= tuple->column( compbase + "P", compmom.R() ); sc &= tuple->column( compbase + "PT", compmom.Pt() ); sc &= tuple->column( compbase + "eta", compmom.Eta() ); sc &= tuple->column( compbase + "phi", compmom.Phi() ); + sc &= tuple->column( compbase + "Mass", compmom.mass() ); + // composite daughters info + using Functors::Common::ImpactParameter; + using Functors::Common::Position; + using Functors::Common::ToLinAlg; + using Sel::getBestPV; + using Sel::Utils::impactParameterChi2; + for ( int i = 0; i < (int)mcps.pions.size(); i++ ) { + std::string pbase = "Pion_" + std::to_string( i ) + "_"; + auto pion = mcps.recpions[i]; + float ip = pion ? ImpactParameter{}( + ToLinAlg{}( Position{}( has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ) ), pion ) + .cast() + : -1.f; + float ipchi2 = + pion ? impactParameterChi2( pion, has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ).cast() : -1.f; + sc &= tuple->column( compbase + pbase + "match", (bool)pion ); + sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. ); + sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. ); + sc &= tuple->column( compbase + pbase + "IP", ip ); + sc &= tuple->column( compbase + pbase + "IPChi2", ipchi2 ); + sc &= tuple->column( compbase + pbase + "PIDK", + pion ? pion->proto()->info( LHCb::ProtoParticle::additionalInfo::CombDLLk, 0. ) : 0. ); + } // track info std::string hfbase = base + "HeavyFlavourTracking_"; - sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks.scalar()[itrack.value()].nHits() : -1 ); + sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks->scalar()[itrack.value()].nHits() : -1 ); // track mc matching info bool fh_from_mcps = false; int fh_trueid = 0; if ( has_hits ) { float max_weight = 0.f; LHCb::MCHit const* mchit = nullptr; - unsigned int hid = tracks.scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID(); + unsigned int hid = tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID(); vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) { if ( weight > max_weight ) { max_weight = weight; @@ -236,13 +362,40 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, } } sc &= tuple->column( hfbase + "Track_first_hit_x", - has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. ); + has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. ); sc &= tuple->column( hfbase + "Track_first_hit_y", - has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. ); + has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. ); sc &= tuple->column( hfbase + "Track_first_hit_z", - has_hits ? tracks.scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. ); + has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. ); + sc &= tuple->column( hfbase + "Track_first_hit_LHCbID", + has_hits ? (unsigned int)tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID() : 0 ); sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps ); sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid ); + // more advanced variables for analysis + auto truedir = mcps.B->momentum().Vect().Unit(); + auto mcorr_true = has_comp ? getMcorr( composite.value(), truedir ) : 0.f; + float mcorr_fh = 0.f; + float mcorr_ps = 0.f; + auto slopes = Gaudi::XYZVector( 0., 0., 1. ); + auto pvpos = Gaudi::XYZPoint( 0., 0., -10 * Gaudi::Units::m ); + if ( has_track ) { + auto track = tracks->scalar()[itrack.value()]; + auto state = track.originState(); + slopes.SetXYZ( state.tx().cast(), state.ty().cast(), 1. ); + auto fhdir = slopes.Unit(); + pvpos.SetXYZ( state.x().cast(), state.y().cast(), state.z().cast() ); + auto pvsv = ( composite.value()->endVertex()->position() - pvpos ).Unit(); + mcorr_fh = getMcorr( composite.value(), fhdir ); + mcorr_ps = getMcorr( composite.value(), pvsv ); + } + sc &= tuple->column( hfbase + "tx", slopes.x() ); + sc &= tuple->column( hfbase + "ty", slopes.y() ); + sc &= tuple->column( hfbase + "true_dir_mcorr", mcorr_true ); + sc &= tuple->column( hfbase + "first_hit_mcorr", mcorr_fh ); + sc &= tuple->column( hfbase + "pvsv_mcorr", mcorr_ps ); + sc &= tuple->column( hfbase + "PV_x", pvpos.x() ); + sc &= tuple->column( hfbase + "PV_y", pvpos.y() ); + sc &= tuple->column( hfbase + "PV_z", pvpos.z() ); // write result sc.andThen( [&] { return tuple->write(); } ).ignore(); } -- GitLab From f6ce5ae8e031b57b16de2b3585014a4872785382 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <mveghel@cern.ch> Date: Wed, 9 Nov 2022 18:49:39 +0100 Subject: [PATCH 07/14] added btracking functor --- Phys/FunctorCore/include/Functors/Composite.h | 39 +++++++++++++++++++ .../include/Functors/JIT_includes.h | 1 + Phys/FunctorCore/python/Functors/__init__.py | 20 ++++++++++ .../tests/src/InstantiateFunctors.cpp | 6 +++ .../PrKernel/PrVeloHeavyFlavourUtils.h | 27 ++++++++++++- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 11 ++---- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 16 ++++---- 7 files changed, 103 insertions(+), 17 deletions(-) diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h index 284ba53cc60..cde1e593fb8 100755 --- a/Phys/FunctorCore/include/Functors/Composite.h +++ b/Phys/FunctorCore/include/Functors/Composite.h @@ -321,4 +321,43 @@ namespace Functors::Composite { } }; + namespace BTracking { + /** @brief Calculate the corrected mass using heavy flavour tracking info. */ + struct CorrectedMass : public Function { + static constexpr auto name() { return "BTracking::CorrectedMass"; } + + template <typename HeavyFlavourTracks, typename Composite> + auto operator()( HeavyFlavourTracks const& tracks, Composite const& composite ) const { + + using FType = SIMDWrapper::scalar::float_v; + using Sel::Utils::mass2; + using Sel::Utils::threeMomentum; + using std::sqrt; + auto const& comp = Sel::Utils::deref_if_ptr( composite ); + + // first get the right heavy flavour track + auto get_track = [&tracks]( int index ) -> std::optional<int> { + for ( auto track : tracks.scalar() ) { + if ( track.compositeIndex() == index ) return track.indices().cast(); + } + return std::nullopt; + }; + auto itrack = get_track( composite.key() ); + if ( !itrack.has_value() ) return FType{0}; + + // direction using first hit (or end vertex if not available, as saved in heavy flavour track) + auto const state = tracks.scalar()[itrack.value()].originState(); + auto const d = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.}; + + // Get the pT variable that we need + auto const mom = threeMomentum( comp ); + auto const perp = mom - d * ( dot( mom, d ) / d.mag2() ); + auto const pt = perp.mag(); + + // Calculate the corrected mass + return pt + sqrt( mass2( comp ) + pt * pt ); + } + }; + } // namespace BTracking + } // namespace Functors::Composite diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h index a83d656b988..8e3554e4d50 100644 --- a/Phys/FunctorCore/include/Functors/JIT_includes.h +++ b/Phys/FunctorCore/include/Functors/JIT_includes.h @@ -65,5 +65,6 @@ #include "SelKernel/VertexRelation.h" // PrKernel #include "PrKernel/PrSelection.h" +#include "PrKernel/PrVeloHeavyFlavourUtils.h" // RelTables #include "Event/TableView.h" diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py index a542a25e872..63e58f5cebc 100644 --- a/Phys/FunctorCore/python/Functors/__init__.py +++ b/Phys/FunctorCore/python/Functors/__init__.py @@ -1372,6 +1372,14 @@ MVA technique to establish pion mass hypothesis which keep in consideration the General information on PID could be found at the following `TWiki page <https://twiki.cern.ch/twiki/bin/view/LHCbPhysics/ChargedPID>`_ """) +_BPVCORRM = Functor( + '_BPVCORRM', 'Composite::CorrectedMass', + 'Compute the corrected mass of the composite using the associated [primary] vertex.' +) +_BTRACKING_BPVCORRM = Functor( + '_BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass', + 'Compute the corrected mass of the composite using the associated heavy flavour track.' +) def BPVCORRM(Vertices: DataHandle): """ Compute the corrected mass of the composite using the :py:func:`~BPV`. @@ -1405,6 +1413,18 @@ def BPVCORRMERR(Vertices: DataHandle): return _BPVCORRMERR.bind(BPV(Vertices), FORWARDARGS) +def BTRACKING_BPVCORRM(HeavyFlavourTracks: DataHandle): + """ Compute the corrected mass of the composite using the associated heavy flavour track. + + Functor's call operator expects a composite like object + + Args: + Vertices: DataHandle of the primary vertices + + """ + return _BTRACKING_BPVCORRM.bind(TES(HeavyFlavourTracks), FORWARDARGS) + + VTX_FDCHI2 = Functor( 'VTX_FDCHI2', 'Composite::FlightDistanceChi2ToVertex', '''Return the flight distance chi2 w.r.t. the given vertex.''') diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp index 7b4b936f371..7abeacb2a0e 100644 --- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp +++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp @@ -208,6 +208,12 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC // FIXME shouldn't these work? // using cmerr6 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHCb::RecVertex, LHCb::Particle const&>; +// +// BTracking::CorrectedMass( heavyflavourtracks, composite ) +// +using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass, + LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks, LHCb::Particle const&>; + // // MassWithHypotheses( vertices, composite ) // diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h index 941c308233b..eb5c7eff50b 100644 --- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h +++ b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h @@ -24,7 +24,7 @@ namespace LHCb::Pr::Velo { namespace Tag { // hit info struct LHCbID : Event::lhcbid_field {}; - struct HitPosition : Event::vec3_field {}; + struct HitPosition : Event::Vec_field<3> {}; struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {}; // crossed sensor info @@ -55,10 +55,22 @@ namespace LHCb::Pr::Velo { using Composite = LHCb::Particle; using PVs = LHCb::Event::PV::PrimaryVertexContainer; - struct HeavyFlavourTracks : Tag::Track_t<HeavyFlavourTracks> { + class HeavyFlavourTracks : public Tag::Track_t<HeavyFlavourTracks> { + public: using base_t = typename Tag::Track_t<HeavyFlavourTracks>; + using base_t::allocator_type; using base_t::base_t; + HeavyFlavourTracks( UniqueIDGenerator const& unique_id_gen, + Zipping::ZipFamilyNumber zip_identifier = Zipping::generateZipIdentifier(), + allocator_type alloc = {} ) + : base_t{std::move( zip_identifier ), std::move( alloc )}, m_unique_id_gen_tag{unique_id_gen.tag()} {} + + // special constructor for zipping + HeavyFlavourTracks( Zipping::ZipFamilyNumber zip_identifier, HeavyFlavourTracks const& other ) + : base_t{std::move( zip_identifier ), other} {} + + // proxy template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType> struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> { using base_t = typename Event::Proxy<simd, behaviour, ContainerType>; @@ -89,6 +101,17 @@ namespace LHCb::Pr::Velo { template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> using proxy_type = TrackProxy<simd, behaviour, ContainerType>; + + template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour> + using reference = LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, HeavyFlavourTracks>>; + template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour> + using const_reference = const LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, const HeavyFlavourTracks>>; + + auto const& unique_id_gen_tag() const { return m_unique_id_gen_tag; } + + private: + // keep the identifier of the generator used to build this container + boost::uuids::uuid m_unique_id_gen_tag; }; } // namespace v1 diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index a7caf773073..94d3418ee9c 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -34,8 +34,8 @@ namespace { // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA) using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks; - using Composites = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composites; using Composite = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite; + using Composites = Composite::Selection; using PVs = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs; using MCParticles = LHCb::MCParticles; using MCHits = LHCb::MCHits; @@ -326,12 +326,9 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, for ( int i = 0; i < (int)mcps.pions.size(); i++ ) { std::string pbase = "Pion_" + std::to_string( i ) + "_"; auto pion = mcps.recpions[i]; - float ip = pion ? ImpactParameter{}( - ToLinAlg{}( Position{}( has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ) ), pion ) - .cast() - : -1.f; - float ipchi2 = - pion ? impactParameterChi2( pion, has_track ? pvs[pv_key] : getBestPV( *composite, pvs ) ).cast() : -1.f; + float ip = + ( pion && has_track ) ? ImpactParameter{}( ToLinAlg{}( Position{}( pvs[pv_key] ) ), pion ).cast() : -1.f; + float ipchi2 = ( pion && has_track ) ? impactParameterChi2( pion, pvs[pv_key] ).cast() : -1.f; sc &= tuple->column( compbase + pbase + "match", (bool)pion ); sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. ); sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. ); diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index 5fab6af5c4b..24c75c8b8fe 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -33,7 +33,7 @@ namespace LHCb::Pr::Velo { // miscellaneous namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag; using FType = SIMDWrapper::scalar::float_v; - using Position = LHCb::LinAlg::Vec3<FType>; + using Position = LHCb::LinAlg::Vec<FType, 3>; using SensorID = LHCb::Detector::VPChannelID::SensorID; using ChannelID = LHCb::Detector::VPChannelID; @@ -82,7 +82,7 @@ namespace LHCb::Pr::Velo { ChannelID id; // sorting in z first, then distance to search window friend bool operator<( HitInfo const& lhs, HitInfo const& rhs ) { - return std::pair{lhs.position.z.cast(), lhs.d.cast()} < std::pair{rhs.position.z.cast(), rhs.d.cast()}; + return std::pair{lhs.position.z().cast(), lhs.d.cast()} < std::pair{rhs.position.z().cast(), rhs.d.cast()}; } // equal just by unique channel id friend bool operator==( HitInfo const& lhs, HitInfo const& rhs ) { return lhs.id == rhs.id; } @@ -119,7 +119,7 @@ namespace LHCb::Pr::Velo { // its associated module (to not be affected by borders between sensors) auto crossed_module = sensor.module(); // go to sensor location - pos = pos + ( sensor.z - pos.z.cast() ) * slopes; + pos = pos + ( sensor.z - pos.z().cast() ) * slopes; // look if this sensor (and associated module) has hits and if they are in search window for ( auto hit : hits.scalar() ) { auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); @@ -127,9 +127,9 @@ namespace LHCb::Pr::Velo { auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); // make sure they are in the same plane up to some tolerance // (not all sensors in a module are exactly) - if ( abs( hitpos.vec3().z.cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) { - auto d = ( hitpos.vec3() - pos ).mag(); - if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos.vec3(), id} ); + if ( abs( hitpos.z().cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) { + auto d = ( hitpos - pos ).mag(); + if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos, id} ); } } } @@ -214,7 +214,7 @@ namespace LHCb::Pr::Velo { // primary vertex using Sel::Utils::endVertexPos; - auto const& PV = Sel::getBestPV( *composite, pvs ); + auto const& PV = Sel::getBestPV( Sel::Utils::deref_if_ptr( composite ), pvs ); auto const& pv = endVertexPos( PV ); track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() ); @@ -245,7 +245,7 @@ namespace LHCb::Pr::Velo { auto direction = found_hits.size() ? found_hits.front().position - pv : fd; auto origin_state = track.field<TrackTag::OriginState>(); origin_state.setPosition( pv ); - origin_state.setDirection( direction.x / direction.z, direction.y / direction.z ); + origin_state.setDirection( direction.x() / direction.z(), direction.y() / direction.z() ); origin_state.setQOverP( 1. / composite->p() ); // fitting -- GitLab From 77bcba1e1c4252553254ff5b192080f3f9bcb550 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Tue, 25 Jul 2023 19:03:06 +0200 Subject: [PATCH 08/14] move to v1 event model --- Phys/FunctorCore/include/Functors/Composite.h | 17 +-- .../include/Functors/JIT_includes.h | 1 + .../tests/src/InstantiateFunctors.cpp | 4 +- .../PrKernel/PrVeloHeavyFlavourUtils.h | 120 ------------------ .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 94 +++++++------- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 116 +++++++++-------- 6 files changed, 121 insertions(+), 231 deletions(-) delete mode 100644 Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h index cde1e593fb8..9a00ace18f7 100755 --- a/Phys/FunctorCore/include/Functors/Composite.h +++ b/Phys/FunctorCore/include/Functors/Composite.h @@ -326,8 +326,8 @@ namespace Functors::Composite { struct CorrectedMass : public Function { static constexpr auto name() { return "BTracking::CorrectedMass"; } - template <typename HeavyFlavourTracks, typename Composite> - auto operator()( HeavyFlavourTracks const& tracks, Composite const& composite ) const { + template <typename Composite2TrackRelations, typename Composite> + auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const { using FType = SIMDWrapper::scalar::float_v; using Sel::Utils::mass2; @@ -336,17 +336,12 @@ namespace Functors::Composite { auto const& comp = Sel::Utils::deref_if_ptr( composite ); // first get the right heavy flavour track - auto get_track = [&tracks]( int index ) -> std::optional<int> { - for ( auto track : tracks.scalar() ) { - if ( track.compositeIndex() == index ) return track.indices().cast(); - } - return std::nullopt; - }; - auto itrack = get_track( composite.key() ); - if ( !itrack.has_value() ) return FType{0}; + auto track_range = relations.relations( &comp ); + if ( track_range.empty() ) return FType{0}; + auto track = track_range.front().to(); // direction using first hit (or end vertex if not available, as saved in heavy flavour track) - auto const state = tracks.scalar()[itrack.value()].originState(); + auto const state = track->firstState(); auto const d = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.}; // Get the pT variable that we need diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h index 8e3554e4d50..b1510005c2d 100644 --- a/Phys/FunctorCore/include/Functors/JIT_includes.h +++ b/Phys/FunctorCore/include/Functors/JIT_includes.h @@ -68,3 +68,4 @@ #include "PrKernel/PrVeloHeavyFlavourUtils.h" // RelTables #include "Event/TableView.h" +#include "Relations/Relation2D.h" diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp index 7abeacb2a0e..6c08a94cff3 100644 --- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp +++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp @@ -209,10 +209,10 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC // using cmerr6 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHCb::RecVertex, LHCb::Particle const&>; // -// BTracking::CorrectedMass( heavyflavourtracks, composite ) +// BTracking::CorrectedMass( heavyflavourtrackrelations, composite ) // using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass, - LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks, LHCb::Particle const&>; + LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>; // // MassWithHypotheses( vertices, composite ) diff --git a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h b/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h deleted file mode 100644 index eb5c7eff50b..00000000000 --- a/Pr/PrKernel/include/PrKernel/PrVeloHeavyFlavourUtils.h +++ /dev/null @@ -1,120 +0,0 @@ -/*****************************************************************************\ -* (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 "Event/Particle.h" -#include "Event/Particle_v2.h" -#include "Event/PrTracksTag.h" -#include "Event/PrimaryVertices.h" -#include "Event/RelationTable.h" -#include "VPDet/DeVP.h" - -namespace LHCb::Pr::Velo { - - namespace HeavyFlavourTracking { - - // main heavy flavour track storage - namespace Tag { - // hit info - struct LHCbID : Event::lhcbid_field {}; - struct HitPosition : Event::Vec_field<3> {}; - struct Hits : Event::vector_field<Event::struct_field<LHCbID, HitPosition>> {}; - - // crossed sensor info - struct nSensors : Event::int_field {}; - - // geometry - struct OriginState : Event::state_field {}; - - // track quality info - struct Chi2 : Event::float_field {}; - - // (hopefully temporary) fix with SOACollection with keys/indices to non SoA objects - struct CompositeIndex : Event::int_field {}; - struct PrimaryVertexIndex : Event::int_field {}; - - template <typename T> - using Track_t = Event::SOACollection<T, CompositeIndex, PrimaryVertexIndex, Hits, nSensors, OriginState, Chi2>; - - /* - // main storage: relating track to associated composite and PV in full SoA event model - template <typename Composites, typename PVs> - using RelationTable2D_t = LHCb::Event::RelationTable2D<Composites, PVs, Hits, nSensors, OriginState, Chi2>; - */ - } // namespace Tag - - namespace v1 { - using Composites = LHCb::Particles; - using Composite = LHCb::Particle; - using PVs = LHCb::Event::PV::PrimaryVertexContainer; - - class HeavyFlavourTracks : public Tag::Track_t<HeavyFlavourTracks> { - public: - using base_t = typename Tag::Track_t<HeavyFlavourTracks>; - using base_t::allocator_type; - using base_t::base_t; - - HeavyFlavourTracks( UniqueIDGenerator const& unique_id_gen, - Zipping::ZipFamilyNumber zip_identifier = Zipping::generateZipIdentifier(), - allocator_type alloc = {} ) - : base_t{std::move( zip_identifier ), std::move( alloc )}, m_unique_id_gen_tag{unique_id_gen.tag()} {} - - // special constructor for zipping - HeavyFlavourTracks( Zipping::ZipFamilyNumber zip_identifier, HeavyFlavourTracks const& other ) - : base_t{std::move( zip_identifier ), other} {} - - // proxy - template <SIMDWrapper::InstructionSet simd, ProxyBehaviour behaviour, typename ContainerType> - struct TrackProxy : Event::Proxy<simd, behaviour, ContainerType> { - using base_t = typename Event::Proxy<simd, behaviour, ContainerType>; - using base_t::Proxy; - - // composite / pv info - [[nodiscard]] auto compositeIndex() const { return this->template get<Tag::CompositeIndex>(); } - [[nodiscard]] auto primaryVertexIndex() const { return this->template get<Tag::PrimaryVertexIndex>(); } - - // hit info - [[nodiscard]] auto nHits() const { return this->template field<Tag::Hits>().size(); } - [[nodiscard]] auto hitID( std::size_t i ) const { - return this->template field<Tag::Hits>()[i].template get<Tag::LHCbID>(); - } - [[nodiscard]] auto hitPosition( std::size_t i ) const { - return this->template field<Tag::Hits>()[i].template get<Tag::HitPosition>(); - } - - // crossed sensor info - [[nodiscard]] auto nSensors() const { return this->template get<Tag::nSensors>(); } - - // geometry - [[nodiscard]] auto originState() const { return this->template get<Tag::OriginState>(); } - - // track quality info - [[nodiscard]] auto chi2() const { return this->template get<Tag::Chi2>(); } - }; - - template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour, typename ContainerType> - using proxy_type = TrackProxy<simd, behaviour, ContainerType>; - - template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour> - using reference = LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, HeavyFlavourTracks>>; - template <SIMDWrapper::InstructionSet simd, LHCb::Pr::ProxyBehaviour behaviour> - using const_reference = const LHCb::Event::ZipProxy<TrackProxy<simd, behaviour, const HeavyFlavourTracks>>; - - auto const& unique_id_gen_tag() const { return m_unique_id_gen_tag; } - - private: - // keep the identifier of the generator used to build this container - boost::uuids::uuid m_unique_id_gen_tag; - }; - } // namespace v1 - - } // namespace HeavyFlavourTracking - -} // namespace LHCb::Pr::Velo diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index 94d3418ee9c..711dbf61ead 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -13,7 +13,7 @@ #include "Event/MCHit.h" #include "Event/MCParticle.h" #include "Event/RawEvent.h" -#include "Event/Track_v3.h" +#include "Event/Track.h" #include "Event/VPDigit.h" #include "Event/VPFullCluster.h" #include "Functors/Common.h" @@ -21,7 +21,7 @@ #include "GaudiAlg/GaudiTupleAlg.h" #include "Linker/LinkedFrom.h" #include "Linker/LinkedTo.h" -#include "PrKernel/PrVeloHeavyFlavourUtils.h" +#include "Relations/Relation2D.h" #include "SelKernel/Utilities.h" #include "SelKernel/VertexRelation.h" @@ -32,14 +32,16 @@ */ namespace { - // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA) - using HeavyFlavourTracks = LHCb::Pr::Velo::HeavyFlavourTracking::v1::HeavyFlavourTracks; - using Composite = LHCb::Pr::Velo::HeavyFlavourTracking::v1::Composite; - using Composites = Composite::Selection; - using PVs = LHCb::Pr::Velo::HeavyFlavourTracking::v1::PVs; - using MCParticles = LHCb::MCParticles; - using MCHits = LHCb::MCHits; - using Links = LHCb::LinksByKey; + // in/outputs + using Composite = LHCb::Particle; + using Composites = LHCb::Particle::Selection; + using PVs = LHCb::Event::PV::PrimaryVertexContainer; + using Track = LHCb::Event::v1::Track; + using Tracks = Track::Container; + using P2TRelations = LHCb::Relation2D<LHCb::Particle, Track>; + using MCParticles = LHCb::MCParticles; + using MCHits = LHCb::MCHits; + using Links = LHCb::LinksByKey; // helper classes struct DecayMCParticles { @@ -165,14 +167,6 @@ namespace { return hits; } - // get track - std::optional<int> getHeavyFlavourTrack( HeavyFlavourTracks const& tracks, int key ) { - for ( auto track : tracks.scalar() ) { - if ( track.compositeIndex() == key ) return track.indices().cast(); - } - return std::nullopt; - } - // composite observables float getMcorr( Composite const* composite, Gaudi::XYZVector const& direction ) { auto mom = composite->momentum().Vect(); @@ -205,8 +199,8 @@ public: private: // data that might not be there - DataObjectReadHandle<Composites> m_composites{this, "Composites", ""}; - DataObjectReadHandle<HeavyFlavourTracks> m_tracks{this, "HeavyFlavourTracks", ""}; + DataObjectReadHandle<Composites> m_composites{this, "Composites", ""}; + DataObjectReadHandle<P2TRelations> m_relations{this, "Composite2HeavyFlavourTrackRelations", ""}; // monitoring mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""}; @@ -242,15 +236,20 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, // look at long tracks and see if pions match to all (as in 'reconstructed') auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 ); // link MCParticle to daughters (pions) of composite - auto composite = - m_composites.exist() ? matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ) : std::nullopt; + std::optional<Composite const*> composite = std::nullopt; + if ( m_composites.exist() ) { + composite = matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ); + } auto has_comp = composite.has_value(); // get heavy flavour track - auto const* tracks = m_tracks.getIfExists(); - auto itrack = ( has_comp && tracks ) ? getHeavyFlavourTrack( *tracks, composite.value()->key() ) : std::nullopt; - auto has_track = has_comp ? itrack.has_value() : false; - auto has_hits = has_track ? (int)tracks->scalar()[itrack.value()].nHits() > 0 : false; - int pv_key = has_track ? (int)tracks->scalar()[itrack.value()].primaryVertexIndex().cast() : -1; + auto const* relations = m_relations.getIfExists(); + Track const* track = nullptr; + if ( has_comp && relations ) { + auto trackrange = relations->relations( composite.value() ); + if ( trackrange.size() ) track = trackrange.front(); + } + auto has_track = track ? true : false; + auto has_hits = has_track ? track->nHits() > 0 : false; //// fill ntuple with info std::string base = "MCP_"; // mc truth @@ -324,11 +323,15 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, using Sel::getBestPV; using Sel::Utils::impactParameterChi2; for ( int i = 0; i < (int)mcps.pions.size(); i++ ) { - std::string pbase = "Pion_" + std::to_string( i ) + "_"; - auto pion = mcps.recpions[i]; - float ip = - ( pion && has_track ) ? ImpactParameter{}( ToLinAlg{}( Position{}( pvs[pv_key] ) ), pion ).cast() : -1.f; - float ipchi2 = ( pion && has_track ) ? impactParameterChi2( pion, pvs[pv_key] ).cast() : -1.f; + std::string pbase = "Pion_" + std::to_string( i ) + "_"; + auto pion = mcps.recpions[i]; + float ip = -1.f; + float ipchi2 = -1.f; + if ( pion && has_track ) { + auto pv = getBestPV( pion, pvs ); + ip = ImpactParameter{}( ToLinAlg{}( Position{}( pv ) ), pion ).cast(); + ipchi2 = impactParameterChi2( pion, pv ).cast(); + } sc &= tuple->column( compbase + pbase + "match", (bool)pion ); sc &= tuple->column( compbase + pbase + "P", pion ? pion->momentum().R() : 0. ); sc &= tuple->column( compbase + pbase + "PT", pion ? pion->momentum().Pt() : 0. ); @@ -339,14 +342,16 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, } // track info std::string hfbase = base + "HeavyFlavourTracking_"; - sc &= tuple->column( hfbase + "Track_nHits", has_track ? (int)tracks->scalar()[itrack.value()].nHits() : -1 ); + sc &= tuple->column( hfbase + "Track_nSensors", + has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::VeloNSensors, 0 ) : -1. ); + sc &= tuple->column( hfbase + "Track_nHits", has_track ? track->nHits() : -1 ); // track mc matching info bool fh_from_mcps = false; int fh_trueid = 0; if ( has_hits ) { float max_weight = 0.f; LHCb::MCHit const* mchit = nullptr; - unsigned int hid = tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID(); + unsigned int hid = track->lhcbIDs().front().vpID(); vplinks.applyToLinks( hid, [&max_weight, &mchits, &mchit]( unsigned int, unsigned int mcHitKey, float weight ) { if ( weight > max_weight ) { max_weight = weight; @@ -358,14 +363,14 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, fh_from_mcps = ( mchit->mcParticle() == mcps.B ) || ( mchit->mcParticle() == mcps.tau ); } } - sc &= tuple->column( hfbase + "Track_first_hit_x", - has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).x().cast() : 0. ); - sc &= tuple->column( hfbase + "Track_first_hit_y", - has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).y().cast() : 0. ); - sc &= tuple->column( hfbase + "Track_first_hit_z", - has_hits ? tracks->scalar()[itrack.value()].hitPosition( 0 ).z().cast() : 0. ); + decltype( track->position() ) fh_pos( 0., 0., 0. ); + LHCb::State const* fh_state = track ? track->stateAt( LHCb::State::Location::FirstMeasurement ) : nullptr; + if ( fh_state ) fh_pos = fh_state->position(); + sc &= tuple->column( hfbase + "Track_first_hit_x", fh_pos.x() ); + sc &= tuple->column( hfbase + "Track_first_hit_y", fh_pos.y() ); + sc &= tuple->column( hfbase + "Track_first_hit_z", fh_pos.z() ); sc &= tuple->column( hfbase + "Track_first_hit_LHCbID", - has_hits ? (unsigned int)tracks->scalar()[itrack.value()].hitID( 0 ).LHCbID().vpID() : 0 ); + has_hits ? (unsigned int)track->lhcbIDs().front().vpID() : 0 ); sc &= tuple->column( hfbase + "Track_first_hit_from_MCPs", fh_from_mcps ); sc &= tuple->column( hfbase + "Track_first_hit_MC_TRUEID", fh_trueid ); // more advanced variables for analysis @@ -376,11 +381,10 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, auto slopes = Gaudi::XYZVector( 0., 0., 1. ); auto pvpos = Gaudi::XYZPoint( 0., 0., -10 * Gaudi::Units::m ); if ( has_track ) { - auto track = tracks->scalar()[itrack.value()]; - auto state = track.originState(); - slopes.SetXYZ( state.tx().cast(), state.ty().cast(), 1. ); + auto state = track->firstState(); + slopes.SetXYZ( state.tx(), state.ty(), 1. ); auto fhdir = slopes.Unit(); - pvpos.SetXYZ( state.x().cast(), state.y().cast(), state.z().cast() ); + pvpos.SetXYZ( state.x(), state.y(), state.z() ); auto pvsv = ( composite.value()->endVertex()->position() - pvpos ).Unit(); mcorr_fh = getMcorr( composite.value(), fhdir ); mcorr_ps = getMcorr( composite.value(), pvsv ); diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index 24c75c8b8fe..ce3f506d4ec 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -12,30 +12,33 @@ #include "DetDesc/GenericConditionAccessorHolder.h" #include "DetDesc/IConditionDerivationMgr.h" #include "Event/PrHits.h" +#include "Event/Track.h" #include "GaudiAlg/Transformer.h" #include "Kernel/LHCbID.h" #include "LHCbMath/MatVec.h" -#include "PrKernel/PrVeloHeavyFlavourUtils.h" +#include "Relations/Relation2D.h" #include "SelKernel/Utilities.h" #include "SelKernel/VertexRelation.h" +#include "VPDet/DeVP.h" #include "VPDet/VPDetPaths.h" -#include "VeloKalmanHelpers.h" namespace LHCb::Pr::Velo { namespace { - // in/outputs (TODO change to RelationTable (v3) when PVs and Composites are all SoA) - using namespace LHCb::Pr::Velo::HeavyFlavourTracking::v1; - using VeloHits = LHCb::Pr::VP::Hits; - using Output = HeavyFlavourTracks; + // in/outputs + using VeloHits = LHCb::Pr::VP::Hits; + using PVs = LHCb::Event::PV::PrimaryVertexContainer; + using Track = LHCb::Event::v1::Track; + using Tracks = Track::Container; + using Relations = LHCb::Relation2D<LHCb::Particle, Track>; + using Output = std::tuple<Tracks, Relations>; // miscellaneous - namespace TrackTag = LHCb::Pr::Velo::HeavyFlavourTracking::Tag; - using FType = SIMDWrapper::scalar::float_v; - using Position = LHCb::LinAlg::Vec<FType, 3>; - using SensorID = LHCb::Detector::VPChannelID::SensorID; - using ChannelID = LHCb::Detector::VPChannelID; + using FType = SIMDWrapper::scalar::float_v; + using Position = LHCb::LinAlg::Vec<FType, 3>; + using SensorID = LHCb::Detector::VPChannelID::SensorID; + using ChannelID = LHCb::Detector::VPChannelID; // to store relevant sensor information struct SensorInfo { @@ -150,9 +153,9 @@ namespace LHCb::Pr::Velo { */ class HeavyFlavourTrackFinder - : public Gaudi::Functional::Transformer<Output( Composites const&, PVs const&, VeloHits const&, DeVP const&, - Sensors const& ), - LHCb::DetDesc::usesConditions<DeVP, Sensors>> { + : public Gaudi::Functional::MultiTransformer<Output( LHCb::Particle::Range const&, PVs const&, VeloHits const&, + DeVP const&, Sensors const& ), + LHCb::DetDesc::usesConditions<DeVP, Sensors>> { public: // standard constructor HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator ); @@ -161,7 +164,8 @@ namespace LHCb::Pr::Velo { StatusCode initialize() override; // main function/operator - Output operator()( Composites const&, PVs const&, VeloHits const&, DeVP const&, Sensors const& ) const override; + Output operator()( LHCb::Particle::Range const&, PVs const&, VeloHits const&, DeVP const&, + Sensors const& ) const override; private: // properties @@ -179,44 +183,39 @@ namespace LHCb::Pr::Velo { // constructor HeavyFlavourTrackFinder::HeavyFlavourTrackFinder( const std::string& name, ISvcLocator* pSvcLocator ) - : Transformer( name, pSvcLocator, - {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ), - KeyValue( "DeVP", LHCb::Det::VP::det_path ), -#ifdef USE_DD4HEP - KeyValue( "Sensors", {"/world:AlgorithmSpecific-" + name + "-sensorinfos"} )}, -#else - KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )}, -#endif - - {KeyValue{"OutputLocation", ""}} ) { - } + : MultiTransformer( name, pSvcLocator, + {KeyValue( "Composites", "" ), KeyValue( "PVs", "" ), KeyValue( "Hits", "" ), + KeyValue( "DeVP", LHCb::Det::VP::det_path ), + KeyValue( "Sensors", {"AlgorithmSpecific-" + name + "-sensorinfos"} )}, + {KeyValue{"OutputTracks", ""}, KeyValue{"OutputRelations", ""}} ) {} // initialization StatusCode LHCb::Pr::Velo::HeavyFlavourTrackFinder::initialize() { - auto sc = Transformer::initialize().andThen( [&] { + auto sc = MultiTransformer::initialize().andThen( [&] { addConditionDerivation<Sensors( DeVP const& )>( {LHCb::Det::VP::det_path}, inputLocation<Sensors>() ); } ); return sc; } // main execution - Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( Composites const& composites, PVs const& pvs, + Output LHCb::Pr::Velo::HeavyFlavourTrackFinder::operator()( LHCb::Particle::Range const& composites, PVs const& pvs, VeloHits const& velohits, DeVP const& velo, Sensors const& sensors ) const { // output - Output output{}; - output.reserve( composites.size() ); + auto result = Output{}; + auto& [tracks, relations] = result; + tracks.reserve( composites.size() ); + relations.reserve( composites.size() ); - for ( auto const composite : composites ) { + for ( LHCb::Particle const* composite : composites ) { // set output - auto track = output.emplace_back<SIMDWrapper::InstructionSet::Scalar>(); - track.field<TrackTag::CompositeIndex>().set( composite->key() ); + using namespace LHCb::Event::Enum::Track; + auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs ); // primary vertex using Sel::Utils::endVertexPos; - auto const& PV = Sel::getBestPV( Sel::Utils::deref_if_ptr( composite ), pvs ); - auto const& pv = endVertexPos( PV ); - track.field<TrackTag::PrimaryVertexIndex>().set( PV.key() ); + auto const& pv_obj = Sel::getBestPV( *composite, pvs ); + auto const& pv = endVertexPos( pv_obj ); // secondary vertex auto const& end = endVertexPos( *composite ); @@ -225,39 +224,50 @@ namespace LHCb::Pr::Velo { // crossed sensors in PV-SV line auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes ); - track.field<TrackTag::nSensors>().set( crossedSensors.size() ); + track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() ); // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window) auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); + auto nhits = found_hits.size(); // store sorted hits - auto saved_hits = track.field<TrackTag::Hits>(); - auto nhits = found_hits.size(); - saved_hits.resize( nhits ); - for ( long unsigned int i = 0; i < nhits; i++ ) { - auto found_hit = found_hits[i]; - saved_hits[i].field<TrackTag::LHCbID>().set( LHCbID( found_hit.id ) ); - saved_hits[i].field<TrackTag::HitPosition>().set( found_hit.position ); - } + std::vector<LHCbID> found_ids; + found_ids.reserve( nhits ); + for ( auto const found_hit : found_hits ) { found_ids.push_back( LHCbID( found_hit.id ) ); } + track->setLhcbIDs( found_ids ); // determine mother particle direction based on closest hit // in first crossed sensor or using SV-PV if not available auto direction = found_hits.size() ? found_hits.front().position - pv : fd; - auto origin_state = track.field<TrackTag::OriginState>(); - origin_state.setPosition( pv ); - origin_state.setDirection( direction.x() / direction.z(), direction.y() / direction.z() ); - origin_state.setQOverP( 1. / composite->p() ); + auto origin_state = LHCb::State(); + origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(), + ( direction.y() / direction.z() ).cast(), composite->charge() / composite->p() ); + origin_state.setLocation( LHCb::State::Location::ClosestToBeam ); + track->addToStates( origin_state ); + + // in case of found hit, add appropriate state + if ( found_hits.size() ) { + auto fh = found_hits.front(); + auto fm_state = + LHCb::State( origin_state.stateVector(), fh.position.z().cast(), LHCb::State::Location::FirstMeasurement ); + fm_state.setX( fh.position.x().cast() ); + fm_state.setY( fh.position.y().cast() ); + track->addToStates( fm_state ); + } - // fitting - // TODO (using VeloKalmanHelpers) - track.field<TrackTag::Chi2>().set( -1.f ); + // save track + tracks.insert( track ); + relations.i_push( composite, track ); // statistics m_nInSensor += crossedSensors.size(); m_nHitCandidates += nhits; } - return output; + // mandatory sorting of relation table + relations.i_sort(); + + return result; } } // namespace LHCb::Pr::Velo -- GitLab From f11f4dda19f1d5789ee2de5c502e257ef8889242 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Thu, 27 Jul 2023 15:31:44 +0200 Subject: [PATCH 09/14] add tolerance to sensor lookup for search cylinder, not line --- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 11 ++++------- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 19 +++++++++++-------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index 711dbf61ead..fb878247a6c 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -34,10 +34,9 @@ namespace { // in/outputs using Composite = LHCb::Particle; - using Composites = LHCb::Particle::Selection; + using Composites = LHCb::Particle::Range; using PVs = LHCb::Event::PV::PrimaryVertexContainer; using Track = LHCb::Event::v1::Track; - using Tracks = Track::Container; using P2TRelations = LHCb::Relation2D<LHCb::Particle, Track>; using MCParticles = LHCb::MCParticles; using MCHits = LHCb::MCHits; @@ -236,11 +235,9 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, // look at long tracks and see if pions match to all (as in 'reconstructed') auto reconstructed = matchToTracks( mcps, tracklinks, mcparts, 0.7 ); // link MCParticle to daughters (pions) of composite - std::optional<Composite const*> composite = std::nullopt; - if ( m_composites.exist() ) { - composite = matchToComposites( mcps, *m_composites.get(), tracklinks, mcparts, 0.7 ); - } - auto has_comp = composite.has_value(); + auto composites = m_composites.get(); + auto composite = matchToComposites( mcps, composites, tracklinks, mcparts, 0.7 ); + auto has_comp = composite.has_value(); // get heavy flavour track auto const* relations = m_relations.getIfExists(); Track const* track = nullptr; diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index ce3f506d4ec..dd7dd0b21e3 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -95,15 +95,18 @@ namespace LHCb::Pr::Velo { // functions template <typename Position3, typename Vector3> SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start, - Position3 const& end, Vector3 const& slopes ) { + Position3 const& end, Vector3 const& slopes, float margin ) { SensorInfos crossedSensors; - auto pos = start; + auto pos = start; + auto margin_normalized = margin / sqrt( slopes.mag2() - 1 ); for ( auto sensor : sensors ) { if ( sensor.z < pos.Z() ) continue; if ( sensor.z > end.Z() ) break; - // if in range look if search window intersects with sensor - pos = pos + ( sensor.z - pos.Z() ) * slopes; - if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) ) + // if in range look if search window intersects with sensor (up to some margin deeper into VELO) + pos = pos + ( sensor.z - pos.Z() ) * slopes; + auto pos_with_margin = pos + margin_normalized * slopes; + if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) | + velo.sensor( sensor.id ).isInsideSensor( pos_with_margin.X().cast(), pos_with_margin.Y().cast() ) ) crossedSensors.push_back( {pos.Z().cast(), sensor.id} ); } return crossedSensors; @@ -222,11 +225,11 @@ namespace LHCb::Pr::Velo { auto const fd = end - pv; auto const slopes = fd / fd.Z(); - // crossed sensors in PV-SV line - auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes ); + // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits) + auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() ); track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() ); - // look for hits in relevant sensors (duplicated remoded and sorted in z and distance to search window) + // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window) auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); auto nhits = found_hits.size(); -- GitLab From 08c83c83f6dbb776dedb0e2cb78d8d0c7094600e Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Thu, 27 Jul 2023 17:37:42 +0200 Subject: [PATCH 10/14] adds Particle maker for containing heavy flavour tracks --- .../include/Functors/JIT_includes.h | 1 - Phys/ParticleMaker/CMakeLists.txt | 1 + .../ParticleWithHeavyFlavourTrackMaker.cpp | 104 ++++++++++++++++++ 3 files changed, 105 insertions(+), 1 deletion(-) create mode 100644 Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp diff --git a/Phys/FunctorCore/include/Functors/JIT_includes.h b/Phys/FunctorCore/include/Functors/JIT_includes.h index b1510005c2d..1b05e840f21 100644 --- a/Phys/FunctorCore/include/Functors/JIT_includes.h +++ b/Phys/FunctorCore/include/Functors/JIT_includes.h @@ -65,7 +65,6 @@ #include "SelKernel/VertexRelation.h" // PrKernel #include "PrKernel/PrSelection.h" -#include "PrKernel/PrVeloHeavyFlavourUtils.h" // RelTables #include "Event/TableView.h" #include "Relations/Relation2D.h" diff --git a/Phys/ParticleMaker/CMakeLists.txt b/Phys/ParticleMaker/CMakeLists.txt index 67420524060..4d39be7209d 100644 --- a/Phys/ParticleMaker/CMakeLists.txt +++ b/Phys/ParticleMaker/CMakeLists.txt @@ -23,6 +23,7 @@ gaudi_add_module(ParticleMaker src/ParticleMassDTFMonitor.cpp src/ParticleMassMonitor.cpp src/ParticleWithBremMaker.cpp + src/ParticleWithHeavyFlavourTrackMaker.cpp src/V0FromDstMaker.cpp src/ParticlesEmptyProducer.cpp LINK diff --git a/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp b/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp new file mode 100644 index 00000000000..3307e55dbc8 --- /dev/null +++ b/Phys/ParticleMaker/src/ParticleWithHeavyFlavourTrackMaker.cpp @@ -0,0 +1,104 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2023 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 "Event/Particle.h" +#include "Event/ProtoParticle.h" +#include "Event/Track.h" +#include "Kernel/IParticlePropertySvc.h" +#include "Kernel/ParticleProperty.h" +#include "LHCbAlgs/Transformer.h" +#include "Relations/Relation2D.h" + +namespace { + using Track = LHCb::Event::v1::Track; + using RelationTable = LHCb::Relation2D<LHCb::Particle, Track>; + using Output = std::tuple<LHCb::Particles, LHCb::ProtoParticles>; +} // namespace + +class ParticleWithHeavyFlavourTrackMaker + : public LHCb::Algorithm::MultiTransformer<Output( LHCb::Particle::Range const&, RelationTable const& )> { +public: + ParticleWithHeavyFlavourTrackMaker( const std::string&, ISvcLocator* ); + + StatusCode initialize() override; + + Output operator()( LHCb::Particle::Range const&, RelationTable const& ) const override; + +private: + // properties + Gaudi::Property<std::string> m_pid{this, "ParticleID", "B+"}; + + // particle property information + LHCb::ParticleProperty const* m_partprop = nullptr; + LHCb::ParticleProperty const* m_antipartprop = nullptr; + + // services + ServiceHandle<LHCb::IParticlePropertySvc> m_particlePropertySvc{this, "ParticleProperty", + "LHCb::ParticlePropertySvc"}; +}; + +DECLARE_COMPONENT( ParticleWithHeavyFlavourTrackMaker ) + +ParticleWithHeavyFlavourTrackMaker::ParticleWithHeavyFlavourTrackMaker( const std::string& name, + ISvcLocator* pSvcLocator ) + : MultiTransformer( name, pSvcLocator, {KeyValue{"InputComposites", ""}, KeyValue{"Composite2TrackRelations", ""}}, + {KeyValue{"OutputParticles", ""}, KeyValue{"OutputProtoParticles", ""}} ) {} + +StatusCode ParticleWithHeavyFlavourTrackMaker::initialize() { + return MultiTransformer::initialize().andThen( [&] { + // Get particle properties from ParticlePropertySvc + m_partprop = m_particlePropertySvc->find( m_pid ); + if ( !m_partprop ) { + throw GaudiException( "Could not find ParticleProperty for " + m_pid.value(), + "ParticleWithHeavyFlavourTrackMaker::initialize", StatusCode::FAILURE ); + } + m_antipartprop = m_partprop->antiParticle(); + return StatusCode::SUCCESS; + } ); +} + +Output ParticleWithHeavyFlavourTrackMaker::operator()( LHCb::Particle::Range const& in_particles, + RelationTable const& relation_table ) const { + auto result = Output{}; + auto& [out_particles, out_protos] = result; + + for ( auto const* in_particle : in_particles ) { + // heavy flavour track + auto track_range = relation_table.relations( in_particle ); + Track const* track = track_range.empty() ? NULL : track_range.front().to(); + + // create proto from heavy flavour track + auto* out_proto = new LHCb::ProtoParticle(); + if ( track ) out_proto->setTrack( track ); + out_protos.add( out_proto ); + + // create new particle containing proto/track and daughters (in_particle) + auto* out_particle = new LHCb::Particle( + ( out_proto->charge() == m_partprop->charge() ? m_partprop : m_antipartprop )->particleID() ); + out_particle->setProto( out_proto ); + + // add information from daughter composite + out_particle->addToDaughters( in_particle ); + + out_particle->setEndVertex( in_particle->endVertex() ); + out_particle->setReferencePoint( in_particle->referencePoint() ); + + out_particle->setMomentum( in_particle->momentum() ); + out_particle->setMeasuredMass( in_particle->measuredMass() ); + + out_particle->setPosCovMatrix( in_particle->posCovMatrix() ); + out_particle->setMomCovMatrix( in_particle->momCovMatrix() ); + out_particle->setPosMomCovMatrix( in_particle->posMomCovMatrix() ); + + out_particles.add( out_particle ); + } + + return result; +} -- GitLab From 33161ecf2583fe1e72e555183ee9a4c4fc38f58f Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Mon, 14 Aug 2023 16:57:42 +0200 Subject: [PATCH 11/14] added lastmeasurement state; added sensor and hit count functors --- Phys/FunctorCore/include/Functors/Composite.h | 33 +++++++--- Phys/FunctorCore/include/Functors/TrackLike.h | 10 +++ Phys/FunctorCore/python/Functors/__init__.py | 66 +++++++++++++++---- .../tests/src/InstantiateFunctors.cpp | 4 +- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 3 +- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 31 ++++++--- 6 files changed, 115 insertions(+), 32 deletions(-) diff --git a/Phys/FunctorCore/include/Functors/Composite.h b/Phys/FunctorCore/include/Functors/Composite.h index 9a00ace18f7..8bd68583260 100755 --- a/Phys/FunctorCore/include/Functors/Composite.h +++ b/Phys/FunctorCore/include/Functors/Composite.h @@ -322,35 +322,48 @@ namespace Functors::Composite { }; namespace BTracking { + + template <typename Composite2TrackRelations, typename Composite> + auto getTrack( Composite2TrackRelations const& relations, Composite const& composite ) { + auto const& comp = Sel::Utils::deref_if_ptr( composite ); + auto track_range = relations.relations( &comp ); + return track_range.empty() ? std::nullopt : Functors::Optional{track_range.front().to()}; + } + + /** @brief Get heavy flavour track associated to composite. */ + struct Track : public Function { + static constexpr auto name() { return "BTracking::Track"; } + + template <typename Composite2TrackRelations, typename Composite> + auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const { + return getTrack( relations, composite ); + } + }; + /** @brief Calculate the corrected mass using heavy flavour tracking info. */ struct CorrectedMass : public Function { static constexpr auto name() { return "BTracking::CorrectedMass"; } template <typename Composite2TrackRelations, typename Composite> auto operator()( Composite2TrackRelations const& relations, Composite const& composite ) const { - using FType = SIMDWrapper::scalar::float_v; using Sel::Utils::mass2; using Sel::Utils::threeMomentum; - using std::sqrt; - auto const& comp = Sel::Utils::deref_if_ptr( composite ); - // first get the right heavy flavour track - auto track_range = relations.relations( &comp ); - if ( track_range.empty() ) return FType{0}; - auto track = track_range.front().to(); + auto track = getTrack( relations, composite ); + if ( !track.has_value() ) return FType{0}; // direction using first hit (or end vertex if not available, as saved in heavy flavour track) - auto const state = track->firstState(); + auto const state = track.value()->firstState(); auto const d = LHCb::LinAlg::Vec<FType, 3>{state.tx(), state.ty(), 1.}; // Get the pT variable that we need - auto const mom = threeMomentum( comp ); + auto const mom = threeMomentum( composite ); auto const perp = mom - d * ( dot( mom, d ) / d.mag2() ); auto const pt = perp.mag(); // Calculate the corrected mass - return pt + sqrt( mass2( comp ) + pt * pt ); + return pt + sqrt( mass2( composite ) + pt * pt ); } }; } // namespace BTracking diff --git a/Phys/FunctorCore/include/Functors/TrackLike.h b/Phys/FunctorCore/include/Functors/TrackLike.h index e955aa16f0b..ea3c1d59a93 100644 --- a/Phys/FunctorCore/include/Functors/TrackLike.h +++ b/Phys/FunctorCore/include/Functors/TrackLike.h @@ -773,6 +773,16 @@ namespace Functors::Track { } }; + /** @brief Number of expected Velo clusters from VELO 3D pattern recognition + */ + struct nPRVelo3DExpect : public Function { + static constexpr auto name() { return "nPRVelo3DExpect"; } + template <typename Data> + auto operator()( Data const& d ) const { + return (int)Sel::Utils::deref_if_ptr( d ).info( LHCb::Event::Enum::Track::AdditionalInfo::nPRVelo3DExpect, -1. ); + } + }; + /** * @brief MC_Reconstructed, return the reconstructed category * for MC Particle. The input of this functor must be reconstructed diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py index 63e58f5cebc..ee5e0ec89a0 100644 --- a/Phys/FunctorCore/python/Functors/__init__.py +++ b/Phys/FunctorCore/python/Functors/__init__.py @@ -1031,6 +1031,13 @@ NFTHITS = Functor( Functor's call operator expects a track like object. """) +NPRVELO3DEXPECT = Functor( + 'NPRVELO3DEXPECT', "Track::nPRVelo3DExpect", + """Number of expected Velo clusters from VELO 3D pattern recognition. + + Functor's call operator expects a track like object. + """) + ### TRACKHISTORY = CAST_TO_INT @ Functor( 'TRACKHISTORY', "Track::History", @@ -1372,14 +1379,6 @@ MVA technique to establish pion mass hypothesis which keep in consideration the General information on PID could be found at the following `TWiki page <https://twiki.cern.ch/twiki/bin/view/LHCbPhysics/ChargedPID>`_ """) -_BPVCORRM = Functor( - '_BPVCORRM', 'Composite::CorrectedMass', - 'Compute the corrected mass of the composite using the associated [primary] vertex.' -) -_BTRACKING_BPVCORRM = Functor( - '_BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass', - 'Compute the corrected mass of the composite using the associated heavy flavour track.' -) def BPVCORRM(Vertices: DataHandle): """ Compute the corrected mass of the composite using the :py:func:`~BPV`. @@ -1413,16 +1412,61 @@ def BPVCORRMERR(Vertices: DataHandle): return _BPVCORRMERR.bind(BPV(Vertices), FORWARDARGS) -def BTRACKING_BPVCORRM(HeavyFlavourTracks: DataHandle): +def BTRACKING_TRACK(HeavyFlavourTrackRelations: DataHandle): + """ Returns heavy flavour track object associated to composite. + + Functor's call operator expects a composite like object + + Args: + HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track + + """ + _BTRACKING_TRACK = Functor('BTRACKING_TRACK', + 'Composite::BTracking::Track', + 'Heavy flavour track associated to composite.') + + return _BTRACKING_TRACK.bind(TES(HeavyFlavourTrackRelations), FORWARDARGS) + + +def BTRACKING_NHITS(HeavyFlavourTrackRelations: DataHandle): + """ Number of hits on heavy flavour track. + + Functor's call operator expects a composite like object + + Args: + HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track + + """ + return NHITS @ BTRACKING_TRACK(HeavyFlavourTrackRelations) + + +def BTRACKING_NPRVELO3DEXPECT(HeavyFlavourTrackRelations: DataHandle): + """ Number of crossed VELO sensors according to pattern recognition window. + + Functor's call operator expects a composite like object + + Args: + HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track + + """ + return NPRVELO3DEXPECT @ BTRACKING_TRACK(HeavyFlavourTrackRelations) + + +def BTRACKING_BPVCORRM(HeavyFlavourTrackRelations: DataHandle): """ Compute the corrected mass of the composite using the associated heavy flavour track. Functor's call operator expects a composite like object Args: - Vertices: DataHandle of the primary vertices + HeavyFlavourTrackRelations: DataHandle of relations of composite to heavy flavour track """ - return _BTRACKING_BPVCORRM.bind(TES(HeavyFlavourTracks), FORWARDARGS) + _BTRACKING_BPVCORRM = Functor( + 'BTRACKING_BPVCORRM', 'Composite::BTracking::CorrectedMass', + 'Compute the corrected mass of the composite using the associated heavy flavour track.' + ) + return _BTRACKING_BPVCORRM.bind( + TES(HeavyFlavourTrackRelations), FORWARDARGS) VTX_FDCHI2 = Functor( diff --git a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp index 6c08a94cff3..451e0856b2e 100644 --- a/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp +++ b/Phys/FunctorCore/tests/src/InstantiateFunctors.cpp @@ -211,7 +211,9 @@ using cmerr5 = std::invoke_result_t<Functors::Composite::CorrectedMassError, LHC // // BTracking::CorrectedMass( heavyflavourtrackrelations, composite ) // -using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass, +using cmbt1 = std::invoke_result_t<Functors::Composite::BTracking::Track, + LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>; +using cmbt2 = std::invoke_result_t<Functors::Composite::BTracking::CorrectedMass, LHCb::Relation2D<LHCb::Particle, LHCb::Event::v1::Track>, LHCb::Particle const&>; // diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index fb878247a6c..73fc4a1cf23 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -340,7 +340,8 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, // track info std::string hfbase = base + "HeavyFlavourTracking_"; sc &= tuple->column( hfbase + "Track_nSensors", - has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::VeloNSensors, 0 ) : -1. ); + has_track ? track->info( LHCb::Event::Enum::Track::AdditionalInfo::nPRVelo3DExpect, -1. ) + : -2. ); sc &= tuple->column( hfbase + "Track_nHits", has_track ? track->nHits() : -1 ); // track mc matching info bool fh_from_mcps = false; diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index dd7dd0b21e3..3a132041878 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -211,10 +211,6 @@ namespace LHCb::Pr::Velo { relations.reserve( composites.size() ); for ( LHCb::Particle const* composite : composites ) { - // set output - using namespace LHCb::Event::Enum::Track; - auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs ); - // primary vertex using Sel::Utils::endVertexPos; auto const& pv_obj = Sel::getBestPV( *composite, pvs ); @@ -227,7 +223,11 @@ namespace LHCb::Pr::Velo { // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits) auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() ); - track->addInfo( AdditionalInfo::VeloNSensors, crossedSensors.size() ); + + // build track (output; also save if no sensors crossed, as search window info is still persisted this way) + using namespace LHCb::Event::Enum::Track; + auto track = new Track( History::PrVeloHeavyFlavour, Type::Velo, PatRecStatus::PatRecIDs ); + track->addInfo( AdditionalInfo::nPRVelo3DExpect, crossedSensors.size() ); // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window) auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); @@ -241,21 +241,34 @@ namespace LHCb::Pr::Velo { // determine mother particle direction based on closest hit // in first crossed sensor or using SV-PV if not available - auto direction = found_hits.size() ? found_hits.front().position - pv : fd; + auto direction = nhits > 0 ? found_hits.front().position - pv : fd; + auto qoverp = composite->charge() / composite->p(); auto origin_state = LHCb::State(); origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(), - ( direction.y() / direction.z() ).cast(), composite->charge() / composite->p() ); + ( direction.y() / direction.z() ).cast(), qoverp ); origin_state.setLocation( LHCb::State::Location::ClosestToBeam ); track->addToStates( origin_state ); - // in case of found hit, add appropriate state - if ( found_hits.size() ) { + // in case of found hit(s), add appropriate states + if ( nhits > 0 ) { auto fh = found_hits.front(); auto fm_state = LHCb::State( origin_state.stateVector(), fh.position.z().cast(), LHCb::State::Location::FirstMeasurement ); fm_state.setX( fh.position.x().cast() ); fm_state.setY( fh.position.y().cast() ); track->addToStates( fm_state ); + + // in case more hits, add last measurement state + if ( nhits > 1 ) { + auto lh = found_hits.back(); + auto lh_direction = lh.position - pv; + auto last_state = LHCb::State(); + last_state.setState( lh.position.x().cast(), lh.position.y().cast(), lh.position.z().cast(), + ( lh_direction.x() / lh_direction.z() ).cast(), + ( lh_direction.y() / lh_direction.z() ).cast(), qoverp ); + last_state.setLocation( LHCb::State::Location::LastMeasurement ); + track->addToStates( last_state ); + } } // save track -- GitLab From 46a62e7c8c0c217e08e0836518a345883a299a18 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Wed, 16 Aug 2023 17:23:43 +0200 Subject: [PATCH 12/14] adding efficiency counters --- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 85 ++++++++++++++++++- 1 file changed, 83 insertions(+), 2 deletions(-) diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index 73fc4a1cf23..ca36d6d97d7 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -202,7 +202,22 @@ private: DataObjectReadHandle<P2TRelations> m_relations{this, "Composite2HeavyFlavourTrackRelations", ""}; // monitoring - mutable Gaudi::Accumulators::StatCounter<float> m_something{this, ""}; + mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_global_eff{this, + "Eff: has HF track | sel'd FS & HF reco'ible"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_pure_eff{this, + "HF track: MC-matched | sel'd FS & HF reco'ible"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructible_hf_track{this, "HF track: reco'ible"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructed_hf_track{this, "HF track: reco'd | reco'ible"}; + + mutable Gaudi::Accumulators::StatCounter<int> m_hftrack_nhits{this, "HF track: #hits | sel'd FS & HF track"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_hftrack_purity{this, "HF track: hit purity | sel'd FS & HF track"}; + + mutable Gaudi::Accumulators::BinomialCounter<> m_vphit_for_mchit_eff{this, "Eff: VP hit for VP MCHit"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_hit_on_track_for_vphit_with_mchit_eff{ + this, "Eff: true hit on track | sel'd FS & HF track"}; + + mutable Gaudi::Accumulators::BinomialCounter<> m_reconstructed_eff{this, "Final state (FS): reco'd"}; + mutable Gaudi::Accumulators::BinomialCounter<> m_final_state_eff{this, "Final state (FS): sel'd | reco'd"}; }; DECLARE_COMPONENT_WITH_ID( PrVeloHeavyFlavourTrackingChecker, "PrVeloHeavyFlavourTrackingChecker" ) @@ -246,7 +261,73 @@ void PrVeloHeavyFlavourTrackingChecker::operator()( MCParticles const& mcparts, if ( trackrange.size() ) track = trackrange.front(); } auto has_track = track ? true : false; - auto has_hits = has_track ? track->nHits() > 0 : false; + int n_hits = has_track ? track->nHits() : 0; + auto has_hits = n_hits > 0; + + //// counters + + // efficiencies of final state reconstruction and selection + m_reconstructed_eff += reconstructed; + if ( reconstructed ) m_final_state_eff += has_comp; + + // check MCHits associated to B+/tau+ if they are reconstructible + // and if so look if they appear on the reconstructed track in case + // the final state is reconstructed, selected and matched and has associated HF track + bool has_reconstructible_hf_track = false; + bool has_reconstructed_hf_track = false; + for ( auto const bhit : bhits ) { + auto ids_for_mchit = channelIDForMCHit.find( bhit ); + auto has_vphit_for_mchit = ids_for_mchit != channelIDForMCHit.end(); + m_vphit_for_mchit_eff += has_vphit_for_mchit; + if ( !has_vphit_for_mchit ) continue; + has_reconstructible_hf_track = true; + auto ids = ids_for_mchit->second; + int n_matches = 0; + if ( has_hits ) { + auto track_lhcbids = track->lhcbIDs(); + n_matches = std::count_if( ids.begin(), ids.end(), [&track_lhcbids]( auto whit ) { + return std::count_if( track_lhcbids.begin(), track_lhcbids.end(), + [&whit]( auto lhcbid ) { return lhcbid.vpID() == whit.id; } ) > 0; + } ); + } + if ( reconstructed && has_comp && has_hits ) { + has_reconstructed_hf_track |= n_matches > 0; + m_hit_on_track_for_vphit_with_mchit_eff += n_matches > 0; + } + } + + float track_purity = -1.f; + if ( reconstructed && has_comp && has_hits ) { + int n_matched_hits = 0; + for ( auto const lhcbid : track->lhcbIDs() ) { + unsigned int channelid = lhcbid.vpID(); + int n_match_per_hit = + bhits.empty() ? 0 + : std::count_if( bhits.begin(), bhits.end(), [&channelIDForMCHit, &channelid]( auto mchit ) { + auto map_iter = channelIDForMCHit.find( mchit ); + int n_match_per_hit_per_mchit = + ( map_iter == channelIDForMCHit.end() ) + ? 0 + : std::count_if( map_iter->second.begin(), map_iter->second.end(), + [&channelid]( auto whit ) { return whit.id == channelid; } ); + return n_match_per_hit_per_mchit > 0; + } ); + if ( n_match_per_hit > 0 ) n_matched_hits++; + m_hftrack_purity += n_match_per_hit > 0; + } + track_purity = n_matched_hits / (float)n_hits; + m_hftrack_nhits += n_hits; + } + + m_reconstructible_hf_track += has_reconstructible_hf_track; + if ( has_reconstructible_hf_track ) { + m_reconstructed_hf_track += has_reconstructed_hf_track; + if ( reconstructed && has_comp ) { + m_hftrack_global_eff += has_hits; + m_hftrack_pure_eff += has_hits && track_purity >= 0.7; + } + } + //// fill ntuple with info std::string base = "MCP_"; // mc truth -- GitLab From f815da3cec8a0f55a55d10ba35af7720a06227be Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Thu, 24 Aug 2023 17:08:22 +0200 Subject: [PATCH 13/14] added control of charge setting --- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index 3a132041878..6f861b214ba 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -174,6 +174,7 @@ namespace LHCb::Pr::Velo { // properties Gaudi::Property<float> m_max_distance{this, "MaxDistanceFromPVSV", 0.5 * Gaudi::Units::mm, "maximum hit distance from PV-SV line (search window)"}; + Gaudi::Property<bool> m_revert_charge{this, "RevertCharge", false, "Revert charge with respect to composite"}; // statistics monitoring mutable Gaudi::Accumulators::StatCounter<> m_nInSensor{this, "Nb of sensors crossed"}; @@ -242,7 +243,8 @@ namespace LHCb::Pr::Velo { // determine mother particle direction based on closest hit // in first crossed sensor or using SV-PV if not available auto direction = nhits > 0 ? found_hits.front().position - pv : fd; - auto qoverp = composite->charge() / composite->p(); + auto charge = composite->charge() != 0 ? ( m_revert_charge ? -1 : 1 ) * composite->charge() : 1; + auto qoverp = charge / composite->p(); auto origin_state = LHCb::State(); origin_state.setState( pv.x().cast(), pv.y().cast(), pv.z().cast(), ( direction.x() / direction.z() ).cast(), ( direction.y() / direction.z() ).cast(), qoverp ); -- GitLab From 4170bbfc06736b09865d6cfd9bf27a8cc765e331 Mon Sep 17 00:00:00 2001 From: Maarten Van Veghel <maarten.vanveghel@cern.ch> Date: Wed, 1 Nov 2023 18:47:29 +0100 Subject: [PATCH 14/14] sensor search with distance to sensor (full cylinder assumption) --- .../tests/options/test_functors_tested.py | 9 ++-- .../src/PrVeloHeavyFlavourTrackingChecker.cpp | 2 +- Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp | 47 +++++++------------ 3 files changed, 23 insertions(+), 35 deletions(-) diff --git a/Phys/FunctorCore/tests/options/test_functors_tested.py b/Phys/FunctorCore/tests/options/test_functors_tested.py index 5eedfb9d66e..8348dc2dbd5 100644 --- a/Phys/FunctorCore/tests/options/test_functors_tested.py +++ b/Phys/FunctorCore/tests/options/test_functors_tested.py @@ -65,10 +65,11 @@ exceptions = [ 'Functors::Track::Flag', 'Functors::Track::TX', 'Functors::Track::TY', 'Functors::Track::history', 'Functors::Track::nDoF', 'Functors::Track::nFTHits', 'Functors::Track::nUTHits', - 'Functors::Track::nVPHits', 'Functors::Common::ForwardArgs<>( ', - 'Functors::Common::ForwardArgs<0>( ', 'Functors::Common::ForwardArgs<1>( ', - 'Functors::Common::ForwardArgs<2>( ', 'Functors::Common::UnitVector', - 'Functors::Common::TES', 'Functors::Common::Call', 'Functors::Common::Dot', + 'Functors::Track::nVPHits', 'Functors::Track::nPRVelo3DExpect', + 'Functors::Common::ForwardArgs<>( ', 'Functors::Common::ForwardArgs<0>( ', + 'Functors::Common::ForwardArgs<1>( ', 'Functors::Common::ForwardArgs<2>( ', + 'Functors::Common::UnitVector', 'Functors::Common::TES', + 'Functors::Common::Call', 'Functors::Common::Dot', 'Functors::Common::NormedDot', 'Functors::Common::BestPV', 'Functors::Common::ToLinAlg', 'Functors::Common::Phi_Coordinate', 'Functors::Examples::PlusN', diff --git a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp index ca36d6d97d7..f270550f811 100644 --- a/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp +++ b/Pr/PrMCTools/src/PrVeloHeavyFlavourTrackingChecker.cpp @@ -116,7 +116,7 @@ namespace { size_t nmatches = 0; for ( LHCb::MCParticle const* pion : mcps.pions ) { map.applyToAllLinks( - [&pion, &mcps, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) { + [&pion, &mcparts, &min_weight, &nmatches]( unsigned int, unsigned int mcPartKey, float weight ) { if ( ( pion == static_cast<LHCb::MCParticle const*>( mcparts.containedObject( mcPartKey ) ) ) && ( weight > min_weight ) ) nmatches += 1; diff --git a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp index 6f861b214ba..986db9a3fc2 100644 --- a/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp +++ b/Pr/PrPixel/src/VeloHeavyFlavourTracking.cpp @@ -48,12 +48,6 @@ namespace LHCb::Pr::Velo { friend bool operator<( SensorInfo const& lhs, SensorInfo const& rhs ) { return std::pair{lhs.z, lhs.id} < std::pair{rhs.z, rhs.id}; } - // module information - auto module() const { - auto cid = ChannelID(); - cid.setSensor( this->id ); - return cid.module(); - } }; using SensorInfos = std::vector<SensorInfo>; @@ -95,25 +89,24 @@ namespace LHCb::Pr::Velo { // functions template <typename Position3, typename Vector3> SensorInfos getCrossedSensors( DeVP const& velo, SensorInfos const& sensors, Position3 const& start, - Position3 const& end, Vector3 const& slopes, float margin ) { + Position3 const& end, Vector3 const& slopes, float max_distance_squared ) { SensorInfos crossedSensors; - auto pos = start; - auto margin_normalized = margin / sqrt( slopes.mag2() - 1 ); + auto pos = start; for ( auto sensor : sensors ) { if ( sensor.z < pos.Z() ) continue; if ( sensor.z > end.Z() ) break; - // if in range look if search window intersects with sensor (up to some margin deeper into VELO) - pos = pos + ( sensor.z - pos.Z() ) * slopes; - auto pos_with_margin = pos + margin_normalized * slopes; - if ( velo.sensor( sensor.id ).isInsideSensor( pos.X().cast(), pos.Y().cast() ) | - velo.sensor( sensor.id ).isInsideSensor( pos_with_margin.X().cast(), pos_with_margin.Y().cast() ) ) - crossedSensors.push_back( {pos.Z().cast(), sensor.id} ); + // if in range look if search window intersects with sensor + pos = pos + ( sensor.z - pos.Z() ) * slopes; + auto dist = + velo.sensor( sensor.id ) + .closestDistanceToSensor( ROOT::Math::XYZPoint( pos.X().cast(), pos.Y().cast(), pos.Z().cast() ) ); + if ( dist.perp2() < max_distance_squared ) crossedSensors.push_back( {pos.Z().cast(), sensor.id} ); } return crossedSensors; } template <typename Position3, typename Vector3> - HitInfos getHitsForCrossedModules( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start, + HitInfos getHitsForCrossedSensors( SensorInfos const& sensors, VeloHits const& hits, Position3 const& start, Vector3 const& slopes, FType max_distance ) { HitInfos poshits; if ( !sensors.size() ) return poshits; @@ -122,25 +115,19 @@ namespace LHCb::Pr::Velo { auto pos = start; // check sensors for ( auto const sensor : sensors ) { - // its associated module (to not be affected by borders between sensors) - auto crossed_module = sensor.module(); // go to sensor location pos = pos + ( sensor.z - pos.z().cast() ) * slopes; - // look if this sensor (and associated module) has hits and if they are in search window + // look if this sensor has hits and if they are in the search window for ( auto hit : hits.scalar() ) { auto id = ChannelID( hit.get<LHCb::Pr::VP::VPHitsTag::ChannelId>().cast() ); - if ( id.module() == crossed_module ) { + if ( id.sensor() == sensor.id ) { auto hitpos = hit.get<LHCb::Pr::VP::VPHitsTag::pos>(); - // make sure they are in the same plane up to some tolerance - // (not all sensors in a module are exactly) - if ( abs( hitpos.z().cast() - sensor.z ) < 0.1 * Gaudi::Units::mm ) { - auto d = ( hitpos - pos ).mag(); - if ( d < max_distance ) poshits.push_back( {d.cast(), hitpos, id} ); - } + auto dist = ( hitpos - pos ).rho(); + if ( dist < max_distance ) poshits.push_back( {dist.cast(), hitpos, id} ); } } } - // sort and remove duplicates (can arise with only module/plane (not sensor) checking) + // sort and remove duplicates (shouldn't be possible in this logic, but leave out of paranoia) std::sort( poshits.begin(), poshits.end() ); poshits.erase( std::unique( poshits.begin(), poshits.end() ), poshits.end() ); return poshits; @@ -222,8 +209,8 @@ namespace LHCb::Pr::Velo { auto const fd = end - pv; auto const slopes = fd / fd.Z(); - // crossed sensors in PV-SV line (or cylinder, with non-zero margin; should be same as for hits) - auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, m_max_distance.value() ); + // crossed sensors in PV-SV cylinder (with m_max_distance as radius) + auto crossedSensors = getCrossedSensors( velo, sensors(), pv, end, slopes, std::pow( m_max_distance, 2 ) ); // build track (output; also save if no sensors crossed, as search window info is still persisted this way) using namespace LHCb::Event::Enum::Track; @@ -231,7 +218,7 @@ namespace LHCb::Pr::Velo { track->addInfo( AdditionalInfo::nPRVelo3DExpect, crossedSensors.size() ); // look for hits in relevant sensors (duplicates removed and sorted in z and distance to search window) - auto found_hits = getHitsForCrossedModules( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); + auto found_hits = getHitsForCrossedSensors( crossedSensors, velohits, pv, slopes, m_max_distance.value() ); auto nhits = found_hits.size(); // store sorted hits -- GitLab