Helper class to group muon sgements that may belong to a muon trajectory. More...

#include <MsTrackSeeder.h>

Inheritance diagram for MuonR4::MsTrackSeeder:

Collaboration diagram for MuonR4::MsTrackSeeder:

Classes
struct	Config
	Configuration object. More...

Public Types
enum class	SeedCoords : std::uint8_t { eDetSection , eSector , ePosOnCylinder }
	Abrivation of the seed coordinates. More...
using	SearchTree_t = Acts::KDTree<3, const xAOD::MuonSegment*, double, std::array, 6>
	Definition of the search tree class.
using	Location = MsTrackSeed::Location
	Enum toggling whether the segment is in the endcap or barrel.
using	SectorProjector = ExpandedSector::SectorProjector
	Recycle the expanded sector.
using	VecOpt_t = std::optional<Amg::Vector3D>
using	PosMomPair_t = std::pair<Amg::Vector3D, Amg::Vector3D>

Public Member Functions
	MsTrackSeeder (const std::string &msgName, Config &&cfg)
	Standard constructor.
SearchTree_t	constructTree (const xAOD::MuonSegmentContainer &segments) const
	Construct a complete search tree from a MuonSegment container.
Amg::Vector2D	expressOnCylinder (const xAOD::MuonSegment &segment, const Location loc, const ExpandedSector sector) const
	Expresses the segment on the cylinder surface.
double	estimateQtimesP (const AtlasFieldCacheCondObj &magField, const Amg::Vector3D &planeNorm, const PosMomPair_t &p1, const PosMomPair_t &p2, const PosMomPair_t &p3) const
	Estimate the charge times momentum of a muon candidate when three points are available.
double	estimateQtimesP (const AtlasFieldCacheCondObj &magField, const Amg::Vector3D &planeNorm, const PosMomPair_t &p1, const PosMomPair_t &p2) const
	Estimate the charge times momentum of a muon candidate when two points are available.
double	estimateQtimesP (const AtlasFieldCacheCondObj &magField, const MsTrackSeed &seed) const
	Estimate the charge times momentum of a muon track candidate from the contained segments.
bool	withinBounds (const Amg::Vector2D &projPos, const Location loc) const
	Returns whether the expression on the cylinder is within the surface bounds.
std::unique_ptr< MsTrackSeedContainer >	findTrackSeeds (const EventContext &ctx, const xAOD::MuonSegmentContainer &segments) const
	Constructs the MS track seeds from the segment container.
bool	msgLvl (const MSG::Level lvl) const
	Test the output level.
MsgStream &	msg () const
	The standard message stream.
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream.
void	setLevel (MSG::Level lvl)
	Change the current logging level.

Static Public Member Functions
static Amg::Vector3D	segPosOntoPhiPlane (const Amg::Vector3D &planeNorm, const int Sector, const Amg::Vector3D &posToProject)
	Projects the segment position onto the plane with global phi = x The local coordinate system is arranged such that the x-axis is co-linear to the phi direction.
static Amg::Vector3D	segDirOntoPhiPlane (const Amg::Vector3D &planeNorm, const Amg::Vector3D &dirToProject)
	Projects the segment direction onto the plane with global phi = x by removing the component orthogonal to the plane.

Private Types
using	TreeRawVec_t = SearchTree_t::vector_t
	Abbrivation of the KDTree raw data vector.

Private Member Functions
double	getPtimesQ (const Amg::Vector3D &forceIntegral, const Amg::Vector3D &deltaDir) const
	Compute the charge times momentum from the integral of lorentz force and the total change in direction.
Amg::Vector3D	forceIntegration (const PosMomPair_t &point1, const PosMomPair_t &point2, const Amg::Vector3D &planeNorm, MagField::AtlasFieldCache &fieldCache) const
	Compute the integral of magnetic force (v x B ) dS along a trajectory, given the initial and final positions and directions of the trajectory.
void	appendSegment (const xAOD::MuonSegment *segment, const Location loc, TreeRawVec_t &outContainer) const
	Append the to the raw data container.
std::unique_ptr< MsTrackSeedContainer >	resolveOverlaps (MsTrackSeedContainer &&unresolved) const
	Removes exact duplciates or partial subsets of the MsTrackSeeds.
void	initMessaging () const
	Initialize our message level and MessageSvc.

Private Attributes
std::set< double >	m_fieldExtpSteps {}
Config	m_cfg {}
std::string	m_nm
	Message source name.
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels).
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer.
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level.
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging).

Detailed Description

Helper class to group muon sgements that may belong to a muon trajectory.

The reconstructed muon segments are projected onto the surface of a cylinder crossing roughly the middle stations of the MS. They are then appended to a 3-dimensional search tree using the extrapolated coordinate on the cylnder, the cylinder surface index and the segment's associated sector number.

Definition at line 29 of file MsTrackSeeder.h.

Member Typedef Documentation

◆ Location

using MuonR4::MsTrackSeeder::Location = MsTrackSeed::Location

Enum toggling whether the segment is in the endcap or barrel.

Definition at line 56 of file MsTrackSeeder.h.

◆ PosMomPair_t

using MuonR4::MsTrackSeeder::PosMomPair_t = std::pair<Amg::Vector3D, Amg::Vector3D>

Definition at line 102 of file MsTrackSeeder.h.

◆ SearchTree_t

using MuonR4::MsTrackSeeder::SearchTree_t = Acts::KDTree<3, const xAOD::MuonSegment*, double, std::array, 6>

Definition of the search tree class.

Definition at line 54 of file MsTrackSeeder.h.

◆ SectorProjector

using MuonR4::MsTrackSeeder::SectorProjector = ExpandedSector::SectorProjector

Recycle the expanded sector.

Definition at line 58 of file MsTrackSeeder.h.

◆ TreeRawVec_t

using MuonR4::MsTrackSeeder::TreeRawVec_t = SearchTree_t::vector_t

private

Abbrivation of the KDTree raw data vector.

Definition at line 171 of file MsTrackSeeder.h.

◆ VecOpt_t

using MuonR4::MsTrackSeeder::VecOpt_t = std::optional<Amg::Vector3D>

Definition at line 59 of file MsTrackSeeder.h.

Member Enumeration Documentation

◆ SeedCoords

enum class MuonR4::MsTrackSeeder::SeedCoords : std::uint8_t

strong

Abrivation of the seed coordinates.

Enumerator
eDetSection	Encode the seed location (-1,1 -> endcaps, 0 -> barrel.
eSector	Sector of the associated spectrometer sector.
ePosOnCylinder	Extrapolation position along the cylinder surface.

Definition at line 61 of file MsTrackSeeder.h.

                                  : std::uint8_t{
                eDetSection,
                eSector,
                ePosOnCylinder
            };

Constructor & Destructor Documentation

◆ MsTrackSeeder()

MuonR4::MsTrackSeeder::MsTrackSeeder	(	const std::string &	msgName,
		Config &&	cfg )

Standard constructor.

Parameters

msgName	Name of the seeder's msgStream
cfg	Configured cylinder dimensions, cuts & selection tool

Initialize the field extraction steps

Definition at line 35 of file MsTrackSeeder.cxx.

                                                                      :
        AthMessaging{msgName},
        m_cfg{std::move(cfg)}{
       
        const double stepSize {1. / static_cast<double>(m_cfg.nFieldSteps)};
        for (std::size_t i = 0; i < m_cfg.nFieldSteps; ++i) {
            m_fieldExtpSteps.insert((static_cast<double>(i) + 0.5) * stepSize);
        }
    }

Member Function Documentation

◆ appendSegment()

void MuonR4::MsTrackSeeder::appendSegment	(	const xAOD::MuonSegment *	segment,
		const Location	loc,
		TreeRawVec_t &	outContainer ) const

private

Append the to the raw data container.

If the projection onto the barrel cylinder / endcap discs exceeds the bounds, the segment is not added. Segments in sector 1/16 are mirrored into sector 0/17 to complete the search range

Parameters

segment	Pointer to the segment to add
loc	Switch whether the segment shall be projected onto barrel/endcap
outContainer	Raw KDTree data vector where the segment is appended

Check whether the segment belongs to the left or right sector as well

Blow-up the number of sectors by a factor of 2. The even numbers represent the segments expressed @ the sector centre. The odd numbers represent the overlap region between two adjacent sectors. For sector 16, the right overlap region is mapped to 1

Enumeration to indicate whether the segment is expressed on the negative endcap (-1), the barrel (0) or the positive endcap

Coordinate on the cylinder

Definition at line 300 of file MsTrackSeeder.cxx.

                                                                        {
        
        const unsigned segSector = segment->sector();
        for (const auto proj : {SectorProjector::leftOverlap, SectorProjector::center, SectorProjector::rightOverlap}) {
            const ExpandedSector projSector{segSector, proj};
            if (segment->nPhiLayers() > 0 &&  projSector != ExpandedSector{segment->position().phi()}) {
                ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Segment @"<<Amg::toString(segment->position())
                    <<" is not in sector "<<projSector);
                continue;
            }
            const Amg::Vector2D refPoint{expressOnCylinder(*segment, loc, projSector)};
            if (!withinBounds(refPoint, loc)) {
                 continue;
            }
            using enum SeedCoords;
            std::array<double, 3> coords{};
            coords[Acts::toUnderlying(eSector)] = projSector.sector();
            coords[Acts::toUnderlying(eDetSection)] =  Acts::copySign(Acts::toUnderlying(loc), refPoint[1]);
            coords[Acts::toUnderlying(ePosOnCylinder)] = refPoint[Location::Barrel == loc];
            ATH_MSG_VERBOSE("Add segment "<<print(*segment)<<", seed quality: "
            <<m_cfg.selector->passSeedingQuality(Gaudi::Hive::currentContext(), *detailedSegment(*segment))
            <<" with "<<coords<<" to the search tree");
            outContainer.emplace_back(std::move(coords), segment);
        }
    }

◆ constructTree()

SearchTree_t MuonR4::MsTrackSeeder::constructTree ( const xAOD::MuonSegmentContainer & segments ) const

Construct a complete search tree from a MuonSegment container.

Parameters

segments Reference to the segment container to construct.

Definition at line 334 of file MsTrackSeeder.cxx.

                                                                                           {
        TreeRawVec_t rawData{};
        rawData.reserve(3*segments.size());
        for (const xAOD::MuonSegment* segment : segments){
            appendSegment(segment, Location::Barrel, rawData);
            appendSegment(segment, Location::Endcap, rawData);
        }
        ATH_MSG_VERBOSE("Create a new tree with "<<rawData.size()<<" entries. ");
        return SearchTree_t{std::move(rawData)};
    }

◆ estimateQtimesP() [1/3]

double MuonR4::MsTrackSeeder::estimateQtimesP	(	const AtlasFieldCacheCondObj &	magField,
		const Amg::Vector3D &	planeNorm,
		const PosMomPair_t &	p1,
		const PosMomPair_t &	p2 ) const

Estimate the charge times momentum of a muon candidate when two points are available.

The position and direction of each point need to be projected onto the given phi plane, and the muon trajectory is approximated as 2D trajectory within this plane.

Parameters

magField	Magnetic field
planeNorm	Normal of the bending plane containing the muon trajectory in this simplified approach
seg1	First segment
seg2	Second segment

Returns: : Estimated Q*P value

Definition at line 175 of file MsTrackSeeder.cxx.

                                                                        {
        MagField::AtlasFieldCache fieldCache{};
        magField.getInitializedCache(fieldCache);
 
        return getPtimesQ(forceIntegration(p1, p2, planeNorm, fieldCache),
                          p2.second - p1.second);
    }

◆ estimateQtimesP() [2/3]

double MuonR4::MsTrackSeeder::estimateQtimesP	(	const AtlasFieldCacheCondObj &	magField,
		const Amg::Vector3D &	planeNorm,
		const PosMomPair_t &	p1,
		const PosMomPair_t &	p2,
		const PosMomPair_t &	p3 ) const

Estimate the charge times momentum of a muon candidate when three points are available.

The given position and direction of each point need to be projected onto the given phi plane, and the muon trajectory is approximated as 2D trajectory within this plane.

Parameters

magField	Magnetic field
planeNorm	Normal of the bending plane containing the muon trajectory in this simplified approach
seg1	First segment
seg2	Second segment
seg3	Third segment

Returns: : Estimated Q*P value

Definition at line 104 of file MsTrackSeeder.cxx.

                                                                        {
        // When 3 points are available, we can use each pair of segments to estimate 
        // the momentum and charge, and then combine the estimates. To make the 
        // combination, we define a struct to hold each PtimesQ estimate
        struct Estimate {
            double PtimesQ{0.};
            // Weighting factor based on the magnitude of the integrated force.
            double weight{0.};
            // Scores as a combination of charge agreement and momentum deviation.
            double score{0.};
        };
        
        MagField::AtlasFieldCache fieldCache{};
        magField.getInitializedCache(fieldCache);
 
        const Amg::Vector3D force12 {forceIntegration(p1, p2, planeNorm, fieldCache)};
        const Amg::Vector3D force23 {forceIntegration(p2, p3, planeNorm, fieldCache)};
        const Amg::Vector3D force13 {force12 + force23};
        
        std::vector<Estimate> estimates{};
        const double weightNorm {force12.mag() + force23.mag()};
        // Pairwise momentum estimates: 12 and 23
        estimates.emplace_back(getPtimesQ(force12, p2.second - p1.second), force12.mag()/weightNorm, 0.);
        estimates.emplace_back(getPtimesQ(force23, p3.second - p2.second), force23.mag()/weightNorm, 0.);
        // Two estimates for pair 13: using segment directions and using position differences
        const double weight13 {force13.mag()/weightNorm};
        estimates.emplace_back(getPtimesQ(force13, p3.second - p1.second), weight13, 0.);
        const Amg::Vector3D t12 {(p2.first - p1.first).unit()};
        const Amg::Vector3D t23 {(p3.first - p2.first).unit()};
        estimates.emplace_back(getPtimesQ(force13, t23 - t12), weight13, 0.);
 
        for (std::size_t i {0}; i < estimates.size(); ++i) {
            Estimate& est1 {estimates[i]};
            for (std::size_t j {i+1}; j < estimates.size(); ++j) {
                Estimate& est2 {estimates[j]};
                // Compute the charge agreement score & momentum deviation
                double w {std::min(est1.weight, est2.weight)};
                double chargeScore {chargeAgree(est1.PtimesQ, est2.PtimesQ) ? 1. : -1.};
                double pDevPenalty {momentumDev(est1.PtimesQ, est2.PtimesQ)};
                // Update the scores
                est1.score += w * (chargeScore - pDevPenalty);
                est2.score += w * (chargeScore - pDevPenalty);
            }
        }
        if (msgLvl(MSG::VERBOSE)) {
            std::vector<std::string> names {"Pair01", "Pair12", "Pair02Seg", "Pair02Pos"};
            for (const auto& [i, est] : Acts::enumerate(estimates)) {
                ATH_MSG_VERBOSE(__func__<<"() Estimate "<<names[i]<<": PtimesQ: "<<est.PtimesQ*1e-3
                    <<", weight: "<<est.weight<<", score: "<<est.score);
            }
        }
        // Find the best charge estimate and estimate the final momentum as a weighted average of the 
        // estimates that agree with the best charge
        const Estimate& bestEstimate {*std::ranges::max_element(estimates, 
            std::ranges::less{}, &Estimate::score)};
        const double charge {std::copysign(1., bestEstimate.PtimesQ)};
 
        double totalSum {0.}, totalWeight {0.};
        for (const Estimate& est : estimates) {
            if (chargeAgree(est.PtimesQ, charge)) {
                totalSum += est.PtimesQ * est.weight;
                totalWeight += est.weight;
            }
        }
        assert(totalWeight > Acts::s_epsilon);
        return totalSum / totalWeight;
    }

◆ estimateQtimesP() [3/3]

double MuonR4::MsTrackSeeder::estimateQtimesP	(	const AtlasFieldCacheCondObj &	magField,
		const MsTrackSeed &	seed ) const

Estimate the charge times momentum of a muon track candidate from the contained segments.

The position and direction of the segments are projected onto a given phi plane, defined by the segments with phi information or the sector plane.
The muon trajectory is approximated as 2D trajectory within this plane to avoid side effects from (non)-present phi measurements.

Parameters

magFiel	Reference to the magnetic field holder
seed	Reference to the seed of interest.

Calculate the averaged phi from the segments

If less than 3 segments are found check whether the track crosses the BEE or EE chamber and use that segment as the third one.

Definition at line 219 of file MsTrackSeeder.cxx.

                                                                         {
        using namespace Muon::MuonStationIndex;
        MagField::AtlasFieldCache fieldCache{};
        magField.getInitializedCache(fieldCache); 
        double deltaPhiAcc {0.};
        std::optional<double> centralPhi {};
        unsigned nSegsWithPhi{0};
        for (const xAOD::MuonSegment* segment : seed.segments()) {
            if (segment->nPhiLayers() > 0) {
                if (!centralPhi) centralPhi = segment->position().phi();
                deltaPhiAcc += P4Helpers::deltaPhi(*centralPhi, segment->position().phi());
                ++nSegsWithPhi;
            }
        }
        const double circPhi {nSegsWithPhi > 0 
            ? P4Helpers::deltaPhi(*centralPhi + deltaPhiAcc / nSegsWithPhi, 0.) 
            : seed.sector().phi()};
 
        std::array<const xAOD::MuonSegment*, 3> segmentsToUse{};
        // Try first to find segments in the inner, middle and outer layers. If both a barrel
        // and endcap segments are present in the same layer, the barrel segment is preferred.
        for (const xAOD::MuonSegment* segment : seed.segments()) {
            ChIndex chIndex {segment->chamberIndex()};
            switch(toLayerIndex(chIndex)) {
                case LayerIndex::Inner:
                    if (!segmentsToUse[0] || isBarrel(chIndex)) {
                        segmentsToUse[0] = segment;
                    }
                    break;
                case LayerIndex::Middle:
                    if (!segmentsToUse[1] || isBarrel(chIndex)) {
                        segmentsToUse[1] = segment;
                    }
                    break;
                case LayerIndex::Outer:
                    if (!segmentsToUse[2] || isBarrel(chIndex)) {
                        segmentsToUse[2] = segment;
                    }
                    break;
                default:
                    break;
            }
        }
        unsigned nSegments = std::ranges::count_if(segmentsToUse, 
            [](const xAOD::MuonSegment* seg) { return seg != nullptr; });
        if (nSegments < 3) {
            auto missingSeg = std::ranges::find(segmentsToUse, nullptr);
            for (const xAOD::MuonSegment* segment : seed.segments()) {
                LayerIndex layIndex {toLayerIndex(segment->chamberIndex())};
                if (layIndex != LayerIndex::Extended && layIndex != LayerIndex::BarrelExtended) {
                    continue;
                }
                assert(missingSeg != segmentsToUse.end());
                *missingSeg = segment;
                nSegments++;
                if (nSegments == 3) {
                    break;
                }
                missingSeg = std::ranges::find(segmentsToUse, nullptr);
            }
            std::ranges::sort(segmentsToUse, [](const xAOD::MuonSegment* seg1, const xAOD::MuonSegment* seg2) {
                if (!seg1 || !seg2) {
                    return seg1 != nullptr;
                }
                return seg1->position().perp() < seg2->position().perp();
            });
        }
        const Amg::Vector3D planeNorm {Acts::makeDirectionFromPhiTheta(circPhi + 90._degree, 90._degree)};
        auto point = [&planeNorm](const xAOD::MuonSegment* seg) {
            return std::make_pair(segPosOntoPhiPlane(planeNorm, seg->sector(), seg->position()),
                                  segDirOntoPhiPlane(planeNorm, seg->direction()));
        };
        return nSegments == 3 
            ? estimateQtimesP(magField, planeNorm, point(segmentsToUse[0]), point(segmentsToUse[1]), point(segmentsToUse[2]))
            : estimateQtimesP(magField, planeNorm, point(segmentsToUse[0]), point(segmentsToUse[1]));
    }

◆ expressOnCylinder()

Amg::Vector2D MuonR4::MsTrackSeeder::expressOnCylinder	(	const xAOD::MuonSegment &	segment,
		const Location	loc,
		const ExpandedSector	sector ) const

Expresses the segment on the cylinder surface.

Parameters

segment	Reference to the segment of consideration
loc	Surface location: [barrel/endcap]
expandedSector	The expanded sector number taking the overlap regions between large and small sectors into account [0;32]

extrapolated position

Definition at line 59 of file MsTrackSeeder.cxx.

                                                                                      {
        const Amg::Vector3D pos{segPosOntoPhiPlane(
            sector.normalDir(), segment.sector(), segment.position())};
        const Amg::Vector3D dir{segment.direction()};
 
        const Amg::Vector2D projPos{pos.perp(), pos.z()};
        const Amg::Vector2D projDir{dir.perp(), dir.z()};
 
        ATH_MSG_VERBOSE( "segment position:" << Amg::toString(segment.position())<< ", direction: " << Amg::toString(segment.direction()) );
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Express segment in @"<<Amg::toString(pos)
                        <<", direction: "<<Amg::toString(dir)<< " sector projector: " << sector 
                        << " location: " << Acts::toUnderlying(loc));
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Projected position onto sector: "<<Amg::toString(projPos)
                        <<", projected direction: "<<Amg::toString(projDir));
 
        double lambda{0.};
        if (Location::Barrel == loc) {
            lambda = Amg::intersect<2>(projPos, projDir, Amg::Vector2D::UnitX(), 
                                        m_cfg.barrelRadius).value_or(10. * Gaudi::Units::km);
                                        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Intersect with barrel at radius: "<<m_cfg.barrelRadius<<" --> "<<Amg::toString(projPos + lambda * projDir));
        } else {
            lambda = Amg::intersect<2>(projPos, projDir, Amg::Vector2D::UnitY(), 
                                       Acts::copySign(m_cfg.endcapDiscZ, projPos[1])).value_or(10. * Gaudi::Units::km);
                                       ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Intersect with endcap at z: "<<Acts::copySign(m_cfg.endcapDiscZ, projPos[1])<<" --> "<<Amg::toString(projPos + lambda * projDir));
        }
        return projPos + lambda * projDir;  
    }

◆ findTrackSeeds()

std::unique_ptr< MsTrackSeedContainer > MuonR4::MsTrackSeeder::findTrackSeeds	(	const EventContext &	ctx,
		const xAOD::MuonSegmentContainer &	segments ) const

Constructs the MS track seeds from the segment container.

Parameters

ctx	EventContext to access conditions / event data
segments	Refrence to the overall event's segment container

Bad segment not suitable for track seeding or the segment coordinates are just mirrored at the overlap between sector 1 -> 16

Define the search range.

Ensure that only endcap / barrel seeds are considered. The values are integers -> add tiny margin

Move 25 cm along the projected plane

Include the neighbouring sectors

Using the cube above, let the tree search for all compatible segments

Ensure that the sector overlap and momentum vectors are compatible with a MS trajectory

No segments were combined

Calculate the seed's position

Definition at line 344 of file MsTrackSeeder.cxx.

                                                                                                                        {
        SearchTree_t orderedSegs{constructTree(segments)};
        MsTrackSeedContainer trackSeeds{};
        using enum SeedCoords;
        for (const auto& [coords, seedCandidate] : orderedSegs) {
            const Segment* recoSeedCandidate = detailedSegment(*seedCandidate);
             if (!m_cfg.selector->passSeedingQuality(ctx, *recoSeedCandidate)){
                ATH_MSG_VERBOSE("Segment "<<print(*seedCandidate)<<" does not pass the seeding quality.");
                continue;
            }    
            SearchTree_t::range_t selectRange{};
            selectRange[Acts::toUnderlying(eDetSection)].shrink(coords[Acts::toUnderlying(eDetSection)] - 0.1, 
                                                                coords[Acts::toUnderlying(eDetSection)] + 0.1);
            selectRange[Acts::toUnderlying(ePosOnCylinder)].shrink(coords[Acts::toUnderlying(ePosOnCylinder)] - m_cfg.seedHalfLength, 
                                                                   coords[Acts::toUnderlying(ePosOnCylinder)] + m_cfg.seedHalfLength);
            selectRange[Acts::toUnderlying(eSector)].shrink(coords[Acts::toUnderlying(eSector)] -0.25, 
                                                            coords[Acts::toUnderlying(eSector)] +0.25);
            
            MsTrackSeed newSeed{static_cast<Location>(std::abs(coords[Acts::toUnderlying(eDetSection)])),
                                ExpandedSector{static_cast<std::int8_t>(coords[Acts::toUnderlying(eSector)])}};
            ATH_MSG_VERBOSE("Search for compatible segments to "<<print(*seedCandidate)<<".");
            orderedSegs.rangeSearchMapDiscard(selectRange, [&](
                    const SearchTree_t::coordinate_t& /*coords*/,
                    const xAOD::MuonSegment* extendWithMe) {
                        const Segment* extendCandidate = detailedSegment(*extendWithMe);
                        if (!m_cfg.selector->compatibleForTrack(ctx, *recoSeedCandidate, *extendCandidate)) {
                            ATH_MSG_VERBOSE("Segment "<<print(*extendWithMe)<<" is not compatible.");
                            return;
                        }
                        auto itr = std::ranges::find_if(newSeed.segments(), [extendWithMe](const xAOD::MuonSegment* onSeed){
                            return extendWithMe->chamberIndex() == onSeed->chamberIndex();
                        });
                        if (itr == newSeed.segments().end()){
                            ATH_MSG_VERBOSE("Add segment "<<print(*extendWithMe)<<" to seed.");
                            newSeed.addSegment(extendWithMe);
                        }
                        else if (reducedChi2(**itr) > reducedChi2(*extendWithMe) &&
                                     (*itr)->nPhiLayers() <= extendWithMe->nPhiLayers()) {
 
                             ATH_MSG_VERBOSE("Replace segment "<<print(**itr)<<" with "<<print(*extendWithMe)
                                             <<" on seed due to better chi2.");
                            newSeed.replaceSegment(*itr, extendWithMe);
                        }
            });
            if (newSeed.segments().empty()) {
                continue;
            }
 
            newSeed.addSegment(seedCandidate);
 
            // Let's check if we build a single station seed and if yes reject it.
            using namespace Muon::MuonStationIndex;
            std::optional<LayerIndex> layerIndex{std::nullopt};
            bool foundSingleStationSeed{true};
 
            for(const xAOD::MuonSegment* seg : newSeed.segments()) {
                if (!layerIndex) {
                    layerIndex = toLayerIndex(seg->chamberIndex());
                    ATH_MSG_DEBUG("First segment is in layer "<<*layerIndex);
                } else if ( (*layerIndex) != toLayerIndex(seg->chamberIndex())) {
                    ATH_MSG_DEBUG("Found segment in layer "<<toLayerIndex(seg->chamberIndex())<<" which is different from the first segment in layer "<<*layerIndex);
                    foundSingleStationSeed = false;
                    break;
                } 
            }
            if(foundSingleStationSeed) {
                continue;
            }
 
            //Check if we have multiple segments from the same station, if so split the seed and create duplicate seeds
            const double r = newSeed.location() == Location::Barrel ? m_cfg.barrelRadius 
                                                                    : coords[Acts::toUnderlying(ePosOnCylinder)];
            const double z = newSeed.location() == Location::Barrel ? coords[Acts::toUnderlying(ePosOnCylinder)] 
                                                                    : coords[Acts::toUnderlying(eDetSection)]* m_cfg.endcapDiscZ;
 
            Amg::Vector3D pos = r * newSeed.sector().radialDir()
                              + z * Amg::Vector3D::UnitZ();
            
            newSeed.setPosition(std::move(pos));
            ATH_MSG_VERBOSE("Add new seed "<<newSeed);
            trackSeeds.emplace_back(std::move(newSeed));
        }
        ATH_MSG_VERBOSE("Found in total "<<trackSeeds.size()<<" before overlap removal");
        return resolveOverlaps(std::move(trackSeeds));
    }

◆ forceIntegration()

Amg::Vector3D MuonR4::MsTrackSeeder::forceIntegration	(	const PosMomPair_t &	point1,
		const PosMomPair_t &	point2,
		const Amg::Vector3D &	planeNorm,
		MagField::AtlasFieldCache &	fieldCache ) const

private

Compute the integral of magnetic force (v x B ) dS along a trajectory, given the initial and final positions and directions of the trajectory.

The trajectory is approximated as a straight line between the two positions, and the magnetic field is evaluated at several points along this line.

Parameters

point1	Initial position and direction
point2	Final position and direction
planeNorm	Normal vector of the bending plane used for the momentum estimation, used to extract the orthogonal component to the field
fieldCache	Magnetic field cache

Returns: : Integral of magnetic force

Definition at line 185 of file MsTrackSeeder.cxx.

                                                                                             {
        const auto& [pos1, dir1] = point1;
        const auto& [pos2, dir2] = point2;
 
        Amg::Vector3D locField{Amg::Vector3D::Zero()};
        Amg::Vector3D accumForce{Amg::Vector3D::Zero()};
        for (double fieldStep : m_fieldExtpSteps) {
            const Amg::Vector3D extPos {(1. - fieldStep) * pos1 + fieldStep * pos2};
            const Amg::Vector3D extDir {((1. - fieldStep) * dir1 + fieldStep * dir2).unit()};
                    
            fieldCache.getField(extPos.data(), locField.data());
            const Amg::Vector3D locForce {locField.dot(planeNorm) * extDir.cross(planeNorm)};
            accumForce += locForce;
 
            ATH_MSG_VERBOSE(__func__<<"() step: "<<fieldStep
                <<", pos: "<<Amg::toString(extPos)<<", dir: "<<Amg::toString(extDir)
                <<" --> local |B|: "<<locField.mag()*1e3 << " [T]"<<", |Bnorm|: "<<locField.dot(planeNorm)*1e3
                <<" [T], local |v x Bnorm|: "<<locForce.mag()*1e3<<" [T].");
        }
        const double dS {(pos2 - pos1).mag() / static_cast<double>(m_fieldExtpSteps.size())};
        return accumForce * dS;
    }

◆ getPtimesQ()

double MuonR4::MsTrackSeeder::getPtimesQ	(	const Amg::Vector3D &	forceIntegral,
		const Amg::Vector3D &	deltaDir ) const

private

Compute the charge times momentum from the integral of lorentz force and the total change in direction.

Parameters

forceIntegral	Cumulative lorentz force
deltaDir	Change in direction

Returns: : Estimated Q*P value

Definition at line 210 of file MsTrackSeeder.cxx.

                                                                        {
        const double PtimesQ {0.3 * Gaudi::Units::GeV * forceIntegral.mag2() / deltaDir.dot(forceIntegral)};
 
        ATH_MSG_VERBOSE("estimateQtimesP() force integral: "<<forceIntegral.mag()<<" [T*m], deltaDir: "
            <<deltaDir.mag()<<", cos: "<<deltaDir.dot(forceIntegral)/ (deltaDir.mag() * forceIntegral.mag())
            <<", PtimesQ: "<<PtimesQ/Gaudi::Units::GeV <<" [GeV].");
        return PtimesQ;                                    
    }

◆ initMessaging()

void AthMessaging::initMessaging ( ) const

privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

{
  m_imsg = Athena::getMessageSvc();
  // If user did not set an explicit level, set a default
  if (m_lvl == MSG::NIL) {
    m_lvl = m_imsg ?
      static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
      MSG::INFO;
  }
}

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 167 of file AthMessaging.h.

{
  MsgStream* ms = m_msg_tls.get();
  if (!ms) {
    if (!m_initialized.test_and_set()) initMessaging();
    ms = new MsgStream(m_imsg,m_nm);
    m_msg_tls.reset( ms );
  }
 
  ms->setLevel (m_lvl);
  return *ms;
}

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 182 of file AthMessaging.h.

183{ return msg() << lvl; }

AthMessaging::msg

MsgStream & msg() const

The standard message stream.

Definition AthMessaging.h:167

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Test the output level.

Parameters

lvl	The message level to test against

Returns: boolean Indicating if messages at given level will be printed

Return values

true	Messages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

{
  // If user did not set explicit message level we have to initialize
  // the messaging and retrieve the default via the MessageSvc.
  if (m_lvl==MSG::NIL && !m_initialized.test_and_set()) initMessaging();
 
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

◆ resolveOverlaps()

std::unique_ptr< MsTrackSeedContainer > MuonR4::MsTrackSeeder::resolveOverlaps ( MsTrackSeedContainer && unresolved ) const

private

Removes exact duplciates or partial subsets of the MsTrackSeeds.

Parameters

unresolved Input MsTrackSeedContainer with duplicates

Resort the seeds starting from the ones with the most segments to the lowest

Definition at line 443 of file MsTrackSeeder.cxx.

                                                                              {
        std::ranges::sort(unresolved, [](const MsTrackSeed& a, const MsTrackSeed&b) {
            return a.segments().size() > b.segments().size();
        });
        auto outputSeeds = std::make_unique<MsTrackSeedContainer>();
        outputSeeds->reserve(unresolved.size());
        std::ranges::copy_if(std::move(unresolved), std::back_inserter(*outputSeeds),
            [&outputSeeds](const MsTrackSeed& testMe) {
                for (const MsTrackSeed&  good : *outputSeeds){
                    if (!testMe.sector().isNeighbour(good.sector())) {
                        continue;
                    }
                    const std::size_t sharedSegs = std::ranges::count_if(testMe.segments(),
                                                                      [&good](const xAOD::MuonSegment* segInTest){
                                                                          return std::ranges::find(good.segments(), segInTest) != 
                                                                                 good.segments().end();
                                                                      });
                    if (sharedSegs == testMe.segments().size()) {
                        return false;
                    }
                }
                return true;
            });
 
        ATH_MSG_VERBOSE("Found in total "<<outputSeeds->size()<<" after overlap removal");
        return outputSeeds;
    } 

◆ segDirOntoPhiPlane()

Amg::Vector3D MuonR4::MsTrackSeeder::segDirOntoPhiPlane	(	const Amg::Vector3D &	planeNorm,
		const Amg::Vector3D &	dirToProject )

static

Projects the segment direction onto the plane with global phi = x by removing the component orthogonal to the plane.

Parameters

planeNorm	Normal of the phi plane onto which the segment is projected
dirToProject	Direction to project

Definition at line 55 of file MsTrackSeeder.cxx.

                                                                                     {
        return (dirToProject - dirToProject.dot(planeNorm) * planeNorm).unit();                                                
    }

◆ segPosOntoPhiPlane()

Amg::Vector3D MuonR4::MsTrackSeeder::segPosOntoPhiPlane	(	const Amg::Vector3D &	planeNorm,
		const int	Sector,
		const Amg::Vector3D &	posToProject )

static

Projects the segment position onto the plane with global phi = x The local coordinate system is arranged such that the x-axis is co-linear to the phi direction.

The segment is moved along the MDT's wire direction in that sector.

Parameters

planeNorm	Normal of the phi plane onto which the segment is projected
Sector	Sector of the segment to be projected, needed to find the wire direction
posToProject	Position to project

Definition at line 45 of file MsTrackSeeder.cxx.

                                                                                     {
        // Find the sensor direction
        Amg::Vector3D projDir {
            ExpandedSector{static_cast<unsigned>(Sector), 
                           ExpandedSector::SectorProjector::center}.normalDir()};
        return Acts::PlanarHelper::intersectPlane(
            posToProject, projDir, planeNorm, Amg::Vector3D::Zero()).position();                                
    }

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level lvl )

inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

{
  m_lvl = lvl;
}

◆ withinBounds()

bool MuonR4::MsTrackSeeder::withinBounds	(	const Amg::Vector2D &	projPos,
		const Location	loc ) const

Returns whether the expression on the cylinder is within the surface bounds.

Parameters

projPos	Projected position on the cylinder
loc	Surface location: [barrel/endcap]

Definition at line 89 of file MsTrackSeeder.cxx.

                                                               {
        using enum Location;
        if (loc == Barrel && std::abs(projPos[1]) > std::min(m_cfg.endcapDiscZ, m_cfg.barrelLength)) {
            ATH_MSG_VERBOSE(__func__<<"()  "<<__LINE__<<" - Position "<<Amg::toString(projPos)<<
                            " exceeds cylinder boundaries ("<<m_cfg.barrelRadius<<", "
                            <<std::min(m_cfg.endcapDiscZ, m_cfg.barrelLength)<<")");
            return false;
        } else if (loc == Endcap && (0 > projPos[0] || projPos[0] > m_cfg.endcapDiscRadius)) {
            ATH_MSG_VERBOSE(__func__<<"()  "<<__LINE__<<" - Position "<<Amg::toString(projPos)<<
            " exceeds endcap boundaries ("<<m_cfg.endcapDiscRadius<<", "<<Acts::copySign(m_cfg.endcapDiscZ, projPos[1])<<")");
            return false;
        }
        return true;
    }

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT

mutableprivateinherited

Messaging initialized (initMessaging).

Definition at line 141 of file AthMessaging.h.

◆ m_cfg

Config MuonR4::MsTrackSeeder::m_cfg {}

private

Definition at line 186 of file MsTrackSeeder.h.

186{};

◆ m_fieldExtpSteps

std::set<double> MuonR4::MsTrackSeeder::m_fieldExtpSteps {}

private

Definition at line 185 of file MsTrackSeeder.h.

185{};

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

135{ nullptr };

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

138{ MSG::NIL };

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls

mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels).

Definition at line 132 of file AthMessaging.h.

◆ m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

The documentation for this class was generated from the following files:

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ Location

◆ PosMomPair_t

◆ SearchTree_t

◆ SectorProjector

◆ TreeRawVec_t

◆ VecOpt_t

Member Enumeration Documentation

◆ SeedCoords

Constructor & Destructor Documentation

◆ MsTrackSeeder()

Member Function Documentation

◆ appendSegment()

◆ constructTree()

◆ estimateQtimesP() [1/3]

◆ estimateQtimesP() [2/3]

◆ estimateQtimesP() [3/3]

◆ expressOnCylinder()

◆ findTrackSeeds()

◆ forceIntegration()

◆ getPtimesQ()

◆ initMessaging()

◆ msg() [1/2]

◆ msg() [2/2]

◆ msgLvl()

◆ resolveOverlaps()

◆ segDirOntoPhiPlane()

◆ segPosOntoPhiPlane()

◆ setLevel()

◆ withinBounds()

Member Data Documentation

◆ ATLAS_THREAD_SAFE

◆ m_cfg

◆ m_fieldExtpSteps

◆ m_imsg

◆ m_lvl

◆ m_msg_tls

◆ m_nm