MuonR4::MdtSegmentSeedGenerator Class Reference

Helper class to generate valid seeds for the segment fit. More...

#include <MdtSegmentSeedGenerator.h>

Inheritance diagram for MuonR4::MdtSegmentSeedGenerator:

Collaboration diagram for MuonR4::MdtSegmentSeedGenerator:

Classes
struct	Config
	Configuration switches of the module More...
struct	DriftCircleSeed
	Helper struct from a generated Mdt seed. More...
struct	SeedFitAuxilliaries
	Auxillary struct to calculate fit constants. More...
struct	SeedFitAuxWithT0
struct	SeedSolution
	Cache of all solutions seen thus far. More...

Public Types
using	HitVec = SpacePointPerLayerSplitter::HitVec

Public Member Functions
	MdtSegmentSeedGenerator (const std::string &name, const SegmentSeed *segmentSeed, const Config &configuration)
	Standard constructor taking the segmentSeed to start with and then few configuration tunes. More...
	~MdtSegmentSeedGenerator ()
std::optional< DriftCircleSeed >	nextSeed (const EventContext &ctx)
	returns the next seed in the row More...
unsigned int	numGenerated () const
	Returns how many seeds have been generated. More...
const Config &	config () const
	Returns the current seed configuration. More...
bool	msgLvl (const MSG::Level lvl) const
	Test the output level. More...
MsgStream &	msg () const
	The standard message stream. More...
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream. More...
void	setLevel (MSG::Level lvl)
	Change the current logging level. More...

Private Types
using	SignComboType = std::array< int, 2 >
	Sign combinations to draw the 4 lines tangent to 2 drift circles The first two are indicating whether the tangent is left/right to the first/second circle. More...

Private Member Functions
std::optional< DriftCircleSeed >	buildSeed (const EventContext &ctx, const HoughHitType &topHit, const HoughHitType &bottomHit, const SignComboType &signs)
	Tries to build the seed from the two hits. More...
SeedFitAuxilliaries	estimateAuxillaries (const DriftCircleSeed &seed) const
	Helper function to estimate the auxillary variables that remain constant during the fit. More...
SeedFitAuxWithT0	estimateAuxillaries (const EventContext &ctx, const DriftCircleSeed &seed) const
	Helper function to estimate the auxillary variables that remain constants during the fit with t0. More...
void	fitDriftCircles (DriftCircleSeed &seed) const
	Refine the seed by performing a fast Mdt segment fit. More...
void	fitDriftCirclesWithT0 (const EventContext &ctx, DriftCircleSeed &seed) const
	Refine the seed by performing a fast Mdt segment fit with t0 constraint. More...
void	moveToNextCandidate ()
	Prepares the generator to generate the seed from the next pair of drift circles. More...
void	initMessaging () const
	Initialize our message level and MessageSvc. More...

Private Attributes
Config	m_cfg {}
const SegmentSeed *	m_segmentSeed {nullptr}
SpacePointPerLayerSplitter	m_hitLayers {m_segmentSeed->getHitsInMax()}
std::size_t	m_upperLayer {0}
	Considered layer to pick the top drift circle from. More...
std::size_t	m_lowerLayer {0}
	Considered layer to pick the bottom drift circle from. More...
std::size_t	m_lowerHitIndex {0}
	Explicit hit to pick in the selected bottom layer. More...
std::size_t	m_upperHitIndex {0}
	Explicit hit to pick in the selected top layer. More...
std::size_t	m_signComboIndex {0}
	Index of the left-right ambiguity between the circles. More...
std::vector< SeedSolution >	m_seenSolutions {}
	Vector caching equivalent solutions to avoid double seeding. More...
unsigned int	m_nGenSeeds {0}
	Counter on how many seeds have been generated. More...
std::string	m_nm
	Message source name. More...
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels) More...
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer. More...
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level. More...
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging) More...

Static Private Attributes
constexpr static std::array< SignComboType, 4 >	s_signCombos

Detailed Description

Helper class to generate valid seeds for the segment fit.

The generator first returns a seed directly made from the patten recogntion. Afterwards it builds seeds by lying tangent lines to a pair of drift circles. The pairing starts from the innermost & outermost layers with tubes. A valid seed must have at least 4 associated hits which are within a chi2 of 5. If two seeds within the parameter resolution are generated, then the latter one is skipped.

Definition at line 24 of file MdtSegmentSeedGenerator.h.

Member Typedef Documentation

HitVec

using MuonR4::MdtSegmentSeedGenerator::HitVec = SpacePointPerLayerSplitter::HitVec

Definition at line 26 of file MdtSegmentSeedGenerator.h.

SignComboType

using MuonR4::MdtSegmentSeedGenerator::SignComboType = std::array<int, 2>

private

Sign combinations to draw the 4 lines tangent to 2 drift circles The first two are indicating whether the tangent is left/right to the first/second circle.

The last sign is picking the sign of the solution arising from the final quadratic equation.

Definition at line 101 of file MdtSegmentSeedGenerator.h.

Constructor & Destructor Documentation

MdtSegmentSeedGenerator()

MuonR4::MdtSegmentSeedGenerator::MdtSegmentSeedGenerator	(	const std::string &	name,
		const SegmentSeed *	segmentSeed,
		const Config &	configuration
	)

Standard constructor taking the segmentSeed to start with and then few configuration tunes.

Parameters

name	Name of the Seed generator's logger
segmentSeed	Seed from which the seeds for the fit are built
configuration	Passed configuration settings of the generator

Check whether the first layer is too busy

Check whether the lower layer is too busy

Definition at line 41 of file MdtSegmentSeedGenerator.cxx.

                                                                                 :
            AthMessaging{name},
            m_cfg{configuration},
            m_segmentSeed{segmentSeed} {
 
        if (m_hitLayers.mdtHits().empty()) return;
        
        if (std::ranges::find_if(m_hitLayers.mdtHits(), [this](const HitVec& vec){
            return vec.size() > m_cfg.busyLayerLimit;
        }) != m_hitLayers.mdtHits().end()) {
            m_cfg.startWithPattern = false;
        }
        // Set the start for the upper layer
        m_upperLayer = m_hitLayers.mdtHits().size()-1; 
 
        while (m_lowerLayer < m_upperLayer && m_hitLayers.mdtHits()[m_lowerLayer].size() >m_cfg.busyLayerLimit){
            ++m_lowerLayer;
        }
        while (m_lowerLayer < m_upperLayer && m_hitLayers.mdtHits()[m_upperLayer].size() > m_cfg.busyLayerLimit) {
            --m_upperLayer;
        }
 
        if (msgLvl(MSG::VERBOSE)) {
            std::stringstream sstr{};
            for (const auto [layCount, layer] : Acts::enumerate(m_hitLayers.mdtHits())) { 
                sstr<<"Mdt-hits in layer "<<layCount<<": "<<layer.size()<<std::endl;
                for (const HoughHitType& hit : layer) {
                    sstr<<"   **** "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<" "
                        <<Amg::toString(hit->positionInChamber())<<", driftRadius: "<<hit->driftRadius()<<std::endl;
                }
            }
            for (const auto [layCount, layer] : Acts::enumerate(m_hitLayers.stripHits())) { 
                sstr<<"Hits in layer "<<layCount<<": "<<layer.size()<<std::endl;
                for (const HoughHitType& hit : layer) {
                    sstr<<"   **** "<<hit->msSector()->idHelperSvc()->toString(hit->identify())<<" "
                        <<Amg::toString(hit->positionInChamber())<<", driftRadius: "<<hit->driftRadius()<<std::endl;
                }
            }
            ATH_MSG_VERBOSE("SeedGenerator - sorting of hits done. Mdt layers: "<<m_hitLayers.mdtHits().size()
                            <<", strip layers: "<<m_hitLayers.stripHits().size()<<std::endl<<sstr.str()<<std::endl<<std::endl);
        }
    }

~MdtSegmentSeedGenerator()

MuonR4::MdtSegmentSeedGenerator::~MdtSegmentSeedGenerator ( )

default

Member Function Documentation

buildSeed()

std::optional< MdtSegmentSeedGenerator::DriftCircleSeed > MuonR4::MdtSegmentSeedGenerator::buildSeed ( const EventContext &  ctx, const HoughHitType &  topHit, const HoughHitType &  bottomHit, const SignComboType &  signs  )

private

Tries to build the seed from the two hits.

Fails if the solution is invalid or if the seed has already been built before

Parameters

topHit	Hit candidate from the upper layer
bottomHit	Hit candidate from the lower layer
sign	Object encoding whether the tangent is left / right

All other sign combinations will also get stuck -> force the measurement selector to get another combo

Create a new line position & direction which also takes the potential phi estimates into account

Recalculate the seed with the calibrated parameters

Check that the is within the predefined window

d/dx asin(x) = 1 / sqrt(1- x*x)

Collect all hits close to the seed line

what ever comes after is not matching onto the segment

Reject seeds with too litle Mdt hit association

Last check wheather another seed with the same left-right combination hasn't already been found

All seed hits are of the same size

Including the places where no seed hit was assigned. Both solutions match in terms of left-right solutions. It's very likely that they're converging to the same segment.

Seed candidate is

Add the solution to the list. That we don't iterate twice over it

If we found a long Mdt seed, then ensure that all subsequent seeds have at least the same amount of Mdt hits.

Associate strip hits

Add all hits with a pull better than the threshold

Definition at line 172 of file MdtSegmentSeedGenerator.cxx.

                                                                       {
        
        const auto* bottomPrd = static_cast<const xAOD::MdtDriftCircle*>(bottomHit->primaryMeasurement()); 
        const auto* topPrd = static_cast<const xAOD::MdtDriftCircle*>(topHit->primaryMeasurement());
        if (bottomPrd->status() != Muon::MdtDriftCircleStatus::MdtStatusDriftTime ||
            topPrd->status() != Muon::MdtDriftCircleStatus::MdtStatusDriftTime) {
                m_signComboIndex = s_signCombos.size();
                return std::nullopt;
        }
        const auto&[signTop, signBot] = signs;
        double R = signBot *bottomHit->driftRadius() - signTop * topHit->driftRadius(); 
        const Amg::Vector3D& bottomPos{bottomHit->positionInChamber()};
        const Amg::Vector3D& topPos{topHit->positionInChamber()};
        const Muon::IMuonIdHelperSvc* idHelperSvc{topHit->msSector()->idHelperSvc()};
        const Amg::Vector3D D = topPos - bottomPos;
        const double thetaTubes = std::atan2(D.y(), D.z()); 
        const double distTubes =  std::hypot(D.y(), D.z());
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Bottom tube "<<idHelperSvc->toString(bottomHit->identify())<<" "<<Amg::toString(bottomPos)
                    <<", driftRadius: "<<bottomHit->driftRadius()<<" - top tube "<<idHelperSvc->toString(topHit->identify())
                    <<" "<<Amg::toString(topPos)<<", driftRadius: "<<topHit->driftRadius()
                    <<", tube distance: "<<Amg::toString(D));
 
        DriftCircleSeed candidateSeed{};
        candidateSeed.parameters = m_segmentSeed->parameters();
        candidateSeed.parentBucket = m_segmentSeed->parentBucket();
        double theta{thetaTubes - std::asin(std::clamp(R / distTubes, -1., 1.))};
        Amg::Vector3D seedDir = Amg::dirFromAngles(90.*Gaudi::Units::deg, theta);
        double Y0 = bottomPos.y()*seedDir.z() - bottomPos.z()*seedDir.y() + signBot*bottomHit->driftRadius();
        double combDriftUncert{std::sqrt(bottomPrd->driftRadiusCov() + topPrd->driftRadiusCov())};
        std::unique_ptr<CalibratedSpacePoint> calibBottom{}, calibTop{};
        if (m_cfg.recalibSeedCircles) {
            candidateSeed.parameters[toInt(ParamDefs::theta)] = theta;
            candidateSeed.parameters[toInt(ParamDefs::y0)] = Y0 / seedDir.z();
            const auto [linePos, lineDir] = makeLine(candidateSeed.parameters);
            calibBottom = m_cfg.calibrator->calibrate(ctx, bottomHit, linePos, lineDir, 
                                                      candidateSeed.parameters[toInt(ParamDefs::time)]);
            calibTop = m_cfg.calibrator->calibrate(ctx, topHit, linePos, lineDir, 
                                                   candidateSeed.parameters[toInt(ParamDefs::time)]);
            R = signBot * calibBottom->driftRadius() - signTop * calibTop->driftRadius();
            theta =  thetaTubes - std::asin(std::clamp(R / distTubes, -1., 1.));
            seedDir = Amg::dirFromAngles(90.*Gaudi::Units::deg, theta);
            Y0 = bottomPos.y()*seedDir.z() - bottomPos.z()*seedDir.y() + signBot*bottomHit->driftRadius();
            combDriftUncert = std::sqrt(driftCov(*calibBottom) + driftCov(*calibTop));
        }
        
        candidateSeed.parameters[toInt(ParamDefs::theta)] = theta;
        candidateSeed.parameters[toInt(ParamDefs::y0)] = Y0 / seedDir.z();
        if (!passRangeCut(m_cfg.thetaRange, theta) || 
            !passRangeCut(m_cfg.interceptRange, candidateSeed.parameters[toInt(ParamDefs::y0)])) {
            return std::nullopt;
        }
 
        const Amg::Vector3D seedPos = Y0 / seedDir.z() * Amg::Vector3D::UnitY();
 
        assert(std::abs(topPos.y()*seedDir.z() - topPos.z() * seedDir.y() + signTop*topHit->driftRadius() - Y0) < std::numeric_limits<float>::epsilon() );
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Candidate seed theta: "<<theta<<", tanTheta: "<<(seedDir.y() / seedDir.z())<<", y0: "<<Y0/seedDir.z());
 
        SeedSolution solCandidate{};
        solCandidate.Y0 = seedPos[toInt(ParamDefs::y0)];
        solCandidate.theta = theta;
        const double denomSquare =  1. - std::pow(R / distTubes, 2); 
        if (denomSquare < std::numeric_limits<double>::epsilon()){
            ATH_MSG_VERBOSE("Invalid seed, rejecting"); 
            return std::nullopt; 
        }
        solCandidate.dTheta =  combDriftUncert / std::sqrt(denomSquare) / distTubes;
        solCandidate.dY0 =  std::hypot(-bottomPos.y()*seedDir.y() + bottomPos.z()*seedDir.z(), 1.) * solCandidate.dTheta;
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Test new "<<solCandidate<<". "<<m_seenSolutions.size());
 
 
        if (std::ranges::find_if(m_seenSolutions,
                        [&solCandidate, this] (const SeedSolution& seen) {
                            const double deltaY = std::abs(seen.Y0 - solCandidate.Y0);
                            const double limitY = std::hypot(seen.dY0, solCandidate.dY0);
                            const double dTheta = std::abs(seen.theta - solCandidate.theta);
                            const double limitTh = std::hypot(seen.dTheta, solCandidate.dTheta);
                            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": "<<seen
                                    <<std::format(" delta Y: {:.2f} {:} {:.2f}", deltaY, deltaY < limitY ? '<' : '>', limitY)
                                    <<std::format(" delta theta: {:.2f} {:} {:.2f}", dTheta, dTheta < limitTh ? '<' : '>', limitTh) );
                            return deltaY < limitY && dTheta < limitTh;;
                        }) != m_seenSolutions.end()){
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Reject due to similarity");
            return std::nullopt;
        }
        for (const auto [layerNr,  hitsInLayer] : Acts::enumerate(m_hitLayers.mdtHits())) {
            ATH_MSG_VERBOSE( __func__<<"() "<<__LINE__<<": "<<hitsInLayer.size()<<" hits in layer "<<(layerNr +1));
            bool hadGoodHit{false};
            for (const HoughHitType testMe : hitsInLayer){
                const double pull = std::sqrt(SegmentFitHelpers::chiSqTermMdt(seedPos, seedDir, *testMe, msg()));            
                ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Test hit "<<idHelperSvc->toString(testMe->identify())
                            <<" "<<Amg::toString(testMe->positionInChamber())<<", pull: "<<pull);              
                if (pull < m_cfg.hitPullCut) {
                    hadGoodHit = true;
                    solCandidate.seedHits.emplace_back(testMe);
                    ++candidateSeed.nMdt;
                }
                else if (hadGoodHit) {
                    break;
                } 
            }
        }
        const unsigned hitCut = std::max(1.*m_cfg.nMdtHitCut, m_cfg.nMdtLayHitCut * m_hitLayers.mdtHits().size()); 
 
        if (1.*candidateSeed.nMdt < hitCut) {
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Too few hits associated "<<candidateSeed.nMdt<<", expect: "<<hitCut<<" hits.");
            return std::nullopt;
        }
        /* Calculate the left-right signs of the used hits */
        if (m_cfg.overlapCorridor) {
            solCandidate.solutionSigns = driftSigns(seedPos, seedDir, solCandidate.seedHits, msg());
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Circle solutions for seed "
                          <<idHelperSvc->toStringChamber(bottomHit->identify())<<" - "<<solCandidate);
            for (unsigned int a = m_cfg.startWithPattern; a< m_seenSolutions.size() ;++a) { 
                const SeedSolution& accepted = m_seenSolutions[a];
                unsigned int nOverlap{0};
                std::vector<int> corridor = driftSigns(seedPos, seedDir, accepted.seedHits,  msg());                
                ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Test seed against accepted "<<accepted<<", updated signs: "<<corridor);
                for (unsigned int l = 0; l < accepted.seedHits.size(); ++l){
                    nOverlap  += corridor[l] == accepted.solutionSigns[l];
                }
                if (nOverlap == corridor.size() && accepted.seedHits.size() >= solCandidate.seedHits.size()) {
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Same set of hits collected within the same corridor");
                    return std::nullopt;
                }
            }
        }
        for (const HoughHitType& hit : solCandidate.seedHits){
            //calibBottom is nullptr after it has been moved, so...
            //cppcheck-suppress accessMoved
            if (hit == bottomHit && calibBottom) {
                candidateSeed.measurements.emplace_back(std::move(calibBottom));
            } 
            //calibTop is nullptr after it has been moved, so...
            //cppcheck-suppress accessMoved
            else if (hit == topHit && calibTop) {
                candidateSeed.measurements.emplace_back(std::move(calibTop));
            } else {
                candidateSeed.measurements.emplace_back(m_cfg.calibrator->calibrate(ctx, hit, seedPos, seedDir, 0.));
            }
        }
        m_seenSolutions.emplace_back(std::move(solCandidate));
        if (m_cfg.tightenHitCut) {
            m_cfg.nMdtHitCut = std::max(m_cfg.nMdtHitCut, candidateSeed.nMdt);
        }
        ++m_nGenSeeds;        
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": In event "<<ctx.eventID().event_number()<<" found new seed solution "<<toString(candidateSeed.parameters));
        if (m_cfg.fastSeedFit) {
            if (!m_cfg.fastSegFitWithT0) {
                fitDriftCircles(candidateSeed);
            } else {
                fitDriftCirclesWithT0(ctx, candidateSeed);
            }
        }
 
        // Combine the seed with the phi estimate
        {
            const Amg::Vector3D parDir = dirFromTangents(m_segmentSeed->tanPhi(), std::tan(theta));
            candidateSeed.parameters[toInt(ParamDefs::theta)] = parDir.theta();
            candidateSeed.parameters[toInt(ParamDefs::phi)] = parDir.phi();
        }
 
        {
            const auto [seedPos, seedDir] = makeLine(candidateSeed.parameters); 
            for (const std::vector<HoughHitType>& hitsInLayer : m_hitLayers.stripHits()) {
                HoughHitType bestHit{nullptr};
                double bestPull{m_cfg.hitPullCut};
                for (const HoughHitType testMe : hitsInLayer){
                    const double pull = std::sqrt(SegmentFitHelpers::chiSqTermStrip(seedPos, seedDir, *testMe, msg())) 
                                      / testMe->dimension();
                    ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<": Test hit "<<idHelperSvc->toString(testMe->identify())
                                <<" "<<Amg::toString(testMe->positionInChamber())<<", pull: "<<pull);
                    if (pull <= bestPull) {
                        bestHit = testMe;
                        bestPull = pull;
                    }
                }
                if (!bestHit) {
                    continue;
                }
                candidateSeed.measurements.push_back(m_cfg.calibrator->calibrate(ctx,bestHit, seedPos, seedDir, 0));           
            }
        }
        return candidateSeed;
    }

config()

const MdtSegmentSeedGenerator::Config & MuonR4::MdtSegmentSeedGenerator::config

(

)

const

Returns the current seed configuration.

Definition at line 37 of file MdtSegmentSeedGenerator.cxx.

                                                                               {
        return m_cfg;
    }

estimateAuxillaries() [1/2]

MdtSegmentSeedGenerator::SeedFitAuxilliaries MuonR4::MdtSegmentSeedGenerator::estimateAuxillaries ( const DriftCircleSeed &  seed ) const

inlineprivate

Helper function to estimate the auxillary variables that remain constant during the fit.

Parameters

seed	Reference to the seed to calculate the variables from

Seed direction vector

y0Prime = y0 * cos(theta)

Calculate the centre of gravity

Calculate the fit constants

Definition at line 378 of file MdtSegmentSeedGenerator.cxx.

                                                                                      {
 
        SeedFitAuxilliaries aux{};
        const Amg::Vector3D seedDir = Amg::dirFromAngles(90.*Gaudi::Units::deg, 
                                                         seed.parameters[toInt(ParamDefs::theta)]);
        const double y0 = seed.parameters[toInt(ParamDefs::y0)] * seedDir.z();
 
        aux.invCovs.reserve(seed.measurements.size());
        aux.driftSigns.reserve(seed.measurements.size());
        double norm{0.};
        for (const std::unique_ptr<CalibratedSpacePoint>& hit : seed.measurements) {
            const double invCov = 1./ driftCov(*hit);
            const Amg::Vector3D& pos{hit->positionInChamber()};
            const int sign = y0  - pos.y() * seedDir.z() + pos.z()* seedDir.y() > 0 ? 1 : -1;
 
            aux.centerOfGrav+= invCov * pos;
            aux.invCovs.push_back(invCov);
            aux.driftSigns.push_back(sign);
 
            norm += invCov;
        }
        aux.covNorm = 1./ norm;
        aux.centerOfGrav *= aux.covNorm;
        for (const auto [covIdx, hit] : Acts::enumerate(seed.measurements)) {
            const double& invCov = aux.invCovs[covIdx];
            const int& sign = aux.driftSigns[covIdx];
            const Amg::Vector3D pos = hit->positionInChamber() - aux.centerOfGrav;
            const double signedCov = invCov * sign;
            aux.T_zzyy += invCov * (std::pow(pos.z(), 2) - std::pow(pos.y(), 2));
            aux.T_yz   += invCov * pos.y()*pos.z();
            aux.T_rz   += signedCov * pos.z() * hit->driftRadius();
            aux.T_ry   += signedCov * pos.y() * hit->driftRadius();
            aux.fitY0  += signedCov * aux.covNorm * hit->driftRadius();
        }
        ATH_MSG_VERBOSE("Estimated T_zzyy: "<<aux.T_zzyy<<", T_yz: "<<aux.T_yz<<", T_rz: "<<aux.T_rz
                     <<", T_ry: "<<aux.T_ry<<", centre "<<Amg::toString(aux.centerOfGrav)<<", y0: "<<aux.fitY0
                     <<", norm: "<<aux.covNorm<<"/"<<norm);
        return aux;
    }

estimateAuxillaries() [2/2]

MdtSegmentSeedGenerator::SeedFitAuxWithT0 MuonR4::MdtSegmentSeedGenerator::estimateAuxillaries ( const EventContext &  ctx, const DriftCircleSeed &  seed  ) const

inlineprivate

Helper function to estimate the auxillary variables that remain constants during the fit with t0.

Parameters

ctx	EventContext to recalibrate the hits
seed	Reference to the seed to calculate the variables from

Definition at line 471 of file MdtSegmentSeedGenerator.cxx.

                                                                                        {
        SeedFitAuxWithT0 aux{estimateAuxillaries(seed)};
        for (const auto [idx, hit] : Acts::enumerate(seed.measurements)){
            const double signedCov = aux.driftSigns[idx] * aux.invCovs[idx];
            const double weight = aux.covNorm * signedCov;
            const double velocity = m_cfg.calibrator->driftVelocity(ctx, *hit); 
            const double acceleration = m_cfg.calibrator->driftAcceleration(ctx, *hit);
            const Amg::Vector3D pos = hit->positionInChamber() - aux.centerOfGrav;
            aux.fitY0Prime+= weight * velocity;
            aux.fitY0TwoPrime+= weight * acceleration;
 
            aux.T_vz += signedCov * pos.z()*velocity;
            aux.T_vy += signedCov * pos.y()*velocity;
            aux.T_az += signedCov * pos.z()*acceleration;
            aux.T_ay += signedCov * pos.y()*acceleration;
 
            aux.R_vr += signedCov * hit->driftRadius() * velocity;
            aux.R_va += signedCov * hit->driftRadius() * acceleration;
            aux.R_vv += signedCov * velocity * velocity;
        }
        ATH_MSG_VERBOSE("Estimated T_vz: "<<aux.T_vz<<", T_vy: "<<aux.T_vy
                     <<", T_az: "<<aux.T_az<<", T_ay: "<<aux.T_ay<<" --- R_vr: "<<aux.R_vr
                     <<", R_va: "<<aux.R_va<<", R_vv: "<<aux.R_vv<<" -- Y0^{'}: "<<aux.fitY0Prime
                     <<", Y0^{''}: "<<aux.fitY0TwoPrime);
        return aux;
    }

fitDriftCircles()

void MuonR4::MdtSegmentSeedGenerator::fitDriftCircles ( DriftCircleSeed &  seed ) const

private

Refine the seed by performing a fast Mdt segment fit.

Parameters

seed	Seed built from the tangent adjacent to the two seed circles

Now it's time to use the guestimate

Definition at line 422 of file MdtSegmentSeedGenerator.cxx.

                                                                               {
 
        const SeedFitAuxilliaries auxVars = estimateAuxillaries(inSeed);
        
        double theta = inSeed.parameters[toInt(ParamDefs::theta)];
        const double thetaGuess = std::atan2( 2.*(auxVars.T_yz - auxVars.T_rz), auxVars.T_zzyy) / 2.;
 
        ATH_MSG_VERBOSE("Start fast fit seed: "<<theta<<", guess: "<<thetaGuess
                    <<", y0: "<<inSeed.parameters[toInt(ParamDefs::y0)]
                    <<", fitY0: "<<auxVars.fitY0<<", centre: "<<Amg::toString(auxVars.centerOfGrav));
        theta = thetaGuess;
        CxxUtils::sincos thetaCS{theta};
        bool converged{false};
        while (!converged && inSeed.nIter++ <= m_cfg.nMaxIter) {
            const CxxUtils::sincos twoTheta{2.*theta};
            const double thetaPrime = 0.5*auxVars.T_zzyy *twoTheta.sn - auxVars.T_yz * twoTheta.cs 
                                    - auxVars.T_rz * thetaCS.cs - auxVars.T_ry * thetaCS.sn;
            if (std::abs(thetaPrime) < m_cfg.precCutOff){
                converged = true;
                break;
            }
 
            const double thetaTwoPrime =  auxVars.T_zzyy * twoTheta.cs + 2.* auxVars.T_yz * twoTheta.sn 
                                       + auxVars.T_rz * thetaCS.sn - auxVars.T_ry * thetaCS.cs;
            const double update = thetaPrime / thetaTwoPrime;
            ATH_MSG_VERBOSE("Fit iteration #"<<inSeed.nIter<<" -- theta: "<<theta<<", thetaPrime: "<<thetaPrime
                        <<", thetaTwoPrime: "<<thetaTwoPrime<<" -- "<<std::format("{:.8f}", update)
                        <<" --> next theta "<<(theta - thetaPrime / thetaTwoPrime));
 
            if (std::abs(update) < m_cfg.precCutOff) {
                converged = true;
                break;
            }
            theta -= update;
            thetaCS = CxxUtils::sincos{theta};
        }
        if (!converged) {
           return;
        }
        double fitY0 = (auxVars.centerOfGrav.y() *thetaCS.cs - auxVars.centerOfGrav.z() * thetaCS.sn + auxVars.fitY0) / thetaCS.cs;
        ATH_MSG_VERBOSE("Drift circle fit converged within "<<inSeed.nIter<<" iterations giving "<<toString(inSeed.parameters)<<", chi2: "<<inSeed.chi2
                        <<" - theta: "<<theta / Gaudi::Units::deg<<", y0: "<<fitY0);
        inSeed.parameters[toInt(ParamDefs::theta)] = theta;
        inSeed.parameters[toInt(ParamDefs::y0)] = fitY0;
    }

fitDriftCirclesWithT0()

void MuonR4::MdtSegmentSeedGenerator::fitDriftCirclesWithT0 ( const EventContext &  ctx, DriftCircleSeed &  seed  ) const

private

Refine the seed by performing a fast Mdt segment fit with t0 constraint.

Parameters

ctx	EventContext to recalibrate the hits
seed	Seed built from the tangent adjacent to the two seed circles

d^{2}chi^{2} / d^{2}theta

Definition at line 499 of file MdtSegmentSeedGenerator.cxx.

                                                                                                              {
        SeedFitAuxWithT0 auxVars{estimateAuxillaries(ctx, inSeed)};
        
 
        bool converged{false};
        AmgSymMatrix(2) cov{AmgSymMatrix(2)::Zero()};
        AmgVector(2) grad{AmgVector(2)::Zero()};
        AmgVector(2) pars{inSeed.parameters[toInt(ParamDefs::theta)], 
                          inSeed.parameters[toInt(ParamDefs::time)]};
        
 
        while (!converged && inSeed.nIter++ <= m_cfg.nMaxIter) {
 
            const CxxUtils::sincos thetaCS{pars[0]};
            const CxxUtils::sincos twoTheta{2.*pars[0]};
            cov(0,0) = auxVars.T_zzyy * twoTheta.cs + 2.* auxVars.T_yz * twoTheta.sn 
                     + auxVars.T_rz * thetaCS.sn - auxVars.T_ry * thetaCS.cs;
 
            // cov(1,0) =  cov(0,1) = auxVars.T_vz * thetaCS.cs + auxVars.T_vy * thetaCS.sn;
            cov(1,1) = - auxVars.fitY0Prime *auxVars.fitY0Prime  - auxVars.fitY0Prime * auxVars.fitY0TwoPrime 
                     - auxVars.T_az * thetaCS.sn - auxVars.T_ay * thetaCS.cs + auxVars.R_vv + auxVars.R_va;
 
            grad[0] =0.5*auxVars.T_zzyy *twoTheta.sn - auxVars.T_yz * twoTheta.cs 
                                    - auxVars.T_rz * thetaCS.cs - auxVars.T_ry * thetaCS.sn;
            grad[1] = auxVars.fitY0 * auxVars.fitY0Prime - auxVars.R_vr  + auxVars.T_vz * thetaCS.sn - auxVars.T_vy * thetaCS.cs;
 
            
            const AmgVector(2) update = cov.inverse()* grad;
            ATH_MSG_VERBOSE("Iteration: "<<inSeed.nIter<<", theta: "<<pars[0] / Gaudi::Units::deg<<", time: "<< 
                        pars[1]<<" gradient: ("<<(grad[0])<<", "<<grad[1]<<"), covariance:" 
                        <<std::endl<<Amg::toString(cov)<<std::endl<<" update: ("<<update[0] / Gaudi::Units::deg
                        <<", "<<update[1]<<").");
            pars -= update;
        }
 
    }

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();
  m_lvl = m_imsg ?
    static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
    MSG::INFO;
}

moveToNextCandidate()

void MuonR4::MdtSegmentSeedGenerator::moveToNextCandidate ( )

inlineprivate

Prepares the generator to generate the seed from the next pair of drift circles.

Vary the left-right solutions

Move to the next hit in the lower layer

Move to the next hit in the upper layer

All combinations of hits & lines in both layers are processed Switch to the next lowerLayer. But skip the busy ones according to the configuration

Abort the loop if we parsed the multi-layer boundary

Definition at line 91 of file MdtSegmentSeedGenerator.cxx.

                                                             {
        const HitVec& lower = m_hitLayers.mdtHits()[m_lowerLayer];
        const HitVec& upper = m_hitLayers.mdtHits()[m_upperLayer];
        if (++m_signComboIndex < s_signCombos.size()) {
            return;
        }
        m_signComboIndex = 0; 
        
         if (++m_lowerHitIndex < lower.size()) {
            return;
        }
        m_lowerHitIndex=0;
        if (++m_upperHitIndex < upper.size()) {
            return;
        }
         m_upperHitIndex = 0;
        while (m_lowerLayer < m_upperLayer && m_hitLayers.mdtHits()[++m_lowerLayer].size() > m_cfg.busyLayerLimit) {
        } 
        
        if (m_lowerLayer < m_upperLayer) {
            return;
        }
        if (m_lowerLayer >= m_hitLayers.firstLayerFrom2ndMl() && numGenerated()){
            m_lowerLayer = m_upperLayer;
            return;
        }
        m_lowerLayer = 0; 
        while (m_lowerLayer < m_upperLayer && m_hitLayers.mdtHits()[--m_upperLayer].size() > m_cfg.busyLayerLimit){
    
        }
    }

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 164 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 179 of file AthMessaging.h.

{ return msg() << lvl; }

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 (!m_initialized.test_and_set()) initMessaging();
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

nextSeed()

std::optional< MdtSegmentSeedGenerator::DriftCircleSeed > MuonR4::MdtSegmentSeedGenerator::nextSeed

(

const EventContext &

ctx

)

returns the next seed in the row

Increment for the next candidate

If a candidate is built return it. Otherwise continue the process

Definition at line 129 of file MdtSegmentSeedGenerator.cxx.

                                                                 {
        std::optional<DriftCircleSeed> found = std::nullopt; 
        if (!m_nGenSeeds && m_cfg.startWithPattern) {
            ++m_nGenSeeds;
            found = std::make_optional<DriftCircleSeed>();
            found->parameters = m_segmentSeed->parameters();
            found->measurements = m_cfg.calibrator->calibrate(ctx,
                                                              m_segmentSeed->getHitsInMax(),
                                                              m_segmentSeed->positionInChamber(),
                                                              m_segmentSeed->directionInChamber(),0.);
            found->parentBucket = m_segmentSeed->parentBucket();
            found->nMdt = std::ranges::count_if(m_segmentSeed->getHitsInMax(),
                                                [](const SpacePoint* hit){
                                                    return hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType;
                                                });
            SeedSolution patternSeed{};
            patternSeed.seedHits.resize(2*m_hitLayers.mdtHits().size());
            patternSeed.solutionSigns.resize(2*m_hitLayers.mdtHits().size());
            patternSeed.Y0 = m_segmentSeed->interceptY();
            patternSeed.theta = m_segmentSeed->tanTheta();
            m_seenSolutions.push_back(std::move(patternSeed));
            return found;
        }
       
        while (m_lowerLayer < m_upperLayer) {
            const HitVec& lower = m_hitLayers.mdtHits().at(m_lowerLayer);
            const HitVec& upper = m_hitLayers.mdtHits().at(m_upperLayer);
            ATH_MSG_VERBOSE("Layers with hits: "<<m_hitLayers.mdtHits().size()
                            <<" -- next bottom hit: "<<m_lowerLayer<<", hit: "<<m_lowerHitIndex
                            <<" ("<<lower.size()<<"), topHit " <<m_upperLayer<<", "<<m_upperHitIndex
                            <<" ("<<upper.size()<<") - ambiguity "<<s_signCombos[m_signComboIndex]);
 
            found = buildSeed(ctx, upper.at(m_upperHitIndex), lower.at(m_lowerHitIndex), s_signCombos.at(m_signComboIndex));
            moveToNextCandidate();
            if (found) {
                return found;
            }
        }
        return std::nullopt; 
    }

numGenerated()

unsigned int MuonR4::MdtSegmentSeedGenerator::numGenerated

(

)

const

Returns how many seeds have been generated.

Definition at line 88 of file MdtSegmentSeedGenerator.cxx.

                                                             {
        return m_nGenSeeds;
    }

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;
}

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::MdtSegmentSeedGenerator::m_cfg {}

private

Definition at line 204 of file MdtSegmentSeedGenerator.h.

m_hitLayers

SpacePointPerLayerSplitter MuonR4::MdtSegmentSeedGenerator::m_hitLayers {m_segmentSeed->getHitsInMax()}

private

Definition at line 208 of file MdtSegmentSeedGenerator.h.

m_imsg

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

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

m_lowerHitIndex

std::size_t MuonR4::MdtSegmentSeedGenerator::m_lowerHitIndex {0}

private

Explicit hit to pick in the selected bottom layer.

Definition at line 215 of file MdtSegmentSeedGenerator.h.

m_lowerLayer

std::size_t MuonR4::MdtSegmentSeedGenerator::m_lowerLayer {0}

private

Considered layer to pick the bottom drift circle from.

Definition at line 213 of file MdtSegmentSeedGenerator.h.

m_lvl

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

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

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_nGenSeeds

unsigned int MuonR4::MdtSegmentSeedGenerator::m_nGenSeeds {0}

private

Counter on how many seeds have been generated.

Definition at line 223 of file MdtSegmentSeedGenerator.h.

m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

m_seenSolutions

std::vector<SeedSolution> MuonR4::MdtSegmentSeedGenerator::m_seenSolutions {}

private

Vector caching equivalent solutions to avoid double seeding.

Definition at line 221 of file MdtSegmentSeedGenerator.h.

m_segmentSeed

const SegmentSeed* MuonR4::MdtSegmentSeedGenerator::m_segmentSeed {nullptr}

private

Definition at line 207 of file MdtSegmentSeedGenerator.h.

m_signComboIndex

std::size_t MuonR4::MdtSegmentSeedGenerator::m_signComboIndex {0}

private

Index of the left-right ambiguity between the circles.

Definition at line 219 of file MdtSegmentSeedGenerator.h.

m_upperHitIndex

std::size_t MuonR4::MdtSegmentSeedGenerator::m_upperHitIndex {0}

private

Explicit hit to pick in the selected top layer.

Definition at line 217 of file MdtSegmentSeedGenerator.h.

m_upperLayer

std::size_t MuonR4::MdtSegmentSeedGenerator::m_upperLayer {0}

private

Considered layer to pick the top drift circle from.

Definition at line 211 of file MdtSegmentSeedGenerator.h.

s_signCombos

constexpr static std::array<SignComboType,4> MuonR4::MdtSegmentSeedGenerator::s_signCombos

staticconstexprprivate

Initial value:

{
                std::array{ 1, 1}, std::array{ 1,-1}, 
                std::array{-1,-1}, std::array{-1, 1},  
            }

Definition at line 102 of file MdtSegmentSeedGenerator.h.

---

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

• MdtSegmentSeedGenerator.h
• MdtSegmentSeedGenerator.cxx