ATLAS Offline Software
Loading...
Searching...
No Matches
Muon::MSVertexTrackletTool Class Reference

#include <MSVertexTrackletTool.h>

Inheritance diagram for Muon::MSVertexTrackletTool:
Collaboration diagram for Muon::MSVertexTrackletTool:

Public Member Functions

 MSVertexTrackletTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~MSVertexTrackletTool ()=default
virtual StatusCode initialize () override
StatusCode findTracklets (std::vector< Tracklet > &tracklets, const EventContext &ctx) const override
 DeclareInterfaceID (Muon::IMSVertexTrackletTool, 1, 0)
 access to tool interface
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysInitialize () override
 Perform system initialization for an algorithm.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Private Types

typedef ServiceHandle< StoreGateSvc > StoreGateSvc_t

Private Member Functions

bool IgnoreMDTChamber (const Muon::MdtPrepData *mdtHit) const
int SortMDThits (std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt, const EventContext &ctx) const
void addMDTHits (std::vector< const Muon::MdtPrepData * > &hits, std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt) const
std::vector< TrackletSegment > TrackletSegmentFitter (const std::vector< const Muon::MdtPrepData * > &mdts) const
std::vector< std::pair< double, double > > SegSeeds (const std::vector< const Muon::MdtPrepData * > &mdts) const
std::vector< TrackletSegment > TrackletSegmentFitterCore (const std::vector< const Muon::MdtPrepData * > &mdts, const std::vector< std::pair< double, double > > &SeedParams) const
std::vector< TrackletSegment > CleanSegments (const std::vector< TrackletSegment > &segs) const
bool DeltabCalc (const TrackletSegment &ML1seg, const TrackletSegment &ML2seg) const
double TrackMomentum (const Identifier trkID, const double deltaAlpha) const
double TrackMomentumError (const TrackletSegment &ml1, const TrackletSegment &ml2) const
double TrackMomentumError (const TrackletSegment &ml1) const
std::vector< Tracklet > ResolveAmbiguousTracklets (std::vector< Tracklet > &tracks) const
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Static Private Member Functions

static double SeedResiduals (const std::vector< const Muon::MdtPrepData * > &mdts, double slope, double inter)
static void convertToTrackParticles (std::vector< Tracklet > &tracklets, SG::WriteHandle< xAOD::TrackParticleContainer > &container)

Private Attributes

SG::ReadHandleKey< Muon::MdtPrepDataContainer > m_mdtTESKey {this, "mdtTES", "MDT_DriftCircles"}
SG::WriteHandleKey< xAOD::TrackParticleContainer > m_TPContainer {this, "xAODTrackParticleContainer", "MSonlyTracklets"}
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}
Gaudi::Property< double > m_d12_max {this, "d12_max", 50.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 2"}
Gaudi::Property< double > m_d13_max {this, "d13_max", 80.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 3"}
Gaudi::Property< double > m_errorCutOff {this, "errorCutOff", 0.001, "minimal hit error"}
Gaudi::Property< double > m_SeedResidual {this, "SeedResidual", 5., "max residual for tracklet seeds"}
Gaudi::Property< double > m_minSegFinderChi2 {this, "MinSegFinderChi2Prob", 0.05, "tracklet segment fitter chi^2 probability cut"}
Gaudi::Property< double > m_BarrelDeltaAlphaCut {this, "BarrelDeltaAlphaCut", 0.2*Gaudi::Units::radian, "maximum delta_alpha allowed in barrel MS chambers"}
Gaudi::Property< double > m_EndcapDeltaAlphaCut {this, "EndcapDeltaAlphaCut", 0.015*Gaudi::Units::radian, "maximum delta_alpha allowed in the endcap MS chambers"}
Gaudi::Property< double > m_maxDeltabCut {this, "maxDeltabCut", 3*Gaudi::Units::millimeter, "maximum delta_b allowed"}
Gaudi::Property< double > m_minpTot {this, "minpTot", 800.*Gaudi::Units::MeV, "minimum measurable total momentum in MeV"}
Gaudi::Property< double > m_maxpTot {this, "maxpTot", 10000.*Gaudi::Units::MeV, "maximum measurable total momentum in MeV beyond which tracklets are assumed to be straight"}
Gaudi::Property< double > m_straightTrackletpTot {this, "straightTrackletpTot", 1.0e5*Gaudi::Units::MeV, "total momentum in MeV assigned to straight tracklets"}
Gaudi::Property< double > m_straightTrackletInvPerr {this, "straightTrackletInvPerr", 5.0e-5/Gaudi::Units::MeV, "error in the inverse momentum in MeV^-1 assigned to straight tracklets"}
Gaudi::Property< bool > m_tightTrackletRequirement {this, "tightTrackletRequirement", false, "tight tracklet requirement (affects efficiency - disabled by default)"}
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default).
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default).
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared

Detailed Description

Definition at line 35 of file MSVertexTrackletTool.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t
privateinherited

Definition at line 376 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ MSVertexTrackletTool()

Muon::MSVertexTrackletTool::MSVertexTrackletTool ( const std::string & type,
const std::string & name,
const IInterface * parent )

Definition at line 27 of file MSVertexTrackletTool.cxx.

27 :
28 AthAlgTool(type, name, parent) {
29 declareInterface<IMSVertexTrackletTool>(this);
30 }
AthAlgTool()
Default constructor:

◆ ~MSVertexTrackletTool()

virtual Muon::MSVertexTrackletTool::~MSVertexTrackletTool ( )
virtualdefault

Member Function Documentation

◆ addMDTHits()

void Muon::MSVertexTrackletTool::addMDTHits ( std::vector< const Muon::MdtPrepData * > & hits,
std::vector< std::vector< const Muon::MdtPrepData * > > & SortedMdt ) const
private

Definition at line 432 of file MSVertexTrackletTool.cxx.

433 {
434 if (hits.empty()) return;
435
436 // calculate number of hits in ML
437 int ntubes = hits.front()->detectorElement()->getNLayers() * hits.front()->detectorElement()->getNtubesperlayer();
438 if (hits.size() > 0.75 * ntubes) return;
439 std::sort(hits.begin(), hits.end(), [this](const Muon::MdtPrepData* mprd1, const Muon::MdtPrepData* mprd2) -> bool {
440 if (m_idHelperSvc->mdtIdHelper().tubeLayer(mprd1->identify()) > m_idHelperSvc->mdtIdHelper().tubeLayer(mprd2->identify()))
441 return false;
442 if (m_idHelperSvc->mdtIdHelper().tubeLayer(mprd1->identify()) < m_idHelperSvc->mdtIdHelper().tubeLayer(mprd2->identify()))
443 return true;
444 if (m_idHelperSvc->mdtIdHelper().tube(mprd1->identify()) < m_idHelperSvc->mdtIdHelper().tube(mprd2->identify())) return true;
445 return false;
446 }); // sort the MDTs by layer and tube number
447
448 SortedMdt.push_back(hits);
449 }
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.

◆ CleanSegments()

std::vector< TrackletSegment > Muon::MSVertexTrackletTool::CleanSegments ( const std::vector< TrackletSegment > & segs) const
private

Definition at line 742 of file MSVertexTrackletTool.cxx.

742 {
743 std::vector<TrackletSegment> CleanSegs;
744 std::vector<TrackletSegment> segs = Segs; // set of segments to perform cleaning on
745 bool keepCleaning(true);
746 int nItr(0);
747
748 while (keepCleaning) {
749 ++nItr;
750 keepCleaning = false;
751
752 for (std::vector<TrackletSegment>::iterator it = segs.begin(); it != segs.end(); ++it) {
753 if (it->isCombined()) continue;
754 std::vector<TrackletSegment> segsToCombine;
755 double tanTh1 = std::tan(it->alpha());
756 double r1 = it->globalPosition().perp();
757 double zi1 = it->globalPosition().z() - r1 / tanTh1;
758 // find all segments with similar parameters & attempt to combine
759 for (std::vector<TrackletSegment>::iterator sit = (it + 1); sit != segs.end(); ++sit) {
760 if (sit->isCombined()) continue;
761 if (it->mdtChamber() != sit->mdtChamber()) continue; // require the segments are in the same chamber
762 if ((it->mdtChEta()) * (sit->mdtChEta()) < 0) continue; // check both segments are on the same side of the detector
763 if (it->mdtChPhi() != sit->mdtChPhi()) continue; // in the same sector
764 if (std::abs(it->alpha() - sit->alpha()) > 0.005) continue; // same trajectory
765 double tanTh2 = std::tan(sit->alpha());
766 double r2 = sit->globalPosition().perp();
767 double zi2 = sit->globalPosition().z() - r2 / tanTh2;
768 // find the distance at the midpoint between the two segments
769 double rmid = (r1 + r2) / 2.;
770 double z1 = rmid / tanTh1 + zi1;
771 double z2 = rmid / tanTh2 + zi2;
772 double zdist = std::abs(z1 - z2);
773 if (zdist < 0.5) {
774 segsToCombine.push_back(*sit);
775 sit->isCombined(true);
776 }
777 } // end sit loop
778
779 // if the segment is unique, keep it
780 if (segsToCombine.empty()) {
781 CleanSegs.push_back(*it);
782 }
783 // else, combine all like segments & refit
784 else if (!segsToCombine.empty()) {
785 // create a vector of all unique MDT hits in the segments
786 std::vector<const Muon::MdtPrepData*> mdts = it->mdtHitsOnTrack();
787 for (const TrackletSegment &seg : segsToCombine) {
788 std::vector<const Muon::MdtPrepData*> tmpmdts = seg.mdtHitsOnTrack();
789 for (const Muon::MdtPrepData *tmpprd : tmpmdts){
790 bool isNewHit(true);
791 for (const Muon::MdtPrepData *tmpprd2 : mdts){
792 if (tmpprd->identify() == tmpprd2->identify()) {
793 isNewHit = false;
794 break;
795 }
796 }
797 if (isNewHit && Amg::error(tmpprd->localCovariance(), Trk::locR) > m_errorCutOff) mdts.push_back(tmpprd);
798 }
799 } // end segsToCombine loop
800
801 // only need to combine if there are extra hits added to the first segment
802 if (mdts.size() > it->mdtHitsOnTrack().size()) {
803 std::vector<TrackletSegment> refitsegs = TrackletSegmentFitter(mdts);
804 // if the refit fails, what to do?
805 if (refitsegs.empty()) {
806 if (segsToCombine.size() == 1) {
807 segsToCombine[0].isCombined(false);
808 CleanSegs.push_back(*it);
809 CleanSegs.push_back(segsToCombine[0]);
810 } else {
811 // loop on the mdts and count the number of segments that share that hit
812 std::vector<int> nSeg;
813 for (unsigned int i = 0; i < mdts.size(); ++i) {
814 nSeg.push_back(0);
815 // hit belongs to the first segment
816 for (unsigned int k = 0; k < it->mdtHitsOnTrack().size(); ++k) {
817 if (it->mdtHitsOnTrack()[k]->identify() == mdts[i]->identify()) {
818 ++nSeg[i];
819 break;
820 }
821 }
822 // hit belongs to one of the duplicate segments
823 for (unsigned int k = 0; k < segsToCombine.size(); ++k) {
824 for (unsigned int m = 0; m < segsToCombine[k].mdtHitsOnTrack().size(); ++m) {
825 if (segsToCombine[k].mdtHitsOnTrack()[m]->identify() == mdts[i]->identify()) {
826 ++nSeg[i];
827 break;
828 }
829 } // end loop on mdtHitsOnTrack
830 } // end loop on segsToCombine
831 } // end loop on mdts
832
833 // loop over the duplicates and remove the MDT used by the fewest segments until the fit converges
834 bool keeprefitting(true);
835 int nItr2(0);
836 while (keeprefitting) {
837 ++nItr2;
838 int nMinSeg(nSeg[0]);
839 const Muon::MdtPrepData* minmdt = mdts[0];
840 std::vector<int> nrfsegs;
841 std::vector<const Muon::MdtPrepData*> refitmdts;
842 // loop on MDTs, identify the overlapping set of hits
843 for (unsigned int i = 1; i < mdts.size(); ++i) {
844 if (nSeg[i] < nMinSeg) {
845 refitmdts.push_back(minmdt);
846 nrfsegs.push_back(nMinSeg);
847 minmdt = mdts[i];
848 nMinSeg = nSeg[i];
849 } else {
850 refitmdts.push_back(mdts[i]);
851 nrfsegs.push_back(nSeg[i]);
852 }
853 }
854 // reset the list of MDTs & the minimum number of segments an MDT must belong to
855 mdts = refitmdts;
856 nSeg = nrfsegs;
857 // try to fit the new set of MDTs
858 refitsegs = TrackletSegmentFitter(mdts);
859 if (!refitsegs.empty()) {
860 for (const TrackletSegment &refitseg : refitsegs) CleanSegs.push_back(refitseg);
861 keeprefitting = false; // stop refitting if segments are found
862 } else if (mdts.size() <= 3) {
863 CleanSegs.push_back(*it);
864 keeprefitting = false;
865 }
866 if (nItr2 > 10) break;
867 } // end while
868 }
869 } else {
870 keepCleaning = true;
871 for (const TrackletSegment &refitseg : refitsegs) CleanSegs.push_back(refitseg);
872 }
873 }
874 // if there are no extra MDT hits, keep only the first segment as unique
875 else
876 CleanSegs.push_back(*it);
877 }
878 } // end it loop
879 if (keepCleaning) {
880 segs = CleanSegs;
881 CleanSegs.clear();
882 }
883 if (nItr > 10) break;
884 } // end while
885
886 return CleanSegs;
887 }
size_t size() const
Number of registered mappings.
std::vector< TrackletSegment > TrackletSegmentFitter(const std::vector< const Muon::MdtPrepData * > &mdts) const
Gaudi::Property< double > m_errorCutOff
double error(const Amg::MatrixX &mat, int index)
return diagonal error of the matrix caller should ensure the matrix is symmetric and the index is in ...
@ locR
Definition ParamDefs.h:44
const Identifier & identify(const UncalibratedMeasurement *meas)
Returns the associated identifier from the muon measurement.

◆ convertToTrackParticles()

void Muon::MSVertexTrackletTool::convertToTrackParticles ( std::vector< Tracklet > & tracklets,
SG::WriteHandle< xAOD::TrackParticleContainer > & container )
staticprivate

Definition at line 330 of file MSVertexTrackletTool.cxx.

331 {
332 // convert tracklets to xAOD::TrackParticle and store in a TrackCollection
333 for (Tracklet &tracklet : tracklets) {
334 xAOD::TrackParticle* trackparticle = new xAOD::TrackParticle();
335 tracklet.setTrackParticle(trackparticle);
336 container->push_back(trackparticle);
337
338 AmgSymMatrix(5) covariance{tracklet.errorMatrix()};
339 auto MyPerigee(std::make_unique<Trk::Perigee>(tracklet.globalPosition(), tracklet.momentum(), tracklet.charge(), Trk::PerigeeSurface(Amg::Vector3D::Zero()), covariance));
340
341 // fill the xAOD::TrackParticle with the tracklet content
342 trackparticle->setDefiningParameters(MyPerigee->parameters()[Trk::d0], MyPerigee->parameters()[Trk::z0],
343 MyPerigee->parameters()[Trk::phi0], MyPerigee->parameters()[Trk::theta],
344 MyPerigee->parameters()[Trk::qOverP]);
345 trackparticle->setFitQuality(1., (float)tracklet.mdtHitsOnTrack().size());
347 std::vector<float> covMatrixVec;
348 Amg::compress(covariance, covMatrixVec);
349 trackparticle->setDefiningParametersCovMatrixVec(covMatrixVec);
350 }
351 return;
352 }
#define AmgSymMatrix(dim)
void setFitQuality(float chiSquared, float numberDoF)
Set the 'Fit Quality' information.
void setDefiningParameters(float d0, float z0, float phi0, float theta, float qOverP)
Set the defining parameters.
void setTrackProperties(const TrackProperties properties)
Methods setting the TrackProperties.
void setDefiningParametersCovMatrixVec(const std::vector< float > &cov)
void compress(const AmgSymMatrix(N) &covMatrix, std::vector< float > &vec)
const SG::AuxVectorData * container() const
Return the container holding this element.
@ phi0
Definition ParamDefs.h:65
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64
TrackParticle_v1 TrackParticle
Reference the current persistent version:
@ LowPtTrack
A LowPt track.

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)

◆ DeclareInterfaceID()

Muon::IMSVertexTrackletTool::DeclareInterfaceID ( Muon::IMSVertexTrackletTool ,
1 ,
0  )
inherited

access to tool interface

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ DeltabCalc()

bool Muon::MSVertexTrackletTool::DeltabCalc ( const TrackletSegment & ML1seg,
const TrackletSegment & ML2seg ) const
private

Definition at line 891 of file MSVertexTrackletTool.cxx.

891 {
892 double ChMid = (ML1seg.getChMidPoint() + ML2seg.getChMidPoint()) / 2.0;
893 // Calculate the Delta b (see http://inspirehep.net/record/1266438)
894 double mid1(100), mid2(1000);
895 double deltab(100);
896 if (m_idHelperSvc->mdtIdHelper().isBarrel(ML1seg.getIdentifier())) {
897 // delta b in the barrel
898 mid1 = (ChMid - ML1seg.globalPosition().perp()) / std::tan(ML1seg.alpha()) + ML1seg.globalPosition().z();
899 mid2 = (ChMid - ML2seg.globalPosition().perp()) / std::tan(ML2seg.alpha()) + ML2seg.globalPosition().z();
900 double r01 = ML1seg.globalPosition().perp() - ML1seg.globalPosition().z() * std::tan(ML1seg.alpha());
901 double r02 = ML2seg.globalPosition().perp() - ML2seg.globalPosition().z() * std::tan(ML2seg.alpha());
902 deltab = (mid2 * std::tan(ML1seg.alpha()) - ChMid + r01) / (std::hypot(1,std::tan(ML1seg.alpha())));
903 double deltab2 = (mid1 * std::tan(ML2seg.alpha()) - ChMid + r02) / (std::hypot(1,std::tan(ML2seg.alpha())));
904 if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
905 } else {
906 // delta b in the endcap
907 mid1 = ML1seg.globalPosition().perp() + std::tan(ML1seg.alpha()) * (ChMid - ML1seg.globalPosition().z());
908 mid2 = ML2seg.globalPosition().perp() + std::tan(ML2seg.alpha()) * (ChMid - ML2seg.globalPosition().z());
909 double z01 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
910 double z02 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
911 deltab = (mid2 / std::tan(ML1seg.alpha()) - ChMid + z01) / (std::hypot(1,1/std::tan(ML1seg.alpha())));
912 double deltab2 = (mid1 / std::tan(ML2seg.alpha()) - ChMid + z02) / (std::hypot(1,1/std::tan(ML2seg.alpha())));
913 if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
914 }
915
916 // calculate the maximum allowed Delta b based on delta alpha uncertainties and ML spacing
917 double dbmax = 5 * std::abs(ChMid - ML1seg.getChMidPoint()) * std::hypot(ML1seg.alphaError(), ML2seg.alphaError());
918 if (dbmax > m_maxDeltabCut) dbmax = m_maxDeltabCut;
919 return std::abs(deltab) < dbmax;
920 }
Scalar perp() const
perp method - perpendicular length
#define z
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Gaudi::Property< double > m_maxDeltabCut
double alpha() const
const Identifier getIdentifier() const
const Amg::Vector3D & globalPosition() const
double alphaError() const
double getChMidPoint() const

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ findTracklets()

StatusCode Muon::MSVertexTrackletTool::findTracklets ( std::vector< Tracklet > & tracklets,
const EventContext & ctx ) const
overridevirtual

Implements Muon::IMSVertexTrackletTool.

Definition at line 44 of file MSVertexTrackletTool.cxx.

44 {
45 // record TrackParticle container in StoreGate
46 SG::WriteHandle<xAOD::TrackParticleContainer> container(m_TPContainer, ctx);
47 ATH_CHECK(container.record(std::make_unique<xAOD::TrackParticleContainer>(), std::make_unique<xAOD::TrackParticleAuxContainer>()));
48
49 // sort the MDT hits into chambers & MLs
50 std::vector<std::vector<const Muon::MdtPrepData*> > SortedMdt;
51
52 int nMDT = SortMDThits(SortedMdt, ctx);
53
54 if (nMDT <= 0) { return StatusCode::SUCCESS; }
55
56 if (msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG("MDT hits are selected and sorted");
57
58 // loop over the MDT hits and find segments
59 // select the tube combinations to be fit
60 /*Select hits in at least 2 layers and require hits be ordered by increasing tube number (see diagrams below).
61 ( )( )(3)( ) ( )(3)( )( ) ( )( )( )( ) ( )( )(3)( ) ( )(2)(3)( )
62 ( )(2)( )( ) ( )(2)( )( ) ( )(2)(3)( ) ( )(1)(2)( ) ( )(1)( )( )
63 (1)( )( )( ) (1)( )( )( ) (1)( )( )( ) ( )( )( )( ) ( )( )( )( )
64 Barrel selection criteria: |z_mdt1 - z_mdt2| < m_d12_max (50 mm), |z_mdt1 - z_mdt3| < m_d13_max (80 mm)
65 Endcap selection criteria: |r_mdt1 - r_mdt2| < m_d12_max (50 mm), |r_mdt1 - r_mdt3| < m_d13_max (80 mm)
66 */
67
68 std::vector<TrackletSegment> segs[6][2][16]; // single ML segment array (indicies [station type][ML][sector]) with station type iterating through barrel inner, middle, outer and then endcap inner, middle, outer
69 std::vector<std::vector<const Muon::MdtPrepData*> >::const_iterator ChamberItr = SortedMdt.begin();
70 for (; ChamberItr != SortedMdt.end(); ++ChamberItr) {
71 std::vector<TrackletSegment> mlsegments;
72 std::vector<const Muon::MdtPrepData*>::const_iterator mdt1 = ChamberItr->begin();
73 std::vector<const Muon::MdtPrepData*>::const_iterator mdtEnd = ChamberItr->end();
74 if (IgnoreMDTChamber(*mdt1)) continue;
75
76 // get information about current chamber
77 Identifier mdt1_ID = (*mdt1)->identify();
78 bool mdt1_isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(mdt1_ID);
79 bool mdt1_isEndcap = m_idHelperSvc->mdtIdHelper().isEndcap(mdt1_ID);
80 int sector = m_idHelperSvc->sector(mdt1_ID);
81 int maxLayer = m_idHelperSvc->mdtIdHelper().tubeLayerMax(mdt1_ID);
82 int ML = m_idHelperSvc->mdtIdHelper().multilayer(mdt1_ID);
83
84 // loop on hits inside the chamber
85 for (; mdt1 != mdtEnd; ++mdt1) {
86 if (Amg::error((*mdt1)->localCovariance(), Trk::locR) < m_errorCutOff) {
87 ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt1_ID) << " with too small error "
88 << Amg::error((*mdt1)->localCovariance(), Trk::locR));
89 continue;
90 }
91
92 int tl1 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt1_ID);
93 if (tl1 == maxLayer) break; // require hits in at least 2 layers
94
95 // loop on second hits
96 std::vector<const Muon::MdtPrepData*>::const_iterator mdt2 = (mdt1 + 1);
97 if (mdt2 == mdtEnd) continue;
98 Identifier mdt2_ID = (*mdt2)->identify();
99 for (; mdt2 != mdtEnd; ++mdt2) {
100 if (Amg::error((*mdt2)->localCovariance(), Trk::locR) < m_errorCutOff) {
101 ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt2_ID)
102 << " with too small error " << Amg::error((*mdt2)->localCovariance(), Trk::locR));
103 continue;
104 }
105
106 // reject the bad tube combinations
107 int tl2 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt2_ID);
108 if (mdt1 == mdt2 || (tl2 - tl1) > 1 || (tl2 - tl1) < 0) continue;
109 if ((tl2 - tl1) == 0 && (m_idHelperSvc->mdtIdHelper().tube(mdt2_ID) -
110 m_idHelperSvc->mdtIdHelper().tube(mdt1_ID)) < 0) continue;
111 // reject bad hit separations
112 if (mdt1_isBarrel && std::abs((*mdt1)->globalPosition().z() - (*mdt2)->globalPosition().z()) > m_d12_max) continue;
113 if (mdt1_isEndcap && std::abs((*mdt1)->globalPosition().perp() - (*mdt2)->globalPosition().perp()) > m_d12_max) continue;
114
115 // loop on third hits
116 std::vector<const Muon::MdtPrepData*>::const_iterator mdt3 = (mdt2 + 1);
117 if (mdt3 == mdtEnd) continue;
118 Identifier mdt3_ID = (*mdt3)->identify();
119 for (; mdt3 != mdtEnd; ++mdt3) {
120 if (Amg::error((*mdt3)->localCovariance(), Trk::locR) < m_errorCutOff) {
121 ATH_MSG_WARNING(" " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt3_ID)
122 << " with too small error " << Amg::error((*mdt3)->localCovariance(), Trk::locR));
123 continue;
124 }
125
126 // reject the bad tube combinations
127 if (mdt1 == mdt3 || mdt2 == mdt3) continue;
128 int tl3 = m_idHelperSvc->mdtIdHelper().tubeLayer(mdt3_ID);
129 if ((tl3 - tl2) > 1 || (tl3 - tl2) < 0 || (tl3 - tl1) <= 0) continue;
130 if ((tl3 - tl2) == 0 && (m_idHelperSvc->mdtIdHelper().tube(mdt3_ID) -
131 m_idHelperSvc->mdtIdHelper().tube(mdt2_ID)) < 0) continue;
132 // reject bad hit separations
133 if (mdt1_isBarrel && std::abs((*mdt1)->globalPosition().z() - (*mdt3)->globalPosition().z()) > m_d13_max) continue;
134 if (mdt1_isEndcap && std::abs((*mdt1)->globalPosition().perp() - (*mdt3)->globalPosition().perp()) > m_d13_max) continue;
135
136 // store and fit the good combinations
137 std::vector<const Muon::MdtPrepData*> mdts;
138 mdts.push_back((*mdt1));
139 mdts.push_back((*mdt2));
140 mdts.push_back((*mdt3));
141 std::vector<TrackletSegment> tmpSegs = TrackletSegmentFitter(mdts);
142 for (const TrackletSegment &tmpSeg : tmpSegs) mlsegments.push_back(tmpSeg);
143 } // end loop on mdt3
144 } // end loop on mdt2
145 } // end loop on mdt1
146
147 // store the reconstructed segments according to station, ML and sector
148 // MS region decoded in MuonIdHelpers/MuonIdHelper.h
149 int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(mdt1_ID);
150 if (mdt1_isBarrel){
151 if (stationRegion == 0)
152 for (const TrackletSegment &mlsegment : mlsegments) segs[0][ML - 1][sector - 1].push_back(mlsegment); // barrel inner
153 else if (stationRegion == 2)
154 for (const TrackletSegment &mlsegment : mlsegments) segs[1][ML - 1][sector - 1].push_back(mlsegment); // barrel middle
155 else if (stationRegion == 3)
156 for (const TrackletSegment &mlsegment : mlsegments) segs[2][ML - 1][sector - 1].push_back(mlsegment); // barrel outer
157 }
158 else if (mdt1_isEndcap){
159 if (stationRegion == 0)
160 for (const TrackletSegment &mlsegment : mlsegments) segs[3][ML - 1][sector - 1].push_back(mlsegment); // endcap inner
161 else if (stationRegion == 2)
162 for (const TrackletSegment &mlsegment : mlsegments) segs[4][ML - 1][sector - 1].push_back(mlsegment); // endcap middle
163 else if (stationRegion == 3)
164 for (const TrackletSegment &mlsegment : mlsegments) segs[5][ML - 1][sector - 1].push_back(mlsegment); // endcap outer
165 }
166 else
167 ATH_MSG_WARNING("Found segments belonging to chamber " << m_idHelperSvc->mdtIdHelper().stationNameString(m_idHelperSvc->mdtIdHelper().stationName(mdt1_ID)) << " that have not been stored");
168 } // end loop on mdt chambers
169
170 // Combine/remove duplicate segments
171 std::vector<TrackletSegment> CleanSegs[6][2][16];
172 for (int st = 0; st < 6; ++st) {
173 for (int ml = 0; ml < 2; ++ml) {
174 for (int sector = 0; sector < 16; ++sector) {
175 if (!segs[st][ml][sector].empty()) {
176 CleanSegs[st][ml][sector] = CleanSegments(segs[st][ml][sector]);
177 }
178 }
179 }
180 }
181
182 // loop over TrackletSegments in barrel inner, middle, outer and endcap inner, middle, outer stations
183 for (int st = 0; st < 6; ++st) {
184 double DeltaAlphaCut = m_BarrelDeltaAlphaCut;
185 for (int sector = 0; sector < 16; ++sector) {
186 for (const TrackletSegment &ML1seg : CleanSegs[st][0][sector]) {
187 // Set the delta alpha cut depending on station type
188 const Identifier trkID = ML1seg.getIdentifier();
189 bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
190 bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
191 int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
192
193 if (isBarrel){
194 if (stationRegion == 0){
195 if (isSmall) DeltaAlphaCut = m_BarrelDeltaAlphaCut; // default value for BIS
196 else DeltaAlphaCut = c_BIL / 750.0;
197 }
198 else if (stationRegion == 2){
199 if (isSmall) DeltaAlphaCut = c_BMS / 750.0;
200 else DeltaAlphaCut = c_BML / 750.0;
201 }
202 else if (stationRegion == 3){
203 if (isSmall) DeltaAlphaCut = m_BarrelDeltaAlphaCut; // default value for BOS
204 else DeltaAlphaCut = c_BOL / 750.0;
205 }
206 }
207 else{
208 DeltaAlphaCut = m_EndcapDeltaAlphaCut;
209 }
210
211 // loop on ML2 segments from same sector
212 for (const TrackletSegment &ML2seg : CleanSegs[st][1][sector]) {
213 if (ML1seg.mdtChamber() != ML2seg.mdtChamber() || ML1seg.mdtChEta() != ML2seg.mdtChEta()) continue;
214
215 double deltaAlpha = ML1seg.alpha() - ML2seg.alpha();
216 bool goodDeltab = DeltabCalc(ML1seg, ML2seg);
217 // select the good combinations
218 if (std::abs(deltaAlpha) < DeltaAlphaCut && goodDeltab) {
219 if (isBarrel) {
220 // barrel chambers
221 double charge_discriminant = deltaAlpha * ML1seg.globalPosition().z() * std::tan(ML1seg.alpha());
222 double charge = charge_discriminant < 0 ? -1 : 1;
223
224 double pTot = TrackMomentum(ML1seg.getIdentifier(), deltaAlpha);
225 if (pTot < m_minpTot) continue;
226 if (pTot > m_maxpTot) {
227 // if we find a straight track, try to do a global refit to minimize the number of duplicates
228 charge = 0;
229 std::vector<const Muon::MdtPrepData*> mdts = ML1seg.mdtHitsOnTrack();
230 std::vector<const Muon::MdtPrepData*> mdts2 = ML2seg.mdtHitsOnTrack();
231 for (const Muon::MdtPrepData *mdt2 : mdts2) mdts.push_back(mdt2);
232 std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
233
234 if (!CombinedSeg.empty()) {
235 // calculate momentum components & uncertainty
236 double Trk1overPErr = TrackMomentumError(CombinedSeg[0]);
237 double pT = pTot * std::sin(CombinedSeg[0].alpha());
238 double pz = pTot * std::cos(CombinedSeg[0].alpha());
239 Amg::Vector3D momentum(pT * std::cos(CombinedSeg[0].globalPosition().phi()),
240 pT * std::sin(CombinedSeg[0].globalPosition().phi()),
241 pz);
242 // create the error matrix
243 AmgSymMatrix(5) matrix;
244 matrix.setIdentity();
245 matrix(0, 0) = std::pow(CombinedSeg[0].rError(),2); // delta locR
246 matrix(1, 1) = std::pow(CombinedSeg[0].zError(),2); // delta locz
247 matrix(2, 2) = std::pow(0.00000000001,2); // delta phi (~0 because we explicitly rotate all tracklets into
248 // the middle of the chamber)
249 matrix(3, 3) = std::pow(CombinedSeg[0].alphaError(),2); // delta theta
250 matrix(4, 4) = std::pow(Trk1overPErr,2); // delta 1/p
251 Tracklet tmpTrk(CombinedSeg[0], momentum, matrix, charge);
252 ATH_MSG_DEBUG("Track " << tracklets.size() << " found with p = (" << momentum.x() << ", "
253 << momentum.y() << ", " << momentum.z()
254 << ") and |p| = " << tmpTrk.momentum().mag() << " MeV");
255 tracklets.push_back(tmpTrk);
256 }
257 } else {
258 // tracklet has a measurable momentum
259 double Trk1overPErr = TrackMomentumError(ML1seg, ML2seg);
260 double pT = pTot * std::sin(ML1seg.alpha());
261 double pz = pTot * std::cos(ML1seg.alpha());
262 Amg::Vector3D momentum(pT * std::cos(ML1seg.globalPosition().phi()),
263 pT * std::sin(ML1seg.globalPosition().phi()),
264 pz);
265 // create the error matrix
266 AmgSymMatrix(5) matrix;
267 matrix.setIdentity();
268 matrix(0, 0) = std::pow(ML1seg.rError(),2); // delta locR
269 matrix(1, 1) = std::pow(ML1seg.zError(),2); // delta locz
270 matrix(2, 2) = std::pow(0.00000000001,2); // delta phi (~0 because we explicitly rotate all tracks into the
271 // middle of the chamber)
272 matrix(3, 3) = std::pow(ML1seg.alphaError(),2); // delta theta
273 matrix(4, 4) = std::pow(Trk1overPErr,2); // delta 1/p
274 Tracklet tmpTrk(ML1seg, ML2seg, momentum, matrix, charge);
275 ATH_MSG_DEBUG("Track " << tracklets.size() << " found with p = (" << momentum.x() << ", "
276 << momentum.y() << ", " << momentum.z()
277 << ") and |p| = " << tmpTrk.momentum().mag() << " MeV");
278 tracklets.push_back(tmpTrk);
279 }
280 } // end barrel chamber selection
281 else if (!isBarrel) {
282 // endcap tracklets
283 // always straight tracklets (no momentum measurement possible)
284 std::vector<const Muon::MdtPrepData*> mdts = ML1seg.mdtHitsOnTrack();
285 std::vector<const Muon::MdtPrepData*> mdts2 = ML2seg.mdtHitsOnTrack();
286 for (const Muon::MdtPrepData *mdt2 : mdts2) mdts.push_back(mdt2);
287 std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
288
289 if (!CombinedSeg.empty()) {
290 double charge = 0;
291 double pTot = m_straightTrackletpTot;
292 double pT = pTot * std::sin(CombinedSeg[0].alpha());
293 double pz = pTot * std::cos(CombinedSeg[0].alpha());
294 Amg::Vector3D momentum(pT * std::cos(CombinedSeg[0].globalPosition().phi()),
295 pT * std::sin(CombinedSeg[0].globalPosition().phi()),
296 pz);
297 // create the error matrix
298 AmgSymMatrix(5) matrix;
299 matrix.setIdentity();
300 matrix(0, 0) = std::pow(CombinedSeg[0].rError(),2); // delta locR
301 matrix(1, 1) = std::pow(CombinedSeg[0].zError(),2); // delta locz
302 matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle
303 // of the chamber)
304 matrix(3, 3) = std::pow(CombinedSeg[0].alphaError(),2); // delta theta
305 matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
306
307 Tracklet tmpTrk(CombinedSeg[0], momentum, matrix, charge);
308 tracklets.push_back(tmpTrk);
309 }
310 } // end endcap tracklet selection
311
312 } // end tracklet selection (delta alpha & delta b)
313
314 } // end loop on ML2 segments
315 } // end loop on ML1 segments
316 } // end loop on sectors
317 } // end loop on stations
318
319 // Resolve any ambiguous tracklets
320 tracklets = ResolveAmbiguousTracklets(tracklets);
321
322 // convert from tracklets to Trk::Tracks
323 convertToTrackParticles(tracklets, container);
324
325 return StatusCode::SUCCESS;
326 }
Scalar phi() const
phi method
Scalar mag() const
mag method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_WARNING(x,...)
double charge(const T &p)
Definition AtlasPID.h:1003
if(pathvar)
#define y
#define x
static const Attributes_t empty
bool msgLvl(const MSG::Level lvl) const
int SortMDThits(std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt, const EventContext &ctx) const
SG::WriteHandleKey< xAOD::TrackParticleContainer > m_TPContainer
double TrackMomentumError(const TrackletSegment &ml1, const TrackletSegment &ml2) const
static void convertToTrackParticles(std::vector< Tracklet > &tracklets, SG::WriteHandle< xAOD::TrackParticleContainer > &container)
bool DeltabCalc(const TrackletSegment &ML1seg, const TrackletSegment &ML2seg) const
double TrackMomentum(const Identifier trkID, const double deltaAlpha) const
bool IgnoreMDTChamber(const Muon::MdtPrepData *mdtHit) const
Gaudi::Property< double > m_d13_max
std::vector< Tracklet > ResolveAmbiguousTracklets(std::vector< Tracklet > &tracks) const
Gaudi::Property< double > m_minpTot
std::vector< TrackletSegment > CleanSegments(const std::vector< TrackletSegment > &segs) const
Gaudi::Property< double > m_EndcapDeltaAlphaCut
Gaudi::Property< double > m_BarrelDeltaAlphaCut
Gaudi::Property< double > m_d12_max
Gaudi::Property< double > m_maxpTot
Eigen::Matrix< double, 3, 1 > Vector3D
bool isSmall(const ChIndex index)
Returns true if the chamber index is in a small sector.
bool isBarrel(const ChIndex index)
Returns true if the chamber index points to a barrel chamber.
constexpr double c_BML
constexpr double c_BIL
constexpr double c_BOL
constexpr double c_BMS

◆ IgnoreMDTChamber()

bool Muon::MSVertexTrackletTool::IgnoreMDTChamber ( const Muon::MdtPrepData * mdtHit) const
private

Definition at line 356 of file MSVertexTrackletTool.cxx.

356 {
357 // return true if the MDT hit is in a chamber to be ignored. These hits are then not used to reconstruct tracklets.
358
359 bool ignore = false;
360 int stName = m_idHelperSvc->mdtIdHelper().stationName(mdtHit->identify());
361 int stEta = m_idHelperSvc->mdtIdHelper().stationEta(mdtHit->identify());
362
363 // Doesn't consider hits belonging to chambers BEE, EEL and EES
364 if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BEE") ||
365 stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("EEL") ||
366 stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("EES")) ignore = true;
367
368 // Doesn't consider hits belonging to chambers BIS7/8
369 if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BIS") && std::abs(stEta) >= 7) ignore = true;
370
371 // Doesn't consider hits belonging to BME or BMG chambers
372 if (stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BME") ||
373 stName == m_idHelperSvc->mdtIdHelper().stationNameIndex("BMG")) ignore = true;
374
375 return ignore;
376 }
Identifier identify() const
return the identifier
const std::string & stName(StIndex index)
convert StIndex into a string

◆ initialize()

StatusCode Muon::MSVertexTrackletTool::initialize ( )
overridevirtual

Definition at line 34 of file MSVertexTrackletTool.cxx.

34 {
35 ATH_CHECK(m_mdtTESKey.initialize());
36 ATH_CHECK(m_TPContainer.initialize());
37 ATH_CHECK(m_idHelperSvc.retrieve());
38
39 return StatusCode::SUCCESS;
40 }
SG::ReadHandleKey< Muon::MdtPrepDataContainer > m_mdtTESKey

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ msg()

MsgStream & AthCommonMsg< AlgTool >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 368 of file AthCommonDataStore.h.

369 {
370 h.renounce();
372 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 352 of file AthCommonDataStore.h.

352 {
354 }

◆ ResolveAmbiguousTracklets()

std::vector< Tracklet > Muon::MSVertexTrackletTool::ResolveAmbiguousTracklets ( std::vector< Tracklet > & tracks) const
private

Definition at line 1012 of file MSVertexTrackletTool.cxx.

1012 {
1013 ATH_MSG_DEBUG("In ResolveAmbiguousTracks");
1014 // considering only tracklets with the number of associated hits
1015 // being more than 3/4 the number of layers in the MS chamber
1016
1018 std::vector<Tracklet> myTracks = tracks;
1019 tracks.clear();
1020 for (const Tracklet &track : myTracks) {
1021 Identifier id1 = track.getML1seg().mdtHitsOnTrack().at(0)->identify();
1022 Identifier id2 = track.getML2seg().mdtHitsOnTrack().at(0)->identify();
1023 int nLayerML1 = m_idHelperSvc->mdtIdHelper().tubeLayerMax(id1);
1024 int nLayerML2 = m_idHelperSvc->mdtIdHelper().tubeLayerMax(id2);
1025 double ratio = (double)(track.mdtHitsOnTrack().size()) / (nLayerML1 + nLayerML2);
1026 if (ratio > 0.75) tracks.push_back(track);
1027 }
1028 }
1029
1030 std::vector<Tracklet> UniqueTracks;
1031 std::vector<unsigned int> AmbigTrks; // indices of ambigious tracklets
1032 for (unsigned int tk1 = 0; tk1 < tracks.size(); ++tk1) {
1033 int nShared = 0;
1034 // check if any Ambiguity has been broken
1035 bool isResolved = false;
1036 for (unsigned int AmbigTrksIdx : AmbigTrks) {
1037 if (tk1 == AmbigTrksIdx) {
1038 isResolved = true;
1039 break;
1040 }
1041 }
1042 if (isResolved) continue;
1043 std::vector<Tracklet> AmbigTracks;
1044 AmbigTracks.push_back(tracks.at(tk1));
1045 // get a point on the track
1046 double Trk1ML1R = tracks.at(tk1).getML1seg().globalPosition().perp();
1047 double Trk1ML1Z = tracks.at(tk1).getML1seg().globalPosition().z();
1048 double Trk1ML2R = tracks.at(tk1).getML2seg().globalPosition().perp();
1049 double Trk1ML2Z = tracks.at(tk1).getML2seg().globalPosition().z();
1050
1051 Identifier tk1ID = tracks.at(tk1).muonIdentifier();
1052 bool tk1_isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(tk1ID);
1053 bool tk1_isEndcap = m_idHelperSvc->mdtIdHelper().isEndcap(tk1ID);
1054
1055 // loop over the rest of the tracks and find any amibuities
1056 for (unsigned int tk2 = (tk1 + 1); tk2 < tracks.size(); ++tk2) {
1057 if (tracks.at(tk1).mdtChamber() == tracks.at(tk2).mdtChamber() && tracks.at(tk1).mdtChPhi() == tracks.at(tk2).mdtChPhi() &&
1058 (tracks.at(tk1).mdtChEta()) * (tracks.at(tk2).mdtChEta()) > 0) {
1059 // check if any Ambiguity has been broken
1060 for (unsigned int AmbigTrksIdx : AmbigTrks) {
1061 if (tk2 == AmbigTrksIdx) {
1062 isResolved = true;
1063 break;
1064 }
1065 }
1066 if (isResolved) continue;
1067 // get a point on the track
1068 double Trk2ML1R = tracks.at(tk2).getML1seg().globalPosition().perp();
1069 double Trk2ML1Z = tracks.at(tk2).getML1seg().globalPosition().z();
1070 double Trk2ML2R = tracks.at(tk2).getML2seg().globalPosition().perp();
1071 double Trk2ML2Z = tracks.at(tk2).getML2seg().globalPosition().z();
1072
1073 // find the distance between the tracks
1074 double DistML1(1000), DistML2(1000);
1075 if (tk1_isBarrel) {
1076 DistML1 = std::abs(Trk1ML1Z - Trk2ML1Z);
1077 DistML2 = std::abs(Trk1ML2Z - Trk2ML2Z);
1078 } else if (tk1_isEndcap) {
1079 DistML1 = std::abs(Trk1ML1R - Trk2ML1R);
1080 DistML2 = std::abs(Trk1ML2R - Trk2ML2R);
1081 }
1082 if (DistML1 < 40 || DistML2 < 40) {
1083 // find how many MDTs the tracks share
1084 std::vector<const Muon::MdtPrepData*> mdts1 = tracks.at(tk1).mdtHitsOnTrack();
1085 std::vector<const Muon::MdtPrepData*> mdts2 = tracks.at(tk2).mdtHitsOnTrack();
1086 nShared = 0;
1087 for (const Muon::MdtPrepData *mdt1 : mdts1) {
1088 for (const Muon::MdtPrepData *mdt2 : mdts2) {
1089 if (mdt1->identify() == mdt2->identify()) {
1090 ++nShared;
1091 break;
1092 }
1093 }
1094 }
1095
1096 if (nShared <= 1) continue; // if the tracks share only 1 hits move to next track
1097 // store the track as ambiguous
1098 AmbigTracks.push_back(tracks.at(tk2));
1099 AmbigTrks.push_back(tk2);
1100 }
1101 } // end chamber match
1102 } // end tk2 loop
1103
1104 if (AmbigTracks.size() == 1) {
1105 UniqueTracks.push_back(tracks.at(tk1));
1106 continue;
1107 }
1108 // Deal with any ambiguities
1109 // Barrel tracks
1110 if (tk1_isBarrel) {
1111 bool hasMomentum = tracks.at(tk1).charge() != 0;
1112 double aveX(0), aveY(0), aveZ(0), aveAlpha(0);
1113 double aveP(0), nAmbigP(0), TrkCharge(tracks.at(tk1).charge());
1114 bool allSameSign(true);
1115
1116 for (const Tracklet &AmbigTrack : AmbigTracks) {
1117 if (!hasMomentum) {
1118 aveX += AmbigTrack.globalPosition().x();
1119 aveY += AmbigTrack.globalPosition().y();
1120 aveZ += AmbigTrack.globalPosition().z();
1121 aveAlpha += AmbigTrack.getML1seg().alpha();
1122 } else {
1123 // check the charge is the same
1124 if (std::abs(AmbigTrack.charge() - TrkCharge) > 0.1) allSameSign = false;
1125 // find the average momentum
1126 aveP += AmbigTrack.momentum().mag();
1127 ++nAmbigP;
1128 aveAlpha += AmbigTrack.alpha();
1129 aveX += AmbigTrack.globalPosition().x();
1130 aveY += AmbigTrack.globalPosition().y();
1131 aveZ += AmbigTrack.globalPosition().z();
1132 }
1133 } // end loop on ambiguous tracks
1134 if (!hasMomentum) {
1135 aveX = aveX / (double)AmbigTracks.size();
1136 aveY = aveY / (double)AmbigTracks.size();
1137 aveZ = aveZ / (double)AmbigTracks.size();
1138 Amg::Vector3D gpos(aveX, aveY, aveZ);
1139 aveAlpha = aveAlpha / (double)AmbigTracks.size();
1140 double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1141 double rErr = tracks.at(tk1).getML1seg().rError();
1142 double zErr = tracks.at(tk1).getML1seg().zError();
1143
1144 TrackletSegment aveSegML1(m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1145 double pT = m_maxpTot * std::sin(aveSegML1.alpha());
1146 double pz = m_maxpTot * std::cos(aveSegML1.alpha());
1147 Amg::Vector3D momentum(pT * std::cos(aveSegML1.globalPosition().phi()),
1148 pT * std::sin(aveSegML1.globalPosition().phi()),
1149 pz);
1150 AmgSymMatrix(5) matrix;
1151 matrix.setIdentity();
1152 matrix(0, 0) = std::pow(tracks.at(tk1).getML1seg().rError(),2); // delta R
1153 matrix(1, 1) = std::pow(tracks.at(tk1).getML1seg().zError(),2); // delta z
1154 matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1155 matrix(3, 3) = std::pow(tracks.at(tk1).getML1seg().alphaError(),2); // delta theta
1156 matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p
1157 Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1158 UniqueTracks.push_back(aveTrack);
1159 } else if (allSameSign) {
1160 aveP = aveP / nAmbigP;
1161 double pT = aveP * std::sin(tracks.at(tk1).getML1seg().alpha());
1162 double pz = aveP * std::cos(tracks.at(tk1).getML1seg().alpha());
1163 Amg::Vector3D momentum(pT * std::cos(tracks.at(tk1).globalPosition().phi()),
1164 pT * std::sin(tracks.at(tk1).globalPosition().phi()),
1165 pz);
1166 Tracklet MyTrack = tracks.at(tk1);
1167 MyTrack.momentum(momentum);
1168 MyTrack.charge(tracks.at(tk1).charge());
1169 UniqueTracks.push_back(MyTrack);
1170 } else {
1171 aveX = aveX / (double)AmbigTracks.size();
1172 aveY = aveY / (double)AmbigTracks.size();
1173 aveZ = aveZ / (double)AmbigTracks.size();
1174 Amg::Vector3D gpos(aveX, aveY, aveZ);
1175 aveAlpha = aveAlpha / (double)AmbigTracks.size();
1176 double alphaErr = tracks.at(tk1).getML1seg().alphaError();
1177 double rErr = tracks.at(tk1).getML1seg().rError();
1178 double zErr = tracks.at(tk1).getML1seg().zError();
1179
1180 TrackletSegment aveSegML1(m_idHelperSvc.get(), tracks.at(tk1).getML1seg().mdtHitsOnTrack(), gpos, aveAlpha, alphaErr, rErr, zErr, 0);
1181 double pT = m_maxpTot * std::sin(aveSegML1.alpha());
1182 double pz = m_maxpTot * std::cos(aveSegML1.alpha());
1183 Amg::Vector3D momentum(pT * std::cos(aveSegML1.globalPosition().phi()),
1184 pT * std::sin(aveSegML1.globalPosition().phi()),
1185 pz);
1186 AmgSymMatrix(5) matrix;
1187 matrix.setIdentity();
1188 matrix(0, 0) = std::pow(tracks.at(tk1).getML1seg().rError(),2); // delta R
1189 matrix(1, 1) = std::pow(tracks.at(tk1).getML1seg().zError(),2); // delta z
1190 matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1191 matrix(3, 3) = std::pow(tracks.at(tk1).getML1seg().alphaError(),2); // delta theta
1192 matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p
1193 Tracklet aveTrack(aveSegML1, momentum, matrix, 0);
1194 UniqueTracks.push_back(aveTrack);
1195 }
1196 } // end barrel Tracks
1197
1198 // Endcap tracks
1199 else if (tk1_isEndcap) {
1200 std::vector<const Muon::MdtPrepData*> AllMdts;
1201 for (Tracklet const &AmbigTrack : AmbigTracks) {
1202 std::vector<const Muon::MdtPrepData*> mdts = AmbigTrack.mdtHitsOnTrack();
1203 std::vector<const Muon::MdtPrepData*> tmpAllMdt = AllMdts;
1204 for (const Muon::MdtPrepData *mdt : mdts) {
1205 bool isNewHit = true;
1206 for (const Muon::MdtPrepData *tmpmdt : tmpAllMdt) {
1207 if (mdt->identify() == tmpmdt->identify()) {
1208 isNewHit = false;
1209 break;
1210 }
1211 }
1212 if (isNewHit) AllMdts.push_back(mdt);
1213 } // end loop on mdts
1214 } // end loop on ambiguous tracks
1215
1216 std::vector<TrackletSegment> MyECsegs = TrackletSegmentFitter(AllMdts);
1217 if (!MyECsegs.empty()) {
1218 TrackletSegment ECseg = MyECsegs.at(0);
1219 ECseg.clearMdt();
1220 double pT = m_maxpTot * std::sin(ECseg.alpha());
1221 double pz = m_maxpTot * std::cos(ECseg.alpha());
1222 Amg::Vector3D momentum(pT * std::cos(ECseg.globalPosition().phi()),
1223 pT * std::sin(ECseg.globalPosition().phi()),
1224 pz);
1225 AmgSymMatrix(5) matrix;
1226 matrix.setIdentity();
1227 matrix(0, 0) = std::pow(ECseg.rError(),2); // delta R
1228 matrix(1, 1) = std::pow(ECseg.zError(),2); // delta z
1229 matrix(2, 2) = std::pow(0.0000001,2); // delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
1230 matrix(3, 3) = std::pow(ECseg.alphaError(),2); // delta theta
1231 matrix(4, 4) = std::pow(m_straightTrackletInvPerr,2); // delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
1232 Tracklet MyCombTrack(MyECsegs.at(0), ECseg, momentum, matrix, 0);
1233 UniqueTracks.push_back(MyCombTrack);
1234 } else
1235 UniqueTracks.push_back(tracks.at(tk1));
1236 } // end endcap tracks
1237
1238 } // end loop on tracks -- tk1
1239
1240 return UniqueTracks;
1241 }
Gaudi::Property< bool > m_tightTrackletRequirement
Gaudi::Property< double > m_straightTrackletInvPerr
void momentum(const Amg::Vector3D &p)
Definition Tracklet.cxx:40
void charge(double charge)
Definition Tracklet.cxx:41

◆ SeedResiduals()

double Muon::MSVertexTrackletTool::SeedResiduals ( const std::vector< const Muon::MdtPrepData * > & mdts,
double slope,
double inter )
staticprivate

Definition at line 573 of file MSVertexTrackletTool.cxx.

573 {
574 // calculate the residual of the MDTs not used to create the seed
575 double resid = 0;
576 for (const Muon::MdtPrepData *mdt : mdts) {
577 double mdtR = mdt->globalPosition().perp();
578 double mdtZ = mdt->globalPosition().z();
579 double res = std::abs((mdt->localPosition()[Trk::locR] - std::abs((mdtR - inter - slope * mdtZ) / std::hypot(slope, 1))) /
580 (Amg::error(mdt->localCovariance(), Trk::locR)));
581 if (res > resid) resid = res;
582 }
583 return resid;
584 }
std::pair< std::vector< unsigned int >, bool > res

◆ SegSeeds()

std::vector< std::pair< double, double > > Muon::MSVertexTrackletTool::SegSeeds ( const std::vector< const Muon::MdtPrepData * > & mdts) const
private

Definition at line 465 of file MSVertexTrackletTool.cxx.

465 {
466 std::vector<std::pair<double, double> > SeedParams;
467 if (mdts.empty())[[unlikely]]{
468 ATH_MSG_DEBUG("SegSeeds called with an empty vector.");
469 return SeedParams;
470 }
471 // create seeds by drawing the 4 possible lines tangent to the two outermost drift circles
472 // see http://cds.cern.ch/record/620198 (section 4.3) for description of the algorithm
473 // keep all seeds which satisfy the criterion: residual(mdt 2) < m_SeedResidual
474 // NOTE: here there is an assumption that each MDT has a radius of 30mm
475 // -- needs to be revisited when the small tubes in sectors 12 & 14 are installed
476 double x1 = mdts.front()->globalPosition().z();
477 double y1 = mdts.front()->globalPosition().perp();
478 double r1 = std::abs(mdts.front()->localPosition()[Trk::locR]);
479
480 double x2 = mdts.back()->globalPosition().z();
481 double y2 = mdts.back()->globalPosition().perp();
482 double r2 = std::abs(mdts.back()->localPosition()[Trk::locR]);
483
484 double DeltaX = x2 - x1;
485 double DeltaY = y2 - y1;
486 double DistanceOfCenters = std::hypot(DeltaX, DeltaY);
487 if (DistanceOfCenters < 30) return SeedParams;
488 double Alpha0 = std::acos(DeltaX / DistanceOfCenters);
489
490 // First seed
491 double phi = mdts.front()->globalPosition().phi();
492 double RSum = r1 + r2;
493 if (RSum > DistanceOfCenters) return SeedParams;
494 double Alpha1 = std::asin(RSum / DistanceOfCenters);
495 double line_theta = Alpha0 + Alpha1;
496 double z_line = x1 + r1 * std::sin(line_theta);
497 double rho_line = y1 - r1 * std::cos(line_theta);
498
499 Amg::Vector3D gPos1(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
500 Amg::Vector3D gDir(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
501 Amg::Vector3D globalDir1(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
502 double gSlope1 = (globalDir1.perp() / globalDir1.z());
503 double gInter1 = gPos1.perp() - gSlope1 * gPos1.z();
504 double resid = SeedResiduals(mdts, gSlope1, gInter1);
505 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope1, gInter1);
506 // Second seed
507 line_theta = Alpha0 - Alpha1;
508 z_line = x1 - r1 * std::sin(line_theta);
509 rho_line = y1 + r1 * std::cos(line_theta);
510 Amg::Vector3D gPos2(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
511 Amg::Vector3D globalDir2(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
512 double gSlope2 = (globalDir2.perp() / globalDir2.z());
513 double gInter2 = gPos2.perp() - gSlope2 * gPos2.z();
514 resid = SeedResiduals(mdts, gSlope2, gInter2);
515 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope2, gInter2);
516
517 double Alpha2 = std::asin(std::abs(r2 - r1) / DistanceOfCenters);
518 if (r1 < r2) {
519 // Third seed
520 line_theta = Alpha0 + Alpha2;
521 z_line = x1 - r1 * std::sin(line_theta);
522 rho_line = y1 + r1 * std::cos(line_theta);
523
524 Amg::Vector3D gPos3(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
525 Amg::Vector3D globalDir3(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
526 double gSlope3 = (globalDir3.perp() / globalDir3.z());
527 double gInter3 = gPos3.perp() - gSlope3 * gPos3.z();
528 resid = SeedResiduals(mdts, gSlope3, gInter3);
529 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope3, gInter3);
530
531 // Fourth seed
532 line_theta = Alpha0 - Alpha2;
533 z_line = x1 + r1 * std::sin(line_theta);
534 rho_line = y1 - r1 * std::cos(line_theta);
535
536 Amg::Vector3D gPos4(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
537 Amg::Vector3D globalDir4(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
538 double gSlope4 = (globalDir4.perp() / globalDir4.z());
539 double gInter4 = gPos4.perp() - gSlope4 * gPos4.z();
540 resid = SeedResiduals(mdts, gSlope4, gInter4);
541 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope4, gInter4);
542 } else {
543 // Third seed
544 line_theta = Alpha0 + Alpha2;
545 z_line = x1 + r1 * std::sin(line_theta);
546 rho_line = y1 - r1 * std::cos(line_theta);
547
548 Amg::Vector3D gPos3(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
549 Amg::Vector3D globalDir3(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
550 double gSlope3 = (globalDir3.perp() / globalDir3.z());
551 double gInter3 = gPos3.perp() - gSlope3 * gPos3.z();
552 resid = SeedResiduals(mdts, gSlope3, gInter3);
553 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope3, gInter3);
554
555 // Fourth seed
556 line_theta = Alpha0 - Alpha2;
557 z_line = x1 - r1 * std::sin(line_theta);
558 rho_line = y1 + r1 * std::cos(line_theta);
559
560 Amg::Vector3D gPos4(rho_line * std::cos(phi), rho_line * std::sin(phi), z_line);
561 Amg::Vector3D globalDir4(std::cos(phi) * std::sin(line_theta), std::sin(phi) * std::sin(line_theta), std::cos(line_theta));
562 double gSlope4 = (globalDir4.perp() / globalDir4.z());
563 double gInter4 = gPos4.perp() - gSlope4 * gPos4.z();
564 resid = SeedResiduals(mdts, gSlope4, gInter4);
565 if (resid < m_SeedResidual) SeedParams.emplace_back(gSlope4, gInter4);
566 }
567
568 return SeedParams;
569 }
static double SeedResiduals(const std::vector< const Muon::MdtPrepData * > &mdts, double slope, double inter)
Gaudi::Property< double > m_SeedResidual
#define unlikely(x)

◆ SortMDThits()

int Muon::MSVertexTrackletTool::SortMDThits ( std::vector< std::vector< const Muon::MdtPrepData * > > & SortedMdt,
const EventContext & ctx ) const
private

Definition at line 380 of file MSVertexTrackletTool.cxx.

380 {
381 SortedMdt.clear();
382 int nMDT(0);
383
384 SG::ReadHandle<Muon::MdtPrepDataContainer> mdtTES(m_mdtTESKey, ctx);
385 if (!mdtTES.isValid()) {
386 if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Muon::MdtPrepDataContainer with key MDT_DriftCircles was not retrieved" << endmsg;
387 return 0;
388 } else {
389 if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Muon::MdtPrepDataContainer with key MDT_DriftCircles retrieved" << endmsg;
390 }
391
392 // iterators over collections, a collection corresponds to a chamber
393 for (const Muon::MdtPrepDataCollection* MDTch : *mdtTES){
394 if (MDTch->empty()) continue;
395 if (IgnoreMDTChamber(*(MDTch->begin()))) continue;
396
397 // sort per multi layer
398 std::vector<const Muon::MdtPrepData*> hitsML1;
399 std::vector<const Muon::MdtPrepData*> hitsML2;
400
401 // loop on mdt hits in the current chamber
402 for (const Muon::MdtPrepData* mdt : *MDTch) {
403 // Removes noisy hits
404 if (mdt->adc() < 50) continue;
405 // Removes dead modules or out of time hits
406 if (mdt->status() != Muon::MdtStatusDriftTime) continue;
407 // Removes tubes out of readout during drift time or with unphysical errors
408 if (mdt->localPosition()[Trk::locR] == 0.) continue;
409 if (mdt->localCovariance()(Trk::locR, Trk::locR) < 1e-6) {
410 ATH_MSG_WARNING("Found MDT with unphysical error " << m_idHelperSvc->mdtIdHelper().print_to_string(mdt->identify())
411 << " cov " << mdt->localCovariance()(Trk::locR, Trk::locR));
412 continue;
413 }
414 ++nMDT;
415
416 // sort per multi layer
417 if (m_idHelperSvc->mdtIdHelper().multilayer(mdt->identify()) == 1)
418 hitsML1.push_back(mdt);
419 else
420 hitsML2.push_back(mdt);
421
422 } // end MdtPrepDataCollection
423
424 // add
425 addMDTHits(hitsML1, SortedMdt);
426 addMDTHits(hitsML2, SortedMdt);
427 } // end MdtPrepDataContainer
428
429 return nMDT;
430 }
#define endmsg
MsgStream & msg() const
void addMDTHits(std::vector< const Muon::MdtPrepData * > &hits, std::vector< std::vector< const Muon::MdtPrepData * > > &SortedMdt) const
@ MdtStatusDriftTime
The tube produced a vaild measurement.
MuonPrepDataCollection< MdtPrepData > MdtPrepDataCollection

◆ sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )
overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, and AthCheckedComponent<::AthAlgTool >.

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ TrackletSegmentFitter()

std::vector< TrackletSegment > Muon::MSVertexTrackletTool::TrackletSegmentFitter ( const std::vector< const Muon::MdtPrepData * > & mdts) const
private

Definition at line 453 of file MSVertexTrackletTool.cxx.

453 {
454 // fits TrackletSegments from an as compatible identified set of MDT hits
455 // create the segment seeds
456 std::vector<std::pair<double, double> > SeedParams = SegSeeds(mdts);
457 // fit the segments
458 std::vector<TrackletSegment> segs = TrackletSegmentFitterCore(mdts, SeedParams);
459
460 return segs;
461 }
std::vector< std::pair< double, double > > SegSeeds(const std::vector< const Muon::MdtPrepData * > &mdts) const
std::vector< TrackletSegment > TrackletSegmentFitterCore(const std::vector< const Muon::MdtPrepData * > &mdts, const std::vector< std::pair< double, double > > &SeedParams) const

◆ TrackletSegmentFitterCore()

std::vector< TrackletSegment > Muon::MSVertexTrackletTool::TrackletSegmentFitterCore ( const std::vector< const Muon::MdtPrepData * > & mdts,
const std::vector< std::pair< double, double > > & SeedParams ) const
private

The AThTh represent the derivative fed in to the Newton minimization. If the value is too small the procedure will diverge anyway let's break the loop

Definition at line 588 of file MSVertexTrackletTool.cxx.

589 {
590 std::vector<TrackletSegment> segs;
591
592 Identifier mdtID = mdts.at(0)->identify();
593 for (const std::pair<double,double> &SeedParam : SeedParams) {
594 // Min chi^2 fit from "Precision of the ATLAS Muon Spectrometer" -- M. Woudstra
595 // http://cds.cern.ch/record/620198?ln=en (section 4.3)
596 double chi2(0);
597 double s(0), sz(0), sy(0);
598 // loop on the mdt hits, find the weighted center
599 for (const Muon::MdtPrepData *prd : mdts) {
600 // Tell clang to optimize assuming that FP exceptions can trap.
601 // Otherwise, it can vectorize the division, which can lead to
602 // spurious division-by-zero traps from unused vector lanes.
604 const double mdt_y = std::hypot(prd->globalPosition().x(), prd->globalPosition().y());
605 const double mdt_z = prd->globalPosition().z();
606 const double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
607 s += 1 / sigma2;
608 sz += mdt_z / sigma2;
609 sy += mdt_y / sigma2;
610 }
611 const double yc = sy / s;
612 const double zc = sz / s;
613
614 // Find the initial parameters of the fit
615 double alpha = std::atan2(SeedParam.first, 1.0);
616 if (alpha < 0) alpha += M_PI;
617 double dalpha = 0;
618 double d = (SeedParam.second - yc + zc * SeedParam.first) * std::cos(alpha);
619 double dd = 0;
620
621 // require segments to point to the second ML
622 if (std::abs(std::cos(alpha)) > 0.97 && (m_idHelperSvc->mdtIdHelper().isBarrel(mdtID))) continue;
623 if (std::abs(std::cos(alpha)) < 0.03 && (m_idHelperSvc->mdtIdHelper().isEndcap(mdtID))) continue;
624
625 // calculate constants used in the fit
626 double sPyy(0), sPyz(0), sPyyzz(0);
627 for (const Muon::MdtPrepData *prd : mdts) {
628 double mdt_y = std::hypot(prd->globalPosition().x(), prd->globalPosition().y());
629 double mdt_z = prd->globalPosition().z();
630 double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
631 sPyy += std::pow(mdt_y-yc,2) / sigma2;
632 sPyz += (mdt_y - yc) * (mdt_z - zc) / sigma2;
633 sPyyzz += ((mdt_y - yc) - (mdt_z - zc)) * ((mdt_y - yc) + (mdt_z - zc)) / sigma2;
634 }
635
636 // iterative fit
637 int Nitr = 0;
638 double deltaAlpha = 0;
639 double deltad = 0;
640 while (true) {
641 double sumRyi(0), sumRzi(0), sumRi(0);
642 chi2 = 0;
643 ++Nitr;
644 const double cos_a = std::cos(alpha);
645 const double sin_a = std::sin(alpha);
646 for (const Muon::MdtPrepData *prd : mdts) {
647 double mdt_y = prd->globalPosition().perp();
648 double mdt_z = prd->globalPosition().z();
649 double yPi = -(mdt_z - zc) * sin_a + (mdt_y - yc) * cos_a - d;
650 double signR = yPi >= 0 ? -1. : 1;
651 double sigma2 = std::pow(Amg::error(prd->localCovariance(), Trk::locR),2);
652 double ri = signR * prd->localPosition()[Trk::locR];
653 sumRyi += ri * (mdt_y - yc) / sigma2;
654 sumRzi += ri * (mdt_z - zc) / sigma2;
655 sumRi += ri / sigma2;
656 chi2 += std::pow(yPi+ri,2) / sigma2;
657 }
658 double bAlpha = -1 * sPyz + cos_a * (sin_a * sPyyzz + 2 * cos_a * sPyz + sumRzi) + sin_a * sumRyi;
659 double AThTh = sPyy + cos_a * (2 * sin_a * sPyz - cos_a * sPyyzz);
663 if (std::abs(AThTh) < 1.e-7) break;
664 // the new alpha & d parameters
665 double alphaNew = alpha + bAlpha / AThTh;
666 double dNew = sumRi / s;
667 // the errors
668 dalpha = std::sqrt(1 / std::abs(AThTh));
669 dd = std::sqrt(1 / s);
670 deltaAlpha = std::abs(alphaNew - alpha);
671 deltad = std::abs(d - dNew);
672 // test if the new segment is different than the previous
673 if (deltaAlpha < 5.e-7 && deltad < 5.e-6) break;
674 alpha = alphaNew;
675 d = dNew;
676 // Guard against infinite loops
677 if (Nitr > 10) break;
678 } // end while loop
679
680 // find the chi^2 probability of the segment
681 double chi2Prob = TMath::Prob(chi2, mdts.size() - 2);
682 // keep only "good" segments
683 if (chi2Prob > m_minSegFinderChi2) {
684 double z0 = zc - d * std::sin(alpha);
685 double dz0 = std::hypot(dd*std::sin(alpha), d*dalpha*std::cos(alpha));
686 double y0 = yc + d * std::cos(alpha);
687 double dy0 = std::hypot(dd*std::cos(alpha), d*dalpha*std::sin(alpha));
688 // find the hit pattern, which side of the wire did the particle pass? (1==Left, 2==Right)
689 /*
690 ( )(/O)( )
691 (./)( )( ) == RRL == 221
692 (O/)( )( )
693 */
694 int pattern(0);
695 if (mdts.size() > 8)
696 pattern = -1; // with more then 8 MDTs the pattern is unique
697 else {
698 for (unsigned int k = 0; k < mdts.size(); ++k) {
699 int base = std::pow(10, k);
700 double mdtR = std::hypot(mdts.at(k)->globalPosition().x(), mdts.at(k)->globalPosition().y());
701 double mdtZ = mdts.at(k)->globalPosition().z();
702 double zTest = (mdtR - y0) / std::tan(alpha) + z0 - mdtZ;
703 if (zTest > 0)
704 pattern += 2 * base;
705 else
706 pattern += base;
707 }
708 }
709
710 // find the position of the tracklet in the global frame
711 double mdtPhi = mdts.at(0)->globalPosition().phi();
712 Amg::Vector3D segpos(y0 * std::cos(mdtPhi), y0 * std::sin(mdtPhi), z0);
713 // create the tracklet
714 TrackletSegment MyTrackletSegment{m_idHelperSvc.get(), mdts, segpos, alpha, dalpha, dy0, dz0, pattern};
715 segs.push_back(MyTrackletSegment);
716 if (pattern == -1) break; // stop if we find a segment with more than 8 hits (guaranteed to be unique!)
717 }
718 } // end loop on segment seeds
719
720 // in case more than 1 segment is reconstructed, check if there are duplicates using the hit patterns
721 if (segs.size() > 1) {
722 std::vector<TrackletSegment> tmpSegs;
723 for (unsigned int i1 = 0; i1 < segs.size(); ++i1) {
724 bool isUnique = true;
725 int pattern1 = segs.at(i1).getHitPattern();
726 for (unsigned int i2 = (i1 + 1); i2 < segs.size(); ++i2) {
727 if (pattern1 == -1) break;
728 int pattern2 = segs.at(i2).getHitPattern();
729 if (pattern1 == pattern2) isUnique = false;
730 }
731 if (isUnique) tmpSegs.push_back(segs.at(i1));
732 }
733 segs = tmpSegs;
734 }
735
736 // return the unique segments
737 return segs;
738 }
#define M_PI
static Double_t sz
Gaudi::Property< double > m_minSegFinderChi2
double chi2(TH1 *h0, TH1 *h1)
std::string base
Definition hcg.cxx:83
#define CXXUTILS_TRAPPING_FP
Definition trapping_fp.h:24

◆ TrackMomentum()

double Muon::MSVertexTrackletTool::TrackMomentum ( const Identifier trkID,
const double deltaAlpha ) const
private

Definition at line 924 of file MSVertexTrackletTool.cxx.

924 {
925 // p = k/delta_alpha
926 bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
927 bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
928 int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
929
930 double dalpha = std::abs(deltaAlpha);
931 double pTot = m_straightTrackletpTot;
932 if (isBarrel){
933 if (stationRegion == 0){
934 if (isSmall) pTot = m_straightTrackletpTot; // default value for BIS
935 else pTot = c_BIL / dalpha;
936 }
937 else if (stationRegion == 2){
938 if (isSmall) pTot = c_BMS / dalpha;
939 else pTot = c_BML / dalpha;
940 }
941 else if (stationRegion == 3){
942 if (isSmall) pTot = m_straightTrackletpTot; // default value for BOS
943 else pTot = c_BOL / dalpha;
944 }
945 }
946
947 if (pTot > m_maxpTot) pTot = m_straightTrackletpTot;
948
949 return pTot;
950 }
Gaudi::Property< double > m_straightTrackletpTot

◆ TrackMomentumError() [1/2]

double Muon::MSVertexTrackletTool::TrackMomentumError ( const TrackletSegment & ml1) const
private

Definition at line 983 of file MSVertexTrackletTool.cxx.

983 {
984 // uncertainty in 1/p
985 const Identifier trkID = ml1.getIdentifier();
986 bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
987 bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
988 int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
989
990 double dalpha = std::abs(ml1.alphaError());
991 double pErr = dalpha / c_BML;
992 if (isBarrel){
993 if (stationRegion == 0){
994 if (isSmall) pErr = dalpha / c_BML; // default value for BIS
995 else pErr = dalpha / c_BIL;
996 }
997 else if (stationRegion == 2){
998 if (isSmall) pErr = dalpha / c_BMS;
999 else pErr = dalpha / c_BML;
1000 }
1001 else if (stationRegion == 3){
1002 if (isSmall) pErr = dalpha / c_BML; // default value for BOS
1003 else pErr = dalpha / c_BOL;
1004 }
1005 }
1006
1007 return pErr;
1008 }

◆ TrackMomentumError() [2/2]

double Muon::MSVertexTrackletTool::TrackMomentumError ( const TrackletSegment & ml1,
const TrackletSegment & ml2 ) const
private

Definition at line 954 of file MSVertexTrackletTool.cxx.

954 {
955 // uncertainty on 1/p
956 const Identifier trkID = ml1.getIdentifier();
957 bool isBarrel = m_idHelperSvc->mdtIdHelper().isBarrel(trkID);
958 bool isSmall = m_idHelperSvc->mdtIdHelper().isSmall(trkID);
959 int stationRegion = m_idHelperSvc->mdtIdHelper().stationRegion(trkID);
960
961 double dalpha = std::hypot(ml1.alphaError(), ml2.alphaError());
962 double pErr = dalpha / c_BML;
963 if (isBarrel){
964 if (stationRegion == 0){
965 if (isSmall) pErr = dalpha / c_BML; // default value for BIS
966 else pErr = dalpha / c_BIL;
967 }
968 else if (stationRegion == 2){
969 if (isSmall) pErr = dalpha / c_BMS;
970 else pErr = dalpha / c_BML;
971 }
972 else if (stationRegion == 3){
973 if (isSmall) pErr = dalpha / c_BML; // default value for BOS
974 else pErr = dalpha / c_BOL;
975 }
976 }
977
978 return pErr;
979 }

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 298 of file AthCommonDataStore.h.

298 {
299 for (auto &a : m_vhka) {
301 for (auto k : keys) {
302 k->setOwner(this);
303 }
304 }
305 }
std::vector< SG::VarHandleKeyArray * > m_vhka

Member Data Documentation

◆ m_BarrelDeltaAlphaCut

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_BarrelDeltaAlphaCut {this, "BarrelDeltaAlphaCut", 0.2*Gaudi::Units::radian, "maximum delta_alpha allowed in barrel MS chambers"}
private

Definition at line 59 of file MSVertexTrackletTool.h.

59{this, "BarrelDeltaAlphaCut", 0.2*Gaudi::Units::radian, "maximum delta_alpha allowed in barrel MS chambers"};

◆ m_d12_max

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_d12_max {this, "d12_max", 50.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 2"}
private

Definition at line 52 of file MSVertexTrackletTool.h.

52{this, "d12_max", 50.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 2"};

◆ m_d13_max

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_d13_max {this, "d13_max", 80.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 3"}
private

Definition at line 53 of file MSVertexTrackletTool.h.

53{this, "d13_max", 80.*Gaudi::Units::millimeter, "max separation in mm between hits in tube 1 and 3"};

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default).

Definition at line 381 of file AthCommonDataStore.h.

◆ m_EndcapDeltaAlphaCut

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_EndcapDeltaAlphaCut {this, "EndcapDeltaAlphaCut", 0.015*Gaudi::Units::radian, "maximum delta_alpha allowed in the endcap MS chambers"}
private

Definition at line 60 of file MSVertexTrackletTool.h.

60{this, "EndcapDeltaAlphaCut", 0.015*Gaudi::Units::radian, "maximum delta_alpha allowed in the endcap MS chambers"};

◆ m_errorCutOff

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_errorCutOff {this, "errorCutOff", 0.001, "minimal hit error"}
private

Definition at line 54 of file MSVertexTrackletTool.h.

54{this, "errorCutOff", 0.001, "minimal hit error"};

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default).

Definition at line 378 of file AthCommonDataStore.h.

◆ m_idHelperSvc

ServiceHandle<Muon::IMuonIdHelperSvc> Muon::MSVertexTrackletTool::m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}
private

Definition at line 48 of file MSVertexTrackletTool.h.

48{this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};

◆ m_maxDeltabCut

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_maxDeltabCut {this, "maxDeltabCut", 3*Gaudi::Units::millimeter, "maximum delta_b allowed"}
private

Definition at line 61 of file MSVertexTrackletTool.h.

61{this, "maxDeltabCut", 3*Gaudi::Units::millimeter, "maximum delta_b allowed"};

◆ m_maxpTot

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_maxpTot {this, "maxpTot", 10000.*Gaudi::Units::MeV, "maximum measurable total momentum in MeV beyond which tracklets are assumed to be straight"}
private

Definition at line 63 of file MSVertexTrackletTool.h.

63{this, "maxpTot", 10000.*Gaudi::Units::MeV, "maximum measurable total momentum in MeV beyond which tracklets are assumed to be straight"};

◆ m_mdtTESKey

SG::ReadHandleKey<Muon::MdtPrepDataContainer> Muon::MSVertexTrackletTool::m_mdtTESKey {this, "mdtTES", "MDT_DriftCircles"}
private

Definition at line 45 of file MSVertexTrackletTool.h.

45{this, "mdtTES", "MDT_DriftCircles"};

◆ m_minpTot

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_minpTot {this, "minpTot", 800.*Gaudi::Units::MeV, "minimum measurable total momentum in MeV"}
private

Definition at line 62 of file MSVertexTrackletTool.h.

62{this, "minpTot", 800.*Gaudi::Units::MeV, "minimum measurable total momentum in MeV"};

◆ m_minSegFinderChi2

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_minSegFinderChi2 {this, "MinSegFinderChi2Prob", 0.05, "tracklet segment fitter chi^2 probability cut"}
private

Definition at line 57 of file MSVertexTrackletTool.h.

57{this, "MinSegFinderChi2Prob", 0.05, "tracklet segment fitter chi^2 probability cut"};

◆ m_SeedResidual

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_SeedResidual {this, "SeedResidual", 5., "max residual for tracklet seeds"}
private

Definition at line 56 of file MSVertexTrackletTool.h.

56{this, "SeedResidual", 5., "max residual for tracklet seeds"};

◆ m_straightTrackletInvPerr

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_straightTrackletInvPerr {this, "straightTrackletInvPerr", 5.0e-5/Gaudi::Units::MeV, "error in the inverse momentum in MeV^-1 assigned to straight tracklets"}
private

Definition at line 65 of file MSVertexTrackletTool.h.

65{this, "straightTrackletInvPerr", 5.0e-5/Gaudi::Units::MeV, "error in the inverse momentum in MeV^-1 assigned to straight tracklets"};

◆ m_straightTrackletpTot

Gaudi::Property<double> Muon::MSVertexTrackletTool::m_straightTrackletpTot {this, "straightTrackletpTot", 1.0e5*Gaudi::Units::MeV, "total momentum in MeV assigned to straight tracklets"}
private

Definition at line 64 of file MSVertexTrackletTool.h.

64{this, "straightTrackletpTot", 1.0e5*Gaudi::Units::MeV, "total momentum in MeV assigned to straight tracklets"};

◆ m_tightTrackletRequirement

Gaudi::Property<bool> Muon::MSVertexTrackletTool::m_tightTrackletRequirement {this, "tightTrackletRequirement", false, "tight tracklet requirement (affects efficiency - disabled by default)"}
private

Definition at line 66 of file MSVertexTrackletTool.h.

66{this, "tightTrackletRequirement", false, "tight tracklet requirement (affects efficiency - disabled by default)"};

◆ m_TPContainer

SG::WriteHandleKey<xAOD::TrackParticleContainer> Muon::MSVertexTrackletTool::m_TPContainer {this, "xAODTrackParticleContainer", "MSonlyTracklets"}
private

Definition at line 46 of file MSVertexTrackletTool.h.

46{this, "xAODTrackParticleContainer", "MSonlyTracklets"};

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared
privateinherited

Definition at line 387 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka
privateinherited

Definition at line 386 of file AthCommonDataStore.h.


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