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< StoreGateSvcStoreGateSvc_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< TrackletSegmentTrackletSegmentFitter (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< TrackletSegmentTrackletSegmentFitterCore (const std::vector< const Muon::MdtPrepData * > &mdts, const std::vector< std::pair< double, double > > &SeedParams) const
std::vector< TrackletSegmentCleanSegments (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< TrackletResolveAmbiguousTracklets (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::MdtPrepDataContainerm_mdtTESKey {this, "mdtTES", "MDT_DriftCircles"}
SG::WriteHandleKey< xAOD::TrackParticleContainerm_TPContainer {this, "xAODTrackParticleContainer", "MSonlyTracklets"}
ServiceHandle< Muon::IMuonIdHelperSvcm_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 388 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 738 of file MSVertexTrackletTool.cxx.

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

887 {
888 double ChMid = (ML1seg.getChMidPoint() + ML2seg.getChMidPoint()) / 2.0;
889 // Calculate the Delta b (see http://inspirehep.net/record/1266438)
890 double mid1(100), mid2(1000);
891 double deltab(100);
892 if (m_idHelperSvc->mdtIdHelper().isBarrel(ML1seg.getIdentifier())) {
893 // delta b in the barrel
894 mid1 = (ChMid - ML1seg.globalPosition().perp()) / std::tan(ML1seg.alpha()) + ML1seg.globalPosition().z();
895 mid2 = (ChMid - ML2seg.globalPosition().perp()) / std::tan(ML2seg.alpha()) + ML2seg.globalPosition().z();
896 double r01 = ML1seg.globalPosition().perp() - ML1seg.globalPosition().z() * std::tan(ML1seg.alpha());
897 double r02 = ML2seg.globalPosition().perp() - ML2seg.globalPosition().z() * std::tan(ML2seg.alpha());
898 deltab = (mid2 * std::tan(ML1seg.alpha()) - ChMid + r01) / (std::hypot(1,std::tan(ML1seg.alpha())));
899 double deltab2 = (mid1 * std::tan(ML2seg.alpha()) - ChMid + r02) / (std::hypot(1,std::tan(ML2seg.alpha())));
900 if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
901 } else {
902 // delta b in the endcap
903 mid1 = ML1seg.globalPosition().perp() + std::tan(ML1seg.alpha()) * (ChMid - ML1seg.globalPosition().z());
904 mid2 = ML2seg.globalPosition().perp() + std::tan(ML2seg.alpha()) * (ChMid - ML2seg.globalPosition().z());
905 double z01 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
906 double z02 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp() / std::tan(ML1seg.alpha());
907 deltab = (mid2 / std::tan(ML1seg.alpha()) - ChMid + z01) / (std::hypot(1,1/std::tan(ML1seg.alpha())));
908 double deltab2 = (mid1 / std::tan(ML2seg.alpha()) - ChMid + z02) / (std::hypot(1,1/std::tan(ML2seg.alpha())));
909 if (std::abs(deltab2) < std::abs(deltab)) deltab = deltab2;
910 }
911
912 // calculate the maximum allowed Delta b based on delta alpha uncertainties and ML spacing
913 double dbmax = 5 * std::abs(ChMid - ML1seg.getChMidPoint()) * std::hypot(ML1seg.alphaError(), ML2seg.alphaError());
914 if (dbmax > m_maxDeltabCut) dbmax = m_maxDeltabCut;
915 return std::abs(deltab) < dbmax;
916 }
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_WARNING(x)
#define ATH_MSG_DEBUG(x)
double charge(const T &p)
Definition AtlasPID.h:997
if(febId1==febId2)
#define y
#define x
static const Attributes_t empty
constexpr int pow(int base, int exp) noexcept
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 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
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 364 of file AthCommonDataStore.h.

364 {
366 }

◆ ResolveAmbiguousTracklets()

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

Definition at line 1008 of file MSVertexTrackletTool.cxx.

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

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

◆ 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 >, AthCheckedComponent<::AthAlgTool >, and DerivationFramework::CfAthAlgTool.

◆ 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 584 of file MSVertexTrackletTool.cxx.

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

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

◆ TrackMomentumError() [1/2]

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

Definition at line 979 of file MSVertexTrackletTool.cxx.

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

◆ TrackMomentumError() [2/2]

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

Definition at line 950 of file MSVertexTrackletTool.cxx.

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

◆ updateVHKA()

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

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }
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 393 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 390 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 399 of file AthCommonDataStore.h.

◆ m_vhka

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

Definition at line 398 of file AthCommonDataStore.h.


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