From 264f5c81679126cbc8311ee5a991d5a33ae2d179 Mon Sep 17 00:00:00 2001
From: Andrea Coccaro <Andrea.Coccaro@cern.ch>
Date: Tue, 16 Sep 2014 13:46:14 +0200
Subject: [PATCH] removing the interface method in the tool header files
(MSVertexTools-00-01-08)
---
.../MSVertexTools/cmt/requirements | 39 +
.../MSVertexTools/src/MSVertexRecoTool.cxx | 1148 +++++++++++++++
.../MSVertexTools/src/MSVertexRecoTool.h | 104 ++
.../src/MSVertexTrackletTool.cxx | 1295 +++++++++++++++++
.../MSVertexTools/src/MSVertexTrackletTool.h | 88 ++
.../src/components/MSVertexTools_entries.cxx | 19 +
.../src/components/MSVertexTools_load.cxx | 5 +
7 files changed, 2698 insertions(+)
create mode 100644 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/cmt/requirements
create mode 100755 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.cxx
create mode 100755 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.h
create mode 100755 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.cxx
create mode 100755 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.h
create mode 100644 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_entries.cxx
create mode 100644 MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_load.cxx
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/cmt/requirements b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/cmt/requirements
new file mode 100644
index 0000000000000..fa7d7598bf4a7
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/cmt/requirements
@@ -0,0 +1,39 @@
+package MSVertexTools
+
+author Daniel Ventura <ventura@cern.ch>
+
+# ============================================================================================
+public
+
+use AtlasPolicy AtlasPolicy-*
+
+# ============================================================================================
+private
+
+use AthenaKernel AthenaKernel-* Control
+use AthenaBaseComps AthenaBaseComps-* Control
+use DataModel DataModel-* Control
+use StoreGate StoreGate-* Control
+use xAODTracking xAODTracking-* Event/xAOD
+use GaudiInterface GaudiInterface-* External
+use AtlasCLHEP AtlasCLHEP-* External
+use Identifier Identifier-* DetectorDescription
+use GeoAdaptors GeoAdaptors-* DetectorDescription/GeoModel
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+use EventPrimitives EventPrimitives-* Event
+use MuonIdHelpers MuonIdHelpers-* MuonSpectrometer
+use MuonReadoutGeometry MuonReadoutGeometry-* MuonSpectrometer/MuonDetDescr
+use MSVertexUtils MSVertexUtils-* MuonSpectrometer/MSVertexReconstruction
+use MSVertexToolInterfaces MSVertexToolInterfaces-* MuonSpectrometer/MSVertexReconstruction
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
+
+# ============================================================================================
+public
+
+apply_pattern component_library
+library MSVertexTools *.cxx components/*.cxx
+
+#private
+#macro cppdebugflags '$(cppdebugflags_s)'
+#macro_remove componentshr_linkopts "-Wl,-s"
\ No newline at end of file
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.cxx b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.cxx
new file mode 100755
index 0000000000000..c19c5ed07a703
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.cxx
@@ -0,0 +1,1148 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MSVertexRecoTool.h"
+#include "GaudiKernel/MsgStream.h"
+#include "EventPrimitives/EventPrimitivesHelpers.h"
+#include "xAODTracking/Vertex.h"
+#include "CLHEP/Random/RandFlat.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/VertexAuxContainer.h"
+#include "xAODTracking/TrackParticle.h"
+
+#define MAXPLANES 100
+
+/*
+ MS vertex reconstruction routine
+ See documentation at https://cds.cern.ch/record/1455664 and https://cds.cern.ch/record/1520894
+ Takes Tracklets as input and creates vertices in the MS volume
+*/
+
+namespace Muon {
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertexRecoTool::MSVertexRecoTool (const std::string& type, const std::string& name,
+ const IInterface* parent)
+ :
+ AthAlgTool(type, name, parent),
+ m_extrapolator("Trk::Extrapolator/AtlasExtrapolator") {
+ declareInterface<IMSVertexRecoTool>(this);
+
+ declareProperty("xAODVertexContainer",xAODContainerName="MSDisplacedVertex");//name of SG container
+
+ // nominal phi angle for tracklets
+ declareProperty("TrackPhiAngle",m_TrackPhiAngle = 0.0);
+ // chi^2 probability cut
+ declareProperty("VxChi2Probability",VxChi2ProbCUT = 0.05);
+ // distance between two adjacent planes
+ declareProperty("VertexPlaneDist",m_VxPlaneDist = 200.);
+ // position of last radial plane
+ declareProperty("VertexMaxRadialPlane",m_VertexMaxRadialPlane = 7000.);
+ // position of first radial plane
+ declareProperty("VertexMinRadialPlane",m_VertexMinRadialPlane = 3500.);
+ // minimum occupancy to be considered "high occupancy"
+ declareProperty("MinimumHighOccupancy",m_ChamberOccupancyMin = 0.25);
+ // number of high occupancy chambers required to be signal like
+ declareProperty("MinimumNumberOfHighOccupancy",m_minHighOccupancyChambers = 2);
+ declareProperty("UseOldMSVxEndcapMethod",m_useOldMSVxEndcapMethod = false);
+ declareProperty("MaxTollDist",m_MaxTollDist = 300);
+
+ // options to calculate vertex systematic uncertainty from the tracklet reco uncertainty
+ declareProperty("DoSystematicUncertainty",m_doSystematics = false);
+ declareProperty("BarrelTrackletUncertainty",m_BarrelTrackletUncert = 0.1);
+ declareProperty("EndcapTrackletUncertainty",m_EndcapTrackletUncert = 0.1);
+ declareProperty("RndmEngine",m_rndmEngineName,"Random track killing parameter: engine used");
+
+ //extrapolator
+ declareProperty("MyExtrapolator", m_extrapolator );
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexRecoTool::initialize() {
+
+ if( AthAlgTool::initialize().isFailure() ) {
+ ATH_MSG_ERROR( "Failed to initialize AthAlgTool" );
+ return StatusCode::FAILURE;
+ }
+
+ if(detStore()->retrieve(m_mdtIdHelper,"MDTIDHELPER").isFailure()) {
+ ATH_MSG_ERROR( "Failed to retrieve the mdtIdHelper" );
+ return StatusCode::FAILURE;
+ }
+
+ if(detStore()->retrieve(m_rpcIdHelper,"RPCIDHELPER").isFailure()) {
+ ATH_MSG_ERROR( "Failed to retrieve the rpcIdHelper" );
+ return StatusCode::FAILURE;
+ }
+
+ if(detStore()->retrieve(m_tgcIdHelper,"TGCIDHELPER").isFailure()) {
+ ATH_MSG_ERROR( "Failed to retrieve the tgcIdHelper" );
+ return StatusCode::FAILURE;
+ }
+
+ if(m_doSystematics) {
+ m_rndmEngine = m_rndmSvc->GetEngine("TrackletKiller");
+ if(m_rndmEngine==0) {
+ ATH_MSG_FATAL( "Could not get RndmEngine" );
+ return StatusCode::FAILURE;
+ }
+ }
+
+ if(m_extrapolator.retrieve().isFailure()) {
+ ATH_MSG_FATAL( "Extrapolator could not be retrieved" );
+ return StatusCode::FAILURE;
+ }
+
+ PI = 3.1415927;
+
+ return StatusCode::SUCCESS;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertexRecoTool::~MSVertexRecoTool() {}
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexRecoTool::findMSvertices(std::vector<Tracklet>& tracklets) {
+
+ //if there are less than 2 tracks, return (vertexing not possible)
+ if(tracklets.size() < 3 || tracklets.size() > 40) return StatusCode::SUCCESS;
+
+ std::vector<MSVertex*> vertices;
+
+ //group the tracks
+ std::vector<Tracklet> BarrelTracklets;
+ std::vector<Tracklet> EndcapTracklets;
+ for(unsigned int i=0; i<tracklets.size(); ++i) {
+ if(tracklets.at(i).mdtChamber() <= 11 || tracklets.at(i).mdtChamber() == 52) BarrelTracklets.push_back(tracklets.at(i));
+ else EndcapTracklets.push_back(tracklets.at(i));
+ }
+
+ //find any clusters of tracks & decide if tracks are from single muon
+ std::vector<Muon::MSVertexRecoTool::TrkCluster> BarrelClusters = findTrackClusters(BarrelTracklets);
+ std::vector<Muon::MSVertexRecoTool::TrkCluster> EndcapClusters = findTrackClusters(EndcapTracklets);
+
+ //if doSystematics, remove tracklets according to the tracklet reco uncertainty and rerun the cluster finder
+ if(m_doSystematics) {
+ std::vector<Tracklet> BarrelSystTracklets,EndcapSystTracklets;
+ for(unsigned int i=0; i<BarrelTracklets.size(); ++i) {
+ float prob = CLHEP::RandFlat::shoot( m_rndmEngine, 0, 1 );
+ if(prob > m_BarrelTrackletUncert) BarrelSystTracklets.push_back(BarrelTracklets.at(i));
+ }
+ if(BarrelSystTracklets.size() >= 3) {
+ std::vector<Muon::MSVertexRecoTool::TrkCluster> BarrelSystClusters = findTrackClusters(BarrelSystTracklets);
+ for(unsigned int i=0; i<BarrelSystClusters.size(); ++i) {
+ BarrelSystClusters.at(i).isSystematic = true;
+ BarrelClusters.push_back(BarrelSystClusters.at(i));
+ }
+ }
+ for(unsigned int i=0; i<EndcapTracklets.size(); ++i) {
+ float prob = CLHEP::RandFlat::shoot( m_rndmEngine, 0, 1 );
+ if(prob > m_EndcapTrackletUncert) EndcapSystTracklets.push_back(EndcapTracklets.at(i));
+ }
+ if(EndcapSystTracklets.size() >= 3) {
+ std::vector<Muon::MSVertexRecoTool::TrkCluster> EndcapSystClusters = findTrackClusters(EndcapSystTracklets);
+ for(unsigned int i=0; i<EndcapSystClusters.size(); ++i) {
+ EndcapSystClusters.at(i).isSystematic = true;
+ EndcapClusters.push_back(EndcapSystClusters.at(i));
+ }
+ }
+ }
+
+ ///////////////////////////////////////////VERTEX ROUTINES/////////////////////////////////////////
+ //find vertices in the barrel MS (vertices using barrel tracklets)
+ for(unsigned int i=0; i<BarrelClusters.size(); ++i) {
+ if(BarrelClusters[i].ntrks < 3) continue;
+ ATH_MSG_DEBUG( "Attempting to build vertex from " << BarrelClusters[i].ntrks << " tracklets in the barrel" );
+ MSVertex* barvertex = MSVxFinder(BarrelClusters[i].tracks);
+ if(!barvertex) continue;
+ if(barvertex->getChi2Probability() > 0.05) {
+ HitCounter(barvertex);
+ if(barvertex->getNMDT() > 250 && (barvertex->getNRPC()+barvertex->getNTGC()) > 200) {
+ ATH_MSG_DEBUG( "Vertex found in the barrel with n_trk = " << barvertex->getNTracks() << " located at (eta,phi) = ("
+ << barvertex->getPosition().eta() << ", " << barvertex->getPosition().phi() << ")" );
+ if(BarrelClusters[i].isSystematic) barvertex->setAuthor(3);
+ vertices.push_back(barvertex);
+ }//end minimum good vertex criteria
+ }
+ }//end loop on barrel tracklet clusters
+
+ //find vertices in the endcap MS (vertices using endcap tracklets)
+ for(unsigned int i=0; i<EndcapClusters.size(); ++i) {
+ if(EndcapClusters[i].ntrks < 3) continue;
+ ATH_MSG_DEBUG( "Attempting to build vertex from " << EndcapClusters[i].ntrks << " tracklets in the endcap" );
+
+ MSVertex* endvertex;
+ if(m_useOldMSVxEndcapMethod)
+ endvertex = MSStraightLineVx_oldMethod(EndcapClusters[i].tracks);
+ else
+ endvertex = MSStraightLineVx(EndcapClusters[i].tracks);
+
+ if(!endvertex) continue;
+ if(endvertex->getPosition().perp() < 10000 && fabs(endvertex->getPosition().z()) < 14000 &&
+ fabs(endvertex->getPosition().z()) > 8000 && endvertex->getNTracks() >= 3) {
+ HitCounter(endvertex);
+ if(endvertex->getNMDT() > 250 && (endvertex->getNRPC()+endvertex->getNTGC()) > 200) {
+ ATH_MSG_DEBUG( "Vertex found in the endcap with n_trk = " << endvertex->getNTracks() << " located at (eta,phi) = ("
+ << endvertex->getPosition().eta() << ", " << endvertex->getPosition().phi() << ")" );
+ if(EndcapClusters[i].isSystematic) endvertex->setAuthor(4);
+ vertices.push_back(endvertex);
+ }//end minimum good vertex criteria
+ }
+ }//end loop on endcap tracklet clusters
+
+ if(vertices.size() == 0) return StatusCode::SUCCESS;
+
+ //convert to xAOD::vertex & cleanup
+ xAOD::VertexContainer* xAODVxContainer = new xAOD::VertexContainer();
+ xAOD::VertexAuxContainer* aux = new xAOD::VertexAuxContainer();
+ xAODVxContainer->setStore( aux );
+ //store the new vertices!
+ CHECK( evtStore()->record( xAODVxContainer, xAODContainerName ) );
+ CHECK( evtStore()->record( aux, xAODContainerName + "Aux" ) );
+ //fill the containers
+ for(std::vector<MSVertex*>::const_iterator vxIt=vertices.begin(); vxIt!=vertices.end(); ++vxIt) {
+ xAOD::Vertex* xAODVx = new xAOD::Vertex();
+ xAODVx->setVertexType( xAOD::VxType::SecVtx );
+ xAODVx->setPosition( (*vxIt)->getPosition() );
+ xAODVx->setFitQuality( (*vxIt)->getChi2(), (*vxIt)->getNTracks()-1 );
+
+ //dress the vertex with the hit counts
+ xAODVx->auxdata< int >( "nMDT" ) = (*vxIt)->getNMDT();
+ xAODVx->auxdata< int >( "nRPC" ) = (*vxIt)->getNRPC();
+ xAODVx->auxdata< int >( "nTGC" ) = (*vxIt)->getNTGC();
+
+ //TODO: Add links to the tracks -- need to think how to proceed .... example code below
+ //unsigned int VTAVsize = vxVertex.vxTrackAtVertex()->size();
+ //for (unsigned int i = 0 ; i < VTAVsize ; ++i) {
+ // Trk::VxTrackAtVertex* VTAV = vxVertex.vxTrackAtVertex()->at(i);
+ // Trk::ITrackLink* trklink = VTAV->trackOrParticleLink();
+ // //do we need really to check this? or can we assume xAOD::TrackParticle objects to be there?
+ // // definitely we don't want those to be Rec::TrackParticle do avoid mixing up things
+ // Trk::LinkToXAODTrackParticle* linkToXAODTP = dynamic_cast<Trk::LinkToXAODTrackParticle*>(trklink);
+ // if (!linkToXAODTP) {
+ // ATH_MSG_WARNING ("Skipping track. Trying to convert to xAOD::Vertex a Trk::VxCandidate with links to something else than xAOD::TrackParticle.");
+ // } else {
+ // //Now set the newlink to the new xAOD vertex
+ // xAODVx->addTrackAtVertex(*linkToXAODTP, VTAV->vtxCompatibility());
+ // }
+ //}
+
+ //store the new xAOD vertex
+ xAODVxContainer->push_back(xAODVx);
+
+ //delete the internal EDM vertex
+ delete (*vxIt);
+ }
+
+
+ return StatusCode::SUCCESS;
+ }//end find vertices
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexRecoTool::finalize() {
+ return StatusCode::SUCCESS;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ Muon::MSVertexRecoTool::TrkCluster MSVertexRecoTool::ClusterizeTracks(std::vector<Tracklet>& tracks) {
+
+ if(tracks.size() > 40) {
+ ATH_MSG_DEBUG( "Too many tracks found, results probably not reliable! returning empty cluster" );
+ TrkCluster emptycluster;
+ emptycluster.ntrks=0;
+ emptycluster.eta=-99999.;
+ emptycluster.phi=-99999.;
+ return emptycluster;
+ }
+ TrkCluster trkClu[50];
+ TrkCluster trkClu0[50];
+ trkClu[0].ntrks = 0;
+ trkClu[0].eta = -10;
+ trkClu[0].phi = -10;
+ int ncluster(0);
+ //use each tracklet as a seed for the clusters
+ for(std::vector<Tracklet>::iterator trkItr=tracks.begin(); trkItr!=tracks.end(); ++trkItr) {
+ TrkCluster clu;
+ clu.eta = trkItr->globalPosition().eta();
+ clu.phi = trkItr->globalPosition().phi();
+ clu.ntrks = 0;
+ for(unsigned int i=0; i<tracks.size(); ++i) clu.trks[i]=0;
+ trkClu[ncluster] = clu;
+ trkClu0[ncluster] = clu;
+ ++ncluster;
+ if(ncluster >= 49) {
+ TrkCluster emptycluster;
+ return emptycluster;
+ }
+ }
+ //loop on the clusters and let the center move to find the optimal cluster centers
+ for(int icl=0; icl<ncluster; ++icl) {
+ bool improvement = true;
+ int nitr(0);
+ while(improvement) {
+ if(nitr > 6) break;
+ int ntracks(0);
+ for(int jcl=0; jcl<ncluster; ++jcl) {
+ float dEta = trkClu[icl].eta - trkClu0[jcl].eta;
+ float dPhi = trkClu[icl].phi - trkClu0[jcl].phi;
+ while(fabs(dPhi) > PI) {
+ if(dPhi < 0) dPhi += 2*PI;
+ else dPhi -= 2*PI;
+ }
+ if(fabs(dEta) < 0.7 && fabs(dPhi) < PI/3.) {
+ ntracks++;
+ trkClu[icl].eta = trkClu[icl].eta - dEta/ntracks;
+ trkClu[icl].phi = trkClu[icl].phi - dPhi/ntracks;
+ while(fabs(trkClu[icl].phi) > PI) {
+ if(trkClu[icl].phi > 0) trkClu[icl].phi -= 2*PI;
+ else trkClu[icl].phi += 2*PI;
+ }
+ }
+ }//end jcl loop
+ //find the number of tracks in the new cluster
+ int itracks[100];
+ for(int k=0; k<ncluster; ++k) itracks[k]=0;
+ int ntracks2(0);
+ for(int jcl=0; jcl<ncluster; ++jcl) {
+ float dEta = fabs(trkClu[icl].eta - trkClu0[jcl].eta);
+ float dPhi = trkClu[icl].phi - trkClu0[jcl].phi;
+ while(fabs(dPhi) > PI) {
+ if(dPhi < 0) dPhi += 2*PI;
+ else dPhi -= 2*PI;
+ }
+ if(dEta < 0.7 && fabs(dPhi) < PI/3.) {
+ ntracks2++;
+ itracks[jcl] = 1;
+ }
+ }//end jcl loop
+ if(ntracks2 > trkClu[icl].ntrks) {
+ trkClu[icl].ntrks = ntracks2;
+ for(int k=0; k<ncluster; ++k) {
+ trkClu[icl].trks[k] = itracks[k];
+ }
+ }
+ else improvement = false;
+ nitr++;
+ }//end while
+ }//end icl loop
+ //find the best cluster
+ TrkCluster BestCluster = trkClu[0];
+ for(int icl=1; icl<ncluster; ++icl) {
+ if(trkClu[icl].ntrks > BestCluster.ntrks) BestCluster = trkClu[icl];
+ }
+ //store the tracks inside the cluster
+ std::vector<Tracklet> unusedTracks;
+ for(std::vector<Tracklet>::iterator trkItr=tracks.begin(); trkItr!=tracks.end(); ++trkItr) {
+ float dEta = fabs(BestCluster.eta - trkItr->globalPosition().eta());
+ float dPhi = BestCluster.phi - trkItr->globalPosition().phi();
+ if(dPhi > PI) dPhi -= 2*PI;
+ else if(dPhi < -PI) dPhi += 2*PI;
+ if(dEta < 0.7 && fabs(dPhi) < PI/3.) BestCluster.tracks.push_back( (*trkItr) );
+ else unusedTracks.push_back( (*trkItr) );
+ }
+ //return the best cluster and the unused tracklets
+ tracks = unusedTracks;
+ return BestCluster;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<Muon::MSVertexRecoTool::TrkCluster> MSVertexRecoTool::findTrackClusters(std::vector<Tracklet>& tracks) {
+ std::vector<Tracklet> trks = tracks;
+ std::vector<TrkCluster> clusters;
+ //keep making clusters until there are no more possible
+ while(true) {
+ TrkCluster clust = ClusterizeTracks(trks);
+ if(clust.ntrks >= 3) clusters.push_back(clust);
+ else break;
+ if(trks.size() < 3) break;
+ }
+
+ if(clusters.size() == 0) {
+ TrkCluster clust;
+ clusters.push_back(clust);
+ }
+
+ return clusters;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertex* MSVertexRecoTool::MSVxFinder(std::vector<Tracklet>& tracklets) {
+ int nTrkToVertex(0);
+ float NominalAngle(m_TrackPhiAngle),MaxOpenAngle(0.20+m_TrackPhiAngle);
+ float aveX(0),aveY(0),aveR(0),aveZ(0);
+ float maxError(200);//maximum allowed uncertainty for tracklets [mm] (above this cut tracklets are discarded)
+
+ //find the average position of the tracklets
+ for(std::vector<Tracklet>::iterator trkItr=tracklets.begin(); trkItr!=tracklets.end(); ++trkItr) {
+ aveR += trkItr->globalPosition().perp();
+ aveX += trkItr->globalPosition().x();
+ aveY += trkItr->globalPosition().y();
+ aveZ += trkItr->globalPosition().z();
+ }
+ aveX = aveX/(float)tracklets.size();
+ aveY = aveY/(float)tracklets.size();
+ aveZ = aveZ/(float)tracklets.size();
+ aveR = aveR/(float)tracklets.size();
+ float avePhi = atan2(aveY,aveX);
+ while(fabs(avePhi) > PI) {
+ if(avePhi < 0) avePhi += 2*PI;
+ else avePhi -= 2*PI;
+ }
+ //calculate the two angles (theta & phi)
+ float LoF = atan2(aveR,aveZ);//Line of Flight (theta)
+ avePhi = vxPhiFinder(LoF,avePhi);
+
+ //find the positions of the radial planes
+ float Rpos[MAXPLANES];
+ float RadialDist = m_VertexMaxRadialPlane - m_VertexMinRadialPlane;
+ float LoFdist = fabs(RadialDist/sin(LoF));
+ int nplanes = LoFdist/m_VxPlaneDist;
+ float PlaneSpacing = 3500/(nplanes-1.0);
+ for(int k=0; k<nplanes; ++k) Rpos[k] = m_VertexMinRadialPlane + PlaneSpacing*k;
+
+ //loop on tracklets and create two types of track parameters -- nominal and phi shifted tracklets
+ std::vector<const Trk::TrackParameters*> TracksForVertexing[MAXPLANES];//vector of tracklets to be used at each vertex plane
+ std::vector<const Trk::TrackParameters*> TracksForErrors[MAXPLANES];//vector of tracklets to be used for uncertainty at each vertex plane
+ std::vector<bool> isNeutralTrack[MAXPLANES];
+ for(unsigned int i=0; i<tracklets.size(); ++i) {
+ if(tracklets.at(i).mdtChamber() <= 11 || tracklets.at(i).mdtChamber() == 52) {//only barrel tracklets
+ nTrkToVertex++;
+ //coordinate transform variables
+ Amg::Vector3D trkgpos(tracklets.at(i).globalPosition().perp()*cos(avePhi),
+ tracklets.at(i).globalPosition().perp()*sin(avePhi),
+ tracklets.at(i).globalPosition().z());
+ float x0 = trkgpos.x();
+ float y0 = trkgpos.y();
+ float r0 = trkgpos.perp();
+
+ ////////Tracks for computing the error//////
+ //////coordinate transform mess
+ //decide which way the tracklet gets rotated -- positive or negative phi
+ float anglesign = 1.0;
+ if((tracklets.at(i).globalPosition().phi() - avePhi) < 0) anglesign = -1.0;
+ float NominalTrkAng = anglesign*NominalAngle;//in case there is a nominal tracklet angle
+ float MaxTrkAng = anglesign*MaxOpenAngle;//the rotated tracklet phi position
+
+ //loop over the planes
+ for(int k=0; k<nplanes; ++k) {
+ if(Rpos[k] > tracklets.at(i).globalPosition().perp()) break; //only use tracklets that start AFTER the vertex plane
+ //nominal tracks for vertexing
+ float Xp = Rpos[k]*cos(avePhi);
+ float Yp = Rpos[k]*sin(avePhi);
+ //in case there is a nominal opening angle, calculate tracklet direction
+ //the tracklet must cross the candidate vertex plane at the correct phi
+ float DelR = sqrt(sq(x0-Xp)+sq(y0-Yp))/cos(NominalAngle);
+ float X1 = DelR*cos(NominalTrkAng+avePhi) + Xp;
+ float Y1 = DelR*sin(NominalTrkAng+avePhi) + Yp;
+ float R1 = sqrt(sq(X1)+sq(Y1));
+ float Norm = r0/R1;
+ X1 = X1*Norm;
+ Y1 = Y1*Norm;
+ float Dirmag = sqrt(sq(X1-Xp)+sq(Y1-Yp));
+ float Xdir = (X1-Xp)/Dirmag;
+ float Ydir = (Y1-Yp)/Dirmag;
+ float trkpx = Xdir*tracklets.at(i).momentum().perp();
+ float trkpy = Ydir*tracklets.at(i).momentum().perp();
+ float trkpz = tracklets.at(i).momentum().z();
+ //check if the tracklet has a charge & momentum measurement -- if not, set charge=1 so extrapolator will work
+ float charge = tracklets.at(i).charge();
+ if(fabs(charge) < 0.1) {
+ charge = 1; //for "straight" tracks, set charge = 1
+ isNeutralTrack[k].push_back(true);
+ }
+ else isNeutralTrack[k].push_back(false);
+
+ //store the tracklet as a Trk::Perigee
+ Amg::Vector3D trkmomentum(trkpx,trkpy,trkpz);
+ Amg::Vector3D trkgpos(X1,Y1,tracklets.at(i).globalPosition().z());
+ AmgSymMatrix(5) * covariance = new AmgSymMatrix(5)(tracklets.at(i).errorMatrix());
+ Trk::Perigee * myPerigee = new Trk::Perigee(0.,0.,trkmomentum.phi(),trkmomentum.theta(),charge/trkmomentum.mag(),Trk::PerigeeSurface(trkgpos),covariance);
+ TracksForVertexing[k].push_back(myPerigee);
+
+ //tracks for errors -- rotate the plane & recalculate the tracklet parameters
+ float xp = Rpos[k]*cos(avePhi);
+ float yp = Rpos[k]*sin(avePhi);
+ float delR = sqrt(sq(x0-xp)+sq(y0-yp))/cos(MaxOpenAngle);
+ float x1 = delR*cos(MaxTrkAng+avePhi) + xp;
+ float y1 = delR*sin(MaxTrkAng+avePhi) + yp;
+ float r1 = sqrt(sq(x1)+sq(y1));
+ float norm = r0/r1;
+ x1 = x1*norm;
+ y1 = y1*norm;
+ float dirmag = sqrt(sq(x1-xp)+sq(y1-yp));
+ float xdir = (x1-xp)/dirmag;
+ float ydir = (y1-yp)/dirmag;
+ float errpx = xdir*tracklets.at(i).momentum().perp();
+ float errpy = ydir*tracklets.at(i).momentum().perp();
+ float errpz = tracklets.at(i).momentum().z();
+
+ //store the tracklet as a Trk::Perigee
+ AmgSymMatrix(5) * covariance2 = new AmgSymMatrix(5)(tracklets.at(i).errorMatrix());
+ Amg::Vector3D trkerrmom(errpx,errpy,errpz);
+ Amg::Vector3D trkerrpos(x1,y1,tracklets.at(i).globalPosition().z());
+ Trk::Perigee * errPerigee = new Trk::Perigee(0.,0.,trkerrmom.phi(),trkerrmom.theta(),charge/trkerrmom.mag(),Trk::PerigeeSurface(trkerrpos),covariance2);
+ TracksForErrors[k].push_back(errPerigee);
+ }//end loop on vertex planes
+ }//end selection of barrel tracks
+ }//end loop on tracks
+
+ //return if there are not enough tracklets
+ if(nTrkToVertex < 3) return NULL;
+
+ //calculate the tracklet positions on each surface
+ bool boundaryCheck = true;
+ std::vector<float> ExtrapZ[MAXPLANES], dlength[MAXPLANES];//extrapolated position & uncertainty
+ std::vector<std::pair<unsigned int, unsigned int> > UsedTracks[MAXPLANES];
+ std::vector<bool> ExtrapSuc[MAXPLANES];//did the extrapolation succeed?
+ std::vector<MSVertex*> m_vertices;
+ m_vertices.reserve(nplanes);
+
+ std::vector<float> sigmaZ[MAXPLANES];//total uncertainty at each plane
+ std::vector<Amg::Vector3D> pAtVx[MAXPLANES];//tracklet momentum expressed at the plane
+ for(int k=0; k<nplanes; ++k) {
+ float rpos = Rpos[k];
+ for(unsigned int i=0; i<TracksForVertexing[k].size(); ++i) {
+ Amg::Transform3D* surfaceTransformMatrix = new Amg::Transform3D;
+ surfaceTransformMatrix->setIdentity();
+ Trk::CylinderSurface cyl(surfaceTransformMatrix, rpos, 10000.);//create the surface
+ //extrapolate to the surface
+ const Trk::AtaCylinder* extrap = dynamic_cast<const Trk::AtaCylinder*>(m_extrapolator->extrapolate(*TracksForVertexing[k].at(i), cyl,Trk::anyDirection,boundaryCheck,Trk::muon));
+ if(extrap) {
+ //if the track is neutral just store the uncertainty due to angular uncertainty of the orignal tracklet
+ if(isNeutralTrack[k].at(i)) {
+ float pTot = sqrt(sq(TracksForVertexing[k].at(i)->momentum().perp())+sq(TracksForVertexing[k].at(i)->momentum().z()));
+ float dirErr = Amg::error(*TracksForVertexing[k].at(i)->covariance(),Trk::theta);
+ float extrapRdist = TracksForVertexing[k].at(i)->position().perp() - Rpos[k];
+ float sz = fabs(20*dirErr*extrapRdist*sq(pTot)/sq(TracksForVertexing[k].at(i)->momentum().perp()));
+ float ExtrapErr = sz;
+ if(ExtrapErr > maxError) ExtrapSuc[k].push_back(false);
+ else {
+ ExtrapSuc[k].push_back(true);
+ std::pair<unsigned int,unsigned int> trkmap(ExtrapZ[k].size(),i);
+ UsedTracks[k].push_back(trkmap);
+ ExtrapZ[k].push_back( extrap->localPosition().y() );
+ sigmaZ[k].push_back(sz);
+ pAtVx[k].push_back(extrap->momentum());
+ dlength[k].push_back(0);
+ }
+ }//end neutral tracklets
+ //if the tracklet has a momentum measurement
+ else {
+ //now extrapolate taking into account the extra path length & differing magnetic field
+ Amg::Transform3D* srfTransMat2 = new Amg::Transform3D;
+ srfTransMat2->setIdentity();
+ Trk::CylinderSurface cyl2(srfTransMat2, rpos, 10000.);
+ const Trk::AtaCylinder* extrap2 = dynamic_cast<const Trk::AtaCylinder*>(m_extrapolator->extrapolate(*TracksForErrors[k].at(i), cyl2,Trk::anyDirection,boundaryCheck,Trk::muon));
+ if(extrap2) {
+ float sz = Amg::error(*extrap->covariance(),Trk::locY);
+ float zdiff = extrap->localPosition().y() - extrap2->localPosition().y();
+ float ExtrapErr = sqrt(sq(sz)+sq(zdiff));
+ if(ExtrapErr > maxError) ExtrapSuc[k].push_back(false);
+ else {
+ //iff both extrapolations succeed && error is acceptable, store the information
+ ExtrapSuc[k].push_back(true);
+ std::pair<unsigned int,unsigned int> trkmap(ExtrapZ[k].size(),i);
+ UsedTracks[k].push_back(trkmap);
+ ExtrapZ[k].push_back( extrap->localPosition().y() );
+ sigmaZ[k].push_back(sz);
+ pAtVx[k].push_back(extrap->momentum());
+ dlength[k].push_back(zdiff);
+ }
+ }
+ else ExtrapSuc[k].push_back(false);//not possible to calculate the uncertainty -- do not use tracklet in vertex
+ delete extrap2;
+ }
+ }//fi extrap
+ else ExtrapSuc[k].push_back(false);//not possible to extrapolate the tracklet
+ delete extrap;
+ }//end loop on perigeebase
+ }//end loop on radial planes
+ //vertex routine
+ std::vector<Amg::Vector3D> trkp[MAXPLANES];
+ //loop on planes
+ for(int k=0; k<nplanes; ++k) {
+ if(ExtrapZ[k].size() < 3) continue;//require at least 3 tracklets to build a vertex
+ //initialize the variables used in the routine
+ float zLoF = Rpos[k]/tan(LoF);
+ float dzLoF(10);
+ float aveZpos(0),posWeight(0);
+ for(unsigned int i=0; i<ExtrapZ[k].size(); ++i) {
+ float ExtrapErr = sqrt(sq(sigmaZ[k][i])+sq(dlength[k][i])+sq(dzLoF));
+ if(isNeutralTrack[k][i]) ExtrapErr = sqrt(sq(sigmaZ[k][i]) + sq(dzLoF));
+ aveZpos += ExtrapZ[k][i]/sq(ExtrapErr);
+ posWeight += 1./sq(ExtrapErr);
+ }
+ //calculate the weighted average position of the tracklets
+ zLoF = aveZpos/posWeight;
+ float zpossigma(dzLoF),Chi2(0),Chi2Prob(-1);
+ int Nitr(0);
+ std::vector<unsigned int> vxtracks;//tracklets to be used in the vertex routine
+ std::vector<bool> blacklist;//tracklets that do not belong to the vertex
+ for(unsigned int i=0; i<ExtrapZ[k].size(); ++i) blacklist.push_back(false);
+ //minimum chi^2 iterative fit
+ while(true) {
+ vxtracks.clear();
+ trkp[k].clear();
+ int tmpnTrks(0);
+ float tmpzLoF(0),tmpzpossigma(0),tmpchi2(0),posWeight(0);
+ float worstdelz(0);
+ unsigned int iworst(0xC0FFEE);
+ //loop on the tracklets, find the chi^2 contribution from each tracklet
+ for(unsigned int i=0; i<ExtrapZ[k].size(); ++i) {
+ if(blacklist[i]) continue;
+ trkp[k].push_back(pAtVx[k][i]);
+ float delz = zLoF - ExtrapZ[k][i];
+ float ExtrapErr = sqrt(sq(sigmaZ[k][i])+sq(dlength[k][i])+sq(dzLoF));
+ float trkchi2 = sq(delz)/sq(ExtrapErr);
+ if(trkchi2 > worstdelz) {
+ iworst = i;
+ worstdelz = trkchi2;
+ }
+ tmpzLoF += ExtrapZ[k][i]/sq(ExtrapErr);
+ posWeight += 1./sq(ExtrapErr);
+ tmpzpossigma += sq(delz);
+ tmpchi2 += trkchi2;
+ tmpnTrks++;
+ }//end loop on tracks
+ if(tmpnTrks < 3) break;//stop searching for a vertex at this plane
+ tmpzpossigma = sqrt(tmpzpossigma/(float)tmpnTrks);
+ zLoF = tmpzLoF/posWeight;
+ zpossigma = tmpzpossigma;
+ float testChi2 = TMath::Prob(tmpchi2,tmpnTrks-1);
+ if(testChi2 < VxChi2ProbCUT) blacklist[iworst] = true;
+ else {
+ Chi2 = tmpchi2;
+ Chi2Prob = testChi2;
+ //loop on the tracklets and find all that belong to the vertex
+ for(unsigned int i=0; i<ExtrapZ[k].size(); ++i) {
+ float delz = zLoF - ExtrapZ[k][i];
+ float ExtrapErr = sqrt(sq(sigmaZ[k][i])+sq(dlength[k][i])+sq(dzLoF));
+ float trkErr = sqrt(sq(ExtrapErr)+sq(zpossigma)) + 0.001;
+ float trkNsigma = fabs(delz/trkErr);
+ if(trkNsigma < 3) vxtracks.push_back(i);
+ }
+ break;//found a vertex, stop removing tracklets!
+ }
+ if( Nitr >= ((int)ExtrapZ[k].size()-3) ) break;//stop searching for a vertex at this plane
+ Nitr++;
+ }//end while
+ if(vxtracks.size() < 3) continue;
+ std::vector<xAOD::TrackParticle*> vxTrackParticles;
+ //create TrackParticle vector for all tracklets used in the vertex fit
+ vxTrackParticles.reserve(vxtracks.size());
+ for(std::vector<unsigned int>::iterator vxtrk=vxtracks.begin(); vxtrk!=vxtracks.end(); ++vxtrk) {
+ for(unsigned int i=0; i<UsedTracks[k].size(); ++i) {
+ if( (*vxtrk) == UsedTracks[k].at(i).first ) {
+ Tracklet trklt = tracklets.at(UsedTracks[k].at(i).second);
+ AmgSymMatrix(5) * covariance = new AmgSymMatrix(5)(trklt.errorMatrix());
+ Trk::Perigee * myPerigee = new Trk::Perigee(0.,0.,trklt.momentum().phi(),trklt.momentum().theta(),trklt.charge()/trklt.momentum().mag(),
+ Trk::PerigeeSurface(trklt.globalPosition()),covariance);
+ xAOD::TrackParticle* trackparticle = new xAOD::TrackParticle();
+ trackparticle->setFitQuality(1.,(int)trklt.mdtHitsOnTrack().size());
+ trackparticle->setTrackProperties(xAOD::TrackProperties::LowPtTrack);
+
+ trackparticle->setDefiningParameters(myPerigee->parameters()[Trk::d0],
+ myPerigee->parameters()[Trk::z0],
+ myPerigee->parameters()[Trk::phi0],
+ myPerigee->parameters()[Trk::theta],
+ myPerigee->parameters()[Trk::qOverP]);
+ std::vector<float> covMatrixVec;
+ Amg::compress(*covariance,covMatrixVec);
+ trackparticle->setDefiningParametersCovMatrixVec(covMatrixVec);
+
+ vxTrackParticles.push_back(trackparticle);
+
+ delete myPerigee;
+ break;
+ }
+ }
+ }
+ Amg::Vector3D m_position(Rpos[k]*cos(avePhi),Rpos[k]*sin(avePhi),zLoF);
+ m_vertices.push_back( new MSVertex(1,m_position,vxTrackParticles,Chi2Prob,Chi2,0,0,0) );
+ }//end loop on Radial planes
+
+ //delete the perigeebase
+ for(int k=0; k<nplanes; ++k) {
+ for(unsigned int i=0; i<TracksForVertexing[k].size(); ++i) delete TracksForVertexing[k].at(i);
+ for(unsigned int i=0; i<TracksForErrors[k].size(); ++i) delete TracksForErrors[k].at(i);
+ }
+
+ //return an empty vertex in case none were reconstructed
+ if(m_vertices.size() == 0) {
+ return NULL;
+ }
+
+ //loop on the vertex candidates and select the best based on max n(tracks) and max chi^2 probability
+ unsigned int bestVx(0);
+ for(unsigned int k=1; k<m_vertices.size(); ++k) {
+ if(m_vertices[k]->getChi2Probability() < VxChi2ProbCUT || m_vertices[k]->getNTracks() < 3) continue;
+ if(m_vertices[k]->getNTracks() < m_vertices[bestVx]->getNTracks()) continue;
+ if(m_vertices[k]->getNTracks() == m_vertices[bestVx]->getNTracks() && m_vertices[k]->getChi2Probability() < m_vertices[bestVx]->getChi2Probability()) continue;
+ bestVx = k;
+ }
+ MSVertex* myVertex = m_vertices[bestVx]->clone();
+ //cleanup
+ for(std::vector<MSVertex*>::iterator it=m_vertices.begin(); it!=m_vertices.end(); ++it) {
+ delete (*it);
+ (*it) = 0;
+ }
+ m_vertices.clear();
+
+ return myVertex;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertex* MSVertexRecoTool::MSStraightLineVx(std::vector<Tracklet> trks) {
+ // Running set of all vertices found. The inner set is the indices of trks that are used to make the vertex
+ std::set<std::set<int> > prelim_vx;
+
+ // Make a list of all 3-tracklet combinations that make vertices
+ for(unsigned int i = 0; i < trks.size() - 2; i++) {
+ for(unsigned int j = i+1; j < trks.size() - 1; j++) {
+ for(unsigned int k = j+1; k < trks.size(); k++) {
+ std::set<int> tmpTracks;
+ tmpTracks.insert(i);
+ tmpTracks.insert(j);
+ tmpTracks.insert(k);
+
+ Amg::Vector3D MyVx;
+ MyVx = VxMinQuad(getTracklets(trks, tmpTracks));
+ if(MyVx.perp() < 10000 && fabs(MyVx.z()) > 7000 && fabs(MyVx.z()) < 15000 && !EndcapHasBadTrack(getTracklets(trks, tmpTracks), MyVx))
+ prelim_vx.insert(tmpTracks);
+ }
+ }
+ }
+
+ // If no preliminary vertices were found from 3 tracklets, then there is no vertex and we are done.
+ if(prelim_vx.size() == 0)
+ return NULL;
+
+ // The remaining algorithm is very time consuming for large numbers of tracklets. To control this,
+ // we run the old algorithm when there are too many tracklets and a vertex is found.
+ if(trks.size() <= 20) {
+ std::set<std::set<int> > new_prelim_vx = prelim_vx;
+ std::set<std::set<int> > old_prelim_vx;
+
+ int foundNewVx = true;
+ while(foundNewVx) {
+ foundNewVx = false;
+
+ old_prelim_vx = new_prelim_vx;
+ new_prelim_vx.clear();
+
+ for(std::set<std::set<int> >::iterator itr = old_prelim_vx.begin(); itr != old_prelim_vx.end(); itr++) {
+ for(unsigned int i_trks = 0; i_trks < trks.size(); i_trks++) {
+ std::set<int> tempCluster = *itr;
+ if(tempCluster.insert(i_trks).second) {
+ Amg::Vector3D MyVx = VxMinQuad(getTracklets(trks, tempCluster));
+ if(MyVx.perp() < 10000 && fabs(MyVx.z()) > 7000 && fabs(MyVx.z()) < 15000 && !EndcapHasBadTrack(getTracklets(trks, tempCluster), MyVx)) {
+ new_prelim_vx.insert(tempCluster);
+ prelim_vx.insert(tempCluster);
+ foundNewVx = true;
+ }
+ }
+ }
+ }
+ }
+ }
+ else {
+ // Since there are 20 or more tracklets, we're going to use the old MSVx finding method. Note that
+ // if the old method fails, we do not return here; in this case a 3-tracklet vertex that was found
+ // earlier in this algorithm will be returned
+ MSVertex* tempVertex = MSStraightLineVx_oldMethod(trks);
+ if(tempVertex)
+ return tempVertex;
+ }
+
+ // Find the preliminary vertex with the maximum number of tracklets - that is the final vertex. If
+ // multiple preliminary vertices with same number of tracklets, the first one found is returned
+ std::set<std::set<int> >::iterator prelim_vx_max = prelim_vx.begin();
+ for(std::set<std::set<int> >::iterator itr = prelim_vx.begin(); itr != prelim_vx.end(); itr++) {
+ if((*itr).size() > (*prelim_vx_max).size())
+ prelim_vx_max = itr;
+ }
+
+ std::vector<Tracklet> tracklets = getTracklets(trks, *prelim_vx_max);
+ // use tracklets to estimate the line of flight of decaying particle
+ float aveX(0),aveY(0);
+ for(std::vector<Tracklet>::iterator trkItr=tracklets.begin(); trkItr!=tracklets.end(); ++trkItr) {
+ aveX += ((Tracklet)*trkItr).globalPosition().x();
+ aveY += ((Tracklet)*trkItr).globalPosition().y();
+ }
+ float tracklet_vxphi = atan2(aveY,aveX);
+ Amg::Vector3D MyVx = VxMinQuad(tracklets);
+ float vxtheta = atan2(MyVx.x(),MyVx.z());
+ float vxphi = vxPhiFinder(vxtheta,tracklet_vxphi);
+
+ Amg::Vector3D vxpos(MyVx.x()*cos(vxphi),MyVx.x()*sin(vxphi),MyVx.z());
+ std::vector<xAOD::TrackParticle*> vxTrkTracks;
+ for(std::vector<Tracklet>::iterator tracklet = tracklets.begin(); tracklet != tracklets.end(); tracklet++) {
+ AmgSymMatrix(5)* covariance = new AmgSymMatrix(5)(((Tracklet)*tracklet).errorMatrix());
+ Trk::Perigee* myPerigee = new Trk::Perigee(vxpos,((Tracklet)*tracklet).momentum(),0,vxpos,covariance);
+ const Trk::TrackParameters * trkpar = dynamic_cast<const Trk::TrackParameters*>(myPerigee);
+ DataVector<const Trk::TrackStateOnSurface>* trackStateOnSurfaces = new DataVector<const Trk::TrackStateOnSurface>;
+ trackStateOnSurfaces->push_back( new Trk::TrackStateOnSurface(0,trkpar) );
+ Trk::TrackInfo info( Trk::TrackInfo::Unknown, Trk::undefined);
+
+ xAOD::TrackParticle* myTrack = new xAOD::TrackParticle();
+ myTrack->setFitQuality(1.,(int)((Tracklet)*tracklet).mdtHitsOnTrack().size());
+ myTrack->setDefiningParameters(myPerigee->parameters()[Trk::d0],
+ myPerigee->parameters()[Trk::z0],
+ myPerigee->parameters()[Trk::phi0],
+ myPerigee->parameters()[Trk::theta],
+ myPerigee->parameters()[Trk::qOverP]);
+
+ vxTrkTracks.push_back(myTrack);
+ }
+
+ MSVertex* vertex = new MSVertex(2,vxpos,vxTrkTracks,1,0,0,0,0);
+ return vertex;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ //vertex finding routine for the endcap
+ MSVertex* MSVertexRecoTool::MSStraightLineVx_oldMethod(std::vector<Tracklet> trks) {
+ //find the line of flight of the vpion
+ float aveX(0),aveY(0),aveR(0),aveZ(0);
+ for(std::vector<Tracklet>::iterator trkItr=trks.begin(); trkItr!=trks.end(); ++trkItr) {
+ aveX += trkItr->globalPosition().x();
+ aveY += trkItr->globalPosition().y();
+ aveR += trkItr->globalPosition().perp();
+ aveZ += trkItr->globalPosition().z();
+ }
+ float vxphi = atan2(aveY,aveX);
+
+ Amg::Vector3D MyVx(0,0,0);
+ std::vector<Tracklet> tracks = RemoveBadTrk(trks,MyVx);
+ if(tracks.size() < 2) return NULL;
+
+ while(true) {
+ MyVx = VxMinQuad(tracks);
+ std::vector<Tracklet> Tracks = RemoveBadTrk(tracks,MyVx);
+ if(tracks.size() == Tracks.size()) break;
+ tracks = Tracks;
+ }
+ if(tracks.size() >= 3 && MyVx.x() > 0) {
+ float vxtheta = atan2(MyVx.x(),MyVx.z());
+ vxphi = vxPhiFinder(vxtheta,vxphi);
+ Amg::Vector3D vxpos(MyVx.x()*cos(vxphi),MyVx.x()*sin(vxphi),MyVx.z());
+ //make Trk::Track for each tracklet used in the vertex fit
+ std::vector<xAOD::TrackParticle*> vxTrackParticles;
+ for(std::vector<Tracklet>::iterator trklt=tracks.begin(); trklt!=tracks.end(); ++trklt) {
+ AmgSymMatrix(5) * covariance = new AmgSymMatrix(5)(trklt->errorMatrix());
+ Trk::Perigee* myPerigee = new Trk::Perigee(0.,0.,trklt->momentum().phi(),trklt->momentum().theta(),trklt->charge()/trklt->momentum().mag(),Trk::PerigeeSurface(trklt->globalPosition()),covariance);
+ xAOD::TrackParticle* trackparticle = new xAOD::TrackParticle();
+ trackparticle->setFitQuality(1.,(int)trklt->mdtHitsOnTrack().size());
+ trackparticle->setTrackProperties(xAOD::TrackProperties::LowPtTrack);
+
+ trackparticle->setDefiningParameters(myPerigee->parameters()[Trk::d0],
+ myPerigee->parameters()[Trk::z0],
+ myPerigee->parameters()[Trk::phi0],
+ myPerigee->parameters()[Trk::theta],
+ myPerigee->parameters()[Trk::qOverP]);
+ std::vector<float> covMatrixVec;
+ Amg::compress(*covariance,covMatrixVec);
+ trackparticle->setDefiningParametersCovMatrixVec(covMatrixVec);
+
+ vxTrackParticles.push_back(trackparticle);
+
+ delete myPerigee;
+
+ }
+
+ MSVertex* vertex = new MSVertex(2,vxpos,vxTrackParticles,1,(float)vxTrackParticles.size(),0,0,0);
+ return vertex;
+ }
+ return NULL;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<Tracklet> MSVertexRecoTool::RemoveBadTrk(std::vector<Tracklet> tracks, Amg::Vector3D Vx) {
+ float MaxTollDist = 300;//max distance between the vertex and tracklet [mm]
+ std::vector<Tracklet> Tracks;
+ if(Vx.x() == 0 && Vx.z() == 0) return tracks;
+ //loop on all tracks and find the worst
+ float WorstTrkDist = MaxTollDist;
+ unsigned int iWorstTrk = 0xC0FFEE;
+ for(unsigned int i=0; i<tracks.size(); ++i) {
+ float TrkSlope = tan(tracks.at(i).getML1seg().alpha());
+ float TrkInter = tracks.at(i).getML1seg().globalPosition().perp() - tracks.at(i).getML1seg().globalPosition().z()*TrkSlope;
+ float dist = fabs((TrkSlope*Vx.z() - Vx.x() + TrkInter)/sqrt(sq(TrkSlope) + 1));
+ if(dist > MaxTollDist && dist > WorstTrkDist) {
+ iWorstTrk = i;
+ WorstTrkDist = dist;
+ }
+ }
+
+ //Remove the worst track from the list
+ for(unsigned int i=0; i<tracks.size(); ++i) {
+ if(i != iWorstTrk) Tracks.push_back(tracks.at(i));
+ }
+ return Tracks;
+ }
+
+ std::vector<Tracklet> MSVertexRecoTool::getTracklets(std::vector<Tracklet> trks, std::set<int> tracklet_subset) {
+ std::vector<Tracklet> returnVal;
+ for(std::set<int>::iterator itr = tracklet_subset.begin(); itr != tracklet_subset.end(); itr++) {
+ if((unsigned int) *itr > trks.size())
+ ATH_MSG_ERROR( "ERROR - Index out of bounds in getTracklets");
+ returnVal.push_back(trks.at(*itr));
+ }
+
+ return returnVal;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ bool MSVertexRecoTool::EndcapHasBadTrack(std::vector<Tracklet> tracks, Amg::Vector3D Vx) {
+ float MaxTollDist = m_MaxTollDist;
+ if(Vx.x() == 0 && Vx.z() == 0) return true;
+ //loop on all tracks and find the worst
+ float WorstTrkDist = MaxTollDist;
+ for(std::vector<Tracklet>::iterator track = tracks.begin(); track != tracks.end(); track++) {
+ float TrkSlope = tan(((Tracklet)*track).getML1seg().alpha());
+ float TrkInter = ((Tracklet)*track).getML1seg().globalPosition().perp() - ((Tracklet)*track).getML1seg().globalPosition().z()*TrkSlope;
+ float dist = fabs((TrkSlope*Vx.z() - Vx.x() + TrkInter)/sqrt(sq(TrkSlope) + 1));
+ if(dist > MaxTollDist && dist > WorstTrkDist) {
+ return true;
+ }
+ }
+
+ // No tracks found that are too far, so it is okay.
+ return false;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ //core algorithm for endcap vertex reconstruction
+ Amg::Vector3D MSVertexRecoTool::VxMinQuad(std::vector<Tracklet> tracks) {
+ float s(0.),sx(0.),sy(0.),sxy(0.),sxx(0.),d(0.);
+ float sigma = 1.;
+ for(unsigned int i=0; i<tracks.size(); ++i) {
+ float TrkSlope = tan(tracks.at(i).getML1seg().alpha());
+ float TrkInter = tracks.at(i).getML1seg().globalPosition().perp() - tracks.at(i).getML1seg().globalPosition().z()*TrkSlope;
+ s += 1./(sq(sigma));
+ sx += TrkSlope/(sq(sigma));
+ sxx += sq(TrkSlope)/sq(sigma);
+ sy += TrkInter/sq(sigma);
+ sxy += (TrkSlope*TrkInter)/sq(sigma);
+ }
+ d = s*sxx - sq(sx);
+ float Rpos = (sxx*sy - sx*sxy)/d;
+ float Zpos = (sx*sy - s*sxy)/d;
+
+ Amg::Vector3D MyVx(Rpos,0,Zpos);
+
+ return MyVx;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ //vertex phi location -- determined from the RPC/TGC hits
+ float MSVertexRecoTool::vxPhiFinder(float theta,float phi) {
+ float nmeas(0);
+ float sinphi(0),cosphi(0);
+ float eta = -1*log(tan(0.5*theta));
+ if(fabs(eta) < 1.5) {
+ const Muon::RpcPrepDataContainer* rpcTES = 0;
+ if(evtStore()->retrieve(rpcTES,"RPC_Measurements").isFailure() || !rpcTES) {
+ ATH_MSG_WARNING( "No RpcPrepDataContainer found in SG!" );
+ return 0;
+ }
+ Muon::RpcPrepDataContainer::const_iterator RpcItr = rpcTES->begin();
+ Muon::RpcPrepDataContainer::const_iterator RpcItrE = rpcTES->end();
+ for(; RpcItr != RpcItrE; ++RpcItr) {
+ Muon::RpcPrepDataCollection::const_iterator rpcItr = (*RpcItr)->begin();
+ Muon::RpcPrepDataCollection::const_iterator rpcItrE = (*RpcItr)->end();
+ for(; rpcItr != rpcItrE; ++rpcItr) {
+ if(m_rpcIdHelper->measuresPhi((*rpcItr)->identify())) {
+ float rpcEta = (*rpcItr)->globalPosition().eta();
+ float rpcPhi = (*rpcItr)->globalPosition().phi();
+ float dphi = phi - rpcPhi;
+ if(dphi > PI) dphi -= 2*PI;
+ else if(dphi < -PI) dphi += 2*PI;
+ float deta = eta - rpcEta;
+ float DR = sqrt(sq(deta)+sq(dphi));
+ if(DR < 0.6) {
+ nmeas++;
+ sinphi += sin(rpcPhi);
+ cosphi += cos(rpcPhi);
+ }
+ }
+ }
+ }
+ }
+ if(fabs(eta) > 0.5) {
+ const Muon::TgcPrepDataContainer* tgcTES = 0;
+ if(evtStore()->retrieve(tgcTES,"TGC_Measurements").isFailure()) {
+ ATH_MSG_WARNING( "No TgcPrepDataContainer found in SG!" );
+ return 0;
+ }
+ Muon::TgcPrepDataContainer::const_iterator TgcItr = tgcTES->begin();
+ Muon::TgcPrepDataContainer::const_iterator TgcItrE = tgcTES->end();
+ for(; TgcItr != TgcItrE; ++TgcItr) {
+ Muon::TgcPrepDataCollection::const_iterator tgcItr = (*TgcItr)->begin();
+ Muon::TgcPrepDataCollection::const_iterator tgcItrE = (*TgcItr)->end();
+ for(; tgcItr != tgcItrE; ++tgcItr) {
+ if(m_tgcIdHelper->isStrip((*tgcItr)->identify())) {
+ float tgcEta = (*tgcItr)->globalPosition().eta();
+ float tgcPhi = (*tgcItr)->globalPosition().phi();
+ float dphi = phi - tgcPhi;
+ if(dphi > PI) dphi -= 2*PI;
+ else if(dphi < -PI) dphi += 2*PI;
+ float deta = eta - tgcEta;
+ float DR = sqrt(sq(deta)+sq(dphi));
+ if(DR < 0.6) {
+ nmeas++;
+ sinphi += sin(tgcPhi);
+ cosphi += cos(tgcPhi);
+ }
+ }
+ }
+ }
+ }
+ float vxphi = phi;
+ if(nmeas > 0) vxphi = atan2(sinphi/nmeas,cosphi/nmeas);
+ return vxphi;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ //count the hits (MDT, RPC & TGC) around the vertex
+ void MSVertexRecoTool::HitCounter(MSVertex* MSRecoVx) {
+ int nHighOccupancy(0);
+ const Muon::MdtPrepDataContainer* mdtTES = 0;
+ if(evtStore()->retrieve(mdtTES,"MDT_DriftCircles").isFailure()) {
+ ATH_MSG_ERROR( "Unable to retrieve the MDT hits" );
+ }
+ const Muon::RpcPrepDataContainer* rpcTES = 0;
+ if(evtStore()->retrieve(rpcTES,"RPC_Measurements").isFailure()) {
+ ATH_MSG_ERROR( "Unable to retrieve the RPC hits" );
+ }
+ const Muon::TgcPrepDataContainer* tgcTES = 0;
+ if(evtStore()->retrieve(tgcTES,"TGC_Measurements").isFailure()) {
+ ATH_MSG_ERROR( "Unable to retrieve the MDT hits" );
+ }
+
+ //MDTs -- count the number around the vertex
+ int nmdt(0);
+ Muon::MdtPrepDataContainer::const_iterator MDTItr = mdtTES->begin();
+ Muon::MdtPrepDataContainer::const_iterator MDTItrE = mdtTES->end();
+ for(; MDTItr != MDTItrE; ++MDTItr) {
+ if( (*MDTItr)->size() == 0) continue;
+ Muon::MdtPrepDataCollection::const_iterator mdt = (*MDTItr)->begin();
+ Muon::MdtPrepDataCollection::const_iterator mdtE = (*MDTItr)->end();
+ Amg::Vector3D ChamberCenter = (*mdt)->detectorElement()->center();
+ float deta = fabs(MSRecoVx->getPosition().eta() - ChamberCenter.eta());
+ if(deta > 0.6) continue;
+ float dphi = MSRecoVx->getPosition().phi() - ChamberCenter.phi();
+ if(dphi > PI) dphi -= 2*PI;
+ else if(dphi < -PI) dphi += 2*PI;
+ if( fabs(dphi) > 0.6 ) continue;
+ int nChHits(0);
+ Identifier id = (*mdt)->identify();
+ float nTubes = (m_mdtIdHelper->tubeLayerMax(id) - m_mdtIdHelper->tubeLayerMin(id) + 1)*(m_mdtIdHelper->tubeMax(id) - m_mdtIdHelper->tubeMin(id) + 1);
+ for(; mdt != mdtE; ++mdt) {
+ if((*mdt)->adc() < 50) continue;
+ if((*mdt)->status() != 1) continue;
+ if((*mdt)->localPosition()[Trk::locR] == 0.) continue;
+ nChHits++;
+ }
+ nmdt += nChHits;
+ double ChamberOccupancy = nChHits/nTubes;
+ if(ChamberOccupancy > m_ChamberOccupancyMin) nHighOccupancy++;
+ }
+
+ ATH_MSG_DEBUG( "Found " << nHighOccupancy << " chambers near the MS vertex with occupancy greater than "
+ << m_ChamberOccupancyMin );
+ if(nHighOccupancy < m_minHighOccupancyChambers) return;
+
+ //RPC -- count the number around the vertex
+ int nrpc(0);
+ Muon::RpcPrepDataContainer::const_iterator RpcItr = rpcTES->begin();
+ Muon::RpcPrepDataContainer::const_iterator RpcItrE = rpcTES->end();
+ for(; RpcItr != RpcItrE; ++RpcItr) {
+ Muon::RpcPrepDataCollection::const_iterator rpcItr = (*RpcItr)->begin();
+ Muon::RpcPrepDataCollection::const_iterator rpcItrE = (*RpcItr)->end();
+ for(; rpcItr != rpcItrE; ++rpcItr) {
+ float rpcEta = (*rpcItr)->globalPosition().eta();
+ float rpcPhi = (*rpcItr)->globalPosition().phi();
+ float dphi = MSRecoVx->getPosition().phi() - rpcPhi;
+ if(dphi > PI) dphi -= 2*PI;
+ else if(dphi < -PI) dphi += 2*PI;
+ float deta = MSRecoVx->getPosition().eta() - rpcEta;
+ float DR = sqrt(sq(deta)+sq(dphi));
+ if(DR < 0.6) nrpc++;
+ if(DR > 1.2) break;
+ }
+ }
+ //TGC -- count the number around the vertex
+ int ntgc(0);
+ Muon::TgcPrepDataContainer::const_iterator TgcItr = tgcTES->begin();
+ Muon::TgcPrepDataContainer::const_iterator TgcItrE = tgcTES->end();
+ for(; TgcItr != TgcItrE; ++TgcItr) {
+ Muon::TgcPrepDataCollection::const_iterator tgcItr = (*TgcItr)->begin();
+ Muon::TgcPrepDataCollection::const_iterator tgcItrE = (*TgcItr)->end();
+ for(; tgcItr != tgcItrE; ++tgcItr) {
+ float tgcEta = (*tgcItr)->globalPosition().eta();
+ float tgcPhi = (*tgcItr)->globalPosition().phi();
+ float dphi = MSRecoVx->getPosition().phi() - tgcPhi;
+ if(dphi > PI) dphi -= 2*PI;
+ else if(dphi < -PI) dphi += 2*PI;
+ float deta = MSRecoVx->getPosition().eta() - tgcEta;
+ float DR = sqrt(sq(deta)+sq(dphi));
+ if(DR < 0.6) ntgc++;
+ if(DR > 1.2) break;
+ }
+ }
+
+ MSRecoVx->setNMDT(nmdt);
+ MSRecoVx->setNRPC(nrpc);
+ MSRecoVx->setNTGC(ntgc);
+
+ return;
+ }
+
+}//end namespace
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.h b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.h
new file mode 100755
index 0000000000000..72e4519717a9c
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexRecoTool.h
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MSVVERTEXRECOTOOL_H
+#define MSVVERTEXRECOTOOL_H
+
+#include "GaudiKernel/AlgTool.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "EventPrimitives/EventPrimitives.h"
+#include "MuonReadoutGeometry/MdtReadoutElement.h"
+#include "GeoAdaptors/GeoMuonHits.h"
+#include "Identifier/Identifier.h"
+#include "TrkExInterfaces/IExtrapolator.h"
+#include "TrkParameters/TrackParameters.h"
+#include "MSVertexToolInterfaces/IMSVertexRecoTool.h"
+#include "MSVertexUtils/Tracklet.h"
+#include "MSVertexUtils/MSVertex.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+
+#include <utility>
+#include <vector>
+
+class StoreGate;
+class MsgStream;
+class Identifier;
+
+namespace Muon {
+
+ class MSVertexRecoTool : virtual public IMSVertexRecoTool, public AthAlgTool
+ {
+
+ public :
+ /** default constructor */
+ MSVertexRecoTool (const std::string& type, const std::string& name, const IInterface* parent);
+ /** destructor */
+ virtual ~MSVertexRecoTool();
+
+ virtual StatusCode initialize(void);
+ virtual StatusCode finalize(void);
+
+ struct TrkCluster {
+ float eta;
+ float phi;
+ int ntrks;
+ int trks[100];
+ bool isSystematic;
+ std::vector<Tracklet> tracks;
+ };
+
+
+ private:
+ //add tool handles, private variables, etc here
+ ToolHandle <Trk::IExtrapolator> m_extrapolator;
+ std::string xAODContainerName;
+ const MdtIdHelper* m_mdtIdHelper;
+ const RpcIdHelper* m_rpcIdHelper;
+ const TgcIdHelper* m_tgcIdHelper;
+ float m_BarrelTrackletUncert;
+ float m_EndcapTrackletUncert;
+ float m_TrackPhiAngle;
+ float VxChi2ProbCUT;
+ float m_VxPlaneDist;
+ float m_VertexMaxRadialPlane;
+ float m_VertexMinRadialPlane;
+ int m_minHighOccupancyChambers;
+ int m_ChamberOccupancyMin;
+ int m_useOldMSVxEndcapMethod;
+ float m_MaxTollDist;
+ float PI;
+ bool m_doSystematics;
+
+ CLHEP::HepRandomEngine* m_rndmEngine;
+ std::string m_rndmEngineName;
+ IAtRndmGenSvc* m_rndmSvc;
+
+
+ public:
+
+ StatusCode findMSvertices(std::vector<Tracklet>& traklets);
+
+ private:
+ MSVertex* MSVxFinder(std::vector<Tracklet>& tracklets);//barrel vertex reco algorithm
+ MSVertex* MSStraightLineVx(std::vector<Tracklet> trks);//endcap vertex reco algorithm
+ MSVertex* MSStraightLineVx_oldMethod(std::vector<Tracklet> trks);
+ float vxPhiFinder(float theta,float phi);//vertex phi location reco algorithm
+ void HitCounter(MSVertex* MSRecoVx);//counts MDT, RPC & TGC around a reco'd vertex
+ std::vector<TrkCluster> findTrackClusters(std::vector<Tracklet>& tracklets);//group tracklets into clusters -- vertex reco runs on each cluster of tracklets
+ TrkCluster ClusterizeTracks(std::vector<Tracklet>& tracks);//core algorithm for creating the clusters
+ Amg::Vector3D VxMinQuad(std::vector<Tracklet> tracks);//endcap vertex reco core
+ std::vector<Tracklet> RemoveBadTrk(std::vector<Tracklet> tracklets, Amg::Vector3D Vx);//endcap vertex algo track selector
+ bool EndcapHasBadTrack(std::vector<Tracklet> tracklets, Amg::Vector3D Vx);
+ std::vector<Tracklet> getTracklets(std::vector<Tracklet> trks, std::set<int> tracklet_subset);
+ float sq(float x) { return (x)*(x); }
+
+ };
+
+
+}
+#endif //MSVERTEXRECOTOOL_H
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.cxx b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.cxx
new file mode 100755
index 0000000000000..15b4109456213
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.cxx
@@ -0,0 +1,1295 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MSVertexTrackletTool.h"
+#include "DataModel/DataVector.h"
+#include "EventPrimitives/EventPrimitivesHelpers.h"
+#include "TrkParameters/TrackParameters.h"
+#include "xAODTracking/TrackParticle.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/TrackParticleAuxContainer.h"
+
+/*
+ Tracklet reconstruction tool
+ See documentation at https://cds.cern.ch/record/1455664 and https://cds.cern.ch/record/1520894
+
+ station name reference:
+ Barrel | Endcap
+ 0 == BIL 6 == BEE | 13 == EIL 21 == EOS
+ 1 == BIS 7 == BIR | 14 == EEL 49 == EIS
+ 2 == BML 8 == BMF | 15 == EES
+ 3 == BMS 9 == BOF | 17 == EML
+ 4 == BOL 10 == BOG | 18 == EMS
+ 5 == BOS 52 == BIM | 20 == EOL
+
+*/
+
+
+namespace Muon {
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ const MdtIdHelper* MSVertexTrackletTool::mdtCompareIdHelper;
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertexTrackletTool::MSVertexTrackletTool (const std::string& type, const std::string& name,
+ const IInterface* parent)
+ :
+ AthAlgTool(type, name, parent) {
+ declareInterface<IMSVertexTrackletTool>(this);
+
+ declareProperty("xAODTrackParticleContainer",m_TPContainer = "MSonlyTracklets");
+ // max residual for tracklet seeds
+ declareProperty("SeedResidual",m_SeedResidual = 5);
+ // segment fitter chi^2 cut
+ declareProperty("MinSegFinderChi2Prob",m_minSegFinderChi2 = 0.05);
+ // maximum delta_alpha allowed in barrel MS chambers
+ declareProperty("BarrelDeltaAlpha",m_BarrelDeltaAlphaCut = 0.100);
+ // maximum delta_b allowed
+ declareProperty("maxDeltab",m_maxDeltabCut = 3);
+ // maximum delta_alpha allowed in the endcap MS chambers
+ declareProperty("EndcapDeltaAlpha",m_EndcapDeltaAlphaCut = 0.015);
+ // tight tracklet requirement (affects efficiency - disabled by default)
+ declareProperty("TightTrackletRequirement",m_tightTrackletRequirement = false);
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexTrackletTool::initialize() {
+
+ if( AthAlgTool::initialize().isFailure() ) {
+ ATH_MSG_ERROR("Failed to initialize AthAlgTool " );
+ return StatusCode::FAILURE;
+ }
+
+ if(detStore()->retrieve(m_mdtIdHelper,"MDTIDHELPER").isFailure()) {
+ ATH_MSG_ERROR("Failed to retrieve the mdtIdHelper");
+ return StatusCode::FAILURE;
+ }
+ mdtCompareIdHelper = m_mdtIdHelper;
+
+ //CONSTANTS
+ PI = 3.1415927;
+ //Delta Alpha Constants -- p = k/(delta_alpha)
+ k_BIL = 28.4366;//MeV*mrad
+ k_BML = 62.8267;//MeV*mrad
+ k_BMS = 53.1259;//MeV*mrad
+ k_BOL = 29.7554;//MeV*mrad
+
+ return StatusCode::SUCCESS;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ MSVertexTrackletTool::~MSVertexTrackletTool() {
+
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexTrackletTool::finalize() {
+ return StatusCode::SUCCESS;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ StatusCode MSVertexTrackletTool::findTracklets(std::vector<Tracklet>& tracklets) {
+
+ //sort the MDT hits into chambers & MLs
+ std::vector<std::vector<Muon::MdtPrepData*> > SortedMdt;
+ nMDT = SortMDThits(SortedMdt);
+
+ //loop over the MDT hits and find segments
+ //select the tube combinations to be fit
+ /*Select hits in at least 2 layers and require hits be ordered by increasing tube number (see diagrams below).
+ ( )( )(3)( ) ( )(3)( )( ) ( )( )( )( ) ( )( )(3)( ) ( )(2)(3)( )
+ ( )(2)( )( ) ( )(2)( )( ) ( )(2)(3)( ) ( )(1)(2)( ) ( )(1)( )( )
+ (1)( )( )( ) (1)( )( )( ) (1)( )( )( ) ( )( )( )( ) ( )( )( )( )
+ Barrel selection criteria: |z_mdt1 - z_mdt2| < 100 mm, |z_mdt1 - z_mdt3| < 160 mm
+ Endcap selection criteria: |r_mdt1 - r_mdt2| < 100 mm, |r_mdt1 - r_mdt3| < 160 mm
+ */
+
+ std::vector<TrackletSegment> segs[6][2][16];//single ML segment array (indicies [station type][ML][sector])
+ std::vector<std::vector<Muon::MdtPrepData*> >::const_iterator ChamberItr = SortedMdt.begin();
+ for(; ChamberItr != SortedMdt.end(); ++ChamberItr) {
+ std::vector<TrackletSegment> mlsegments;
+ //loop on hits inside the chamber
+ std::vector<Muon::MdtPrepData*>::const_iterator mdt1 = ChamberItr->begin();
+ std::vector<Muon::MdtPrepData*>::const_iterator mdtEnd = ChamberItr->end();
+ int stName = m_mdtIdHelper->stationName((*mdt1)->identify());
+ int stEta = m_mdtIdHelper->stationEta((*mdt1)->identify());
+ if(stName == 6 || stName == 14 || stName == 15) continue; //ignore hits from BEE, EEL and EES
+ if(stName == 1 && fabs(stEta) >= 7) continue;//ignore hits from BIS7/8
+ //convert to the hardware sector [1-16]
+ int sector = 2*(m_mdtIdHelper->stationPhi((*mdt1)->identify()));
+ if(stName == 0 || stName == 2 || stName == 4 || stName == 13 || stName == 17 || stName == 20) sector -= 1;
+ //information about the chamber we are looking at
+ int maxLayer = m_mdtIdHelper->tubeLayerMax((*mdt1)->identify());
+ int ML = m_mdtIdHelper->multilayer((*mdt1)->identify());
+ for(; mdt1 != mdtEnd; ++mdt1) {
+ int tl1 = m_mdtIdHelper->tubeLayer((*mdt1)->identify());
+ if(tl1 == maxLayer) break;//require hits in at least 2 layers
+ std::vector<Muon::MdtPrepData*>::const_iterator mdt2 = (mdt1+1);
+ for(; mdt2 != mdtEnd; ++mdt2) {
+ //select the correct combinations
+ int tl2 = m_mdtIdHelper->tubeLayer((*mdt2)->identify());
+ if(mdt1 == mdt2 || (tl2 - tl1) > 1 || (tl2 - tl1) < 0 ) continue;
+ if(fabs((*mdt2)->globalPosition().z() - (*mdt1)->globalPosition().z()) > 100 && (stName <= 11 || stName == 52)) continue;
+ //if chamber is endcap, use distance in r
+ if( (stName > 11 && stName <=21) || stName == 49 ) {
+ float mdt1R = (*mdt1)->globalPosition().perp();
+ float mdt2R = (*mdt2)->globalPosition().perp();
+ if(fabs(mdt1R-mdt2R) > 100) continue;
+ }
+ if( (tl2-tl1) == 0 && (m_mdtIdHelper->tube((*mdt2)->identify()) - m_mdtIdHelper->tube((*mdt1)->identify())) < 0) continue;
+ //find the third hit
+ std::vector<Muon::MdtPrepData*>::const_iterator mdt3 = (mdt2+1);
+ for(; mdt3 != mdtEnd; ++mdt3) {
+ //reject the bad tube combinations
+ int tl3 = m_mdtIdHelper->tubeLayer((*mdt3)->identify());
+ if(mdt1 == mdt3 || mdt2 == mdt3) continue;
+ if( (tl3-tl2) > 1 || (tl3-tl2) < 0 || (tl3-tl1) <= 0) continue;
+ if( (tl3-tl2) == 0 && (m_mdtIdHelper->tube((*mdt3)->identify()) - m_mdtIdHelper->tube((*mdt2)->identify())) < 0) continue;
+ if( fabs((*mdt3)->globalPosition().z() - (*mdt1)->globalPosition().z()) > 160 && (stName <= 11 || stName == 52) ) continue;
+ //if chamber is endcap, use distance in r
+ if( (stName > 11 && stName <=21) || stName == 49 ) {
+ float mdt1R = (*mdt1)->globalPosition().perp();
+ float mdt3R = (*mdt3)->globalPosition().perp();
+ if(fabs(mdt1R-mdt3R) > 160) continue;
+ }
+ //store and fit the good combinations
+ std::vector<Muon::MdtPrepData*> mdts;
+ mdts.push_back((*mdt1));
+ mdts.push_back((*mdt2));
+ mdts.push_back((*mdt3));
+ std::vector<TrackletSegment> tmpSegs = TrackletSegmentFitter(mdts);
+ for(unsigned int i=0; i<tmpSegs.size(); ++i) mlsegments.push_back(tmpSegs.at(i));
+ }//end loop on mdt3
+ }//end loop on mdt2
+ }//end loop on mdt1
+
+ //store the reconstructed segments according to station, ML and sector
+ //station identifiers:
+ // 0 == BI 1 == BM 2 == BO
+ // 3 == EI 4 == EM 5 == EO
+ if(stName == 0 || stName == 1 || stName == 7 || stName == 52) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[0][ML-1][sector-1].push_back(mlsegments.at(k));
+ else if(stName == 2 || stName == 3 || stName == 8) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[1][ML-1][sector-1].push_back(mlsegments.at(k));
+ else if(stName == 4 || stName == 5 || stName == 9 || stName == 10) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[2][ML-1][sector-1].push_back(mlsegments.at(k));
+ else if(stName == 13 || stName == 49) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[3][ML-1][sector-1].push_back(mlsegments.at(k));
+ else if(stName == 17 || stName == 18) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[4][ML-1][sector-1].push_back(mlsegments.at(k));
+ else if(stName == 20 || stName == 21) for(unsigned int k=0; k<mlsegments.size(); ++k) segs[5][ML-1][sector-1].push_back(mlsegments.at(k));
+ else ATH_MSG_WARNING( "Found segments belonging to chamber " << stName << " that have not been stored" );
+ }//end loop on mdt chambers
+
+ //Combine/remove duplicate segments
+ std::vector<TrackletSegment> CleanSegs[6][2][16];
+ int nSegs(0);
+ for(int st=0; st<6; ++st) {
+ for(int ml=0; ml<2; ++ml) {
+ for(int sector=0; sector<16; ++sector) {
+ if(segs[st][ml][sector].size() > 0) {
+ CleanSegs[st][ml][sector] = CleanSegments(segs[st][ml][sector]);
+ nSegs += CleanSegs[st][ml][sector].size();
+ }
+ }
+ }
+ }
+
+ for(int st=0; st<6; ++st) {
+ m_DeltaAlphaCut = m_BarrelDeltaAlphaCut;
+ if(st >= 3) m_DeltaAlphaCut = m_EndcapDeltaAlphaCut;
+ for(int sector=0; sector<16; ++sector) {
+ for(unsigned int i1=0; i1<CleanSegs[st][0][sector].size(); ++i1) {
+ //Set the delta alpha cut depending on station type
+ int stName = CleanSegs[st][0][sector].at(i1).mdtChamber();
+ if( stName == 0) m_DeltaAlphaCut = k_BIL/750.0;
+ else if(stName == 2) m_DeltaAlphaCut = k_BML/750.0;
+ else if(stName == 3) m_DeltaAlphaCut = k_BMS/750.0;
+ else if(stName == 4) m_DeltaAlphaCut = k_BOL/750.0;
+ else if(stName == 7) m_DeltaAlphaCut = k_BIL/750.0;
+ else if(stName == 8) m_DeltaAlphaCut = k_BML/750.0;
+ else if(stName == 9) m_DeltaAlphaCut = k_BOL/750.0;
+ else if(stName == 10) m_DeltaAlphaCut = k_BOL/750.0;
+ else if(stName == 52) m_DeltaAlphaCut = k_BIL/750.0;
+ else if(stName <= 11) m_DeltaAlphaCut = 0.02;
+ else m_DeltaAlphaCut = m_EndcapDeltaAlphaCut;
+ //loop on ML2 segments from same sector
+ for(unsigned int i2=0; i2<CleanSegs[st][1][sector].size(); ++i2) {
+ if(CleanSegs[st][0][sector].at(i1).mdtChamber()!=CleanSegs[st][1][sector].at(i2).mdtChamber() ||
+ CleanSegs[st][0][sector].at(i1).mdtChEta()!=CleanSegs[st][1][sector].at(i2).mdtChEta()) continue;
+ float deltaAlpha = CleanSegs[st][0][sector].at(i1).alpha()-CleanSegs[st][1][sector].at(i2).alpha();
+ bool goodDeltab = DeltabCalc(CleanSegs[st][0][sector].at(i1),CleanSegs[st][1][sector].at(i2));
+ //select the good combinations
+ if(fabs(deltaAlpha) < m_DeltaAlphaCut && goodDeltab) {
+ if(st < 3) {//barrel chambers
+ float charge = 1;
+ if( deltaAlpha*CleanSegs[st][0][sector].at(i1).globalPosition().z()*tan(CleanSegs[st][0][sector].at(i1).alpha()) < 0 ) charge = -1;
+ float pTot = TrackMomentum(CleanSegs[st][0][sector].at(i1).mdtChamber(),deltaAlpha);
+ if(pTot < 800.) continue;
+ if(pTot >= 9999.) {
+ //if we find a straight track, try to do a global refit to minimize the number of duplicates
+ charge = 0;
+ std::vector<Muon::MdtPrepData*> mdts = CleanSegs[st][0][sector].at(i1).mdtHitsOnTrack();
+ std::vector<Muon::MdtPrepData*> mdts2 = CleanSegs[st][1][sector].at(i2).mdtHitsOnTrack();
+ for(unsigned int k=0; k<mdts2.size(); ++k) mdts.push_back(mdts2.at(k));
+ std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
+ if(CombinedSeg.size() > 0) {
+ //calculate momentum components & uncertainty
+ float Trk1overPErr = TrackMomentumError(CombinedSeg[0]);
+ float pT = pTot*sin(CombinedSeg[0].alpha());
+ float pz = pTot*cos(CombinedSeg[0].alpha());
+ Amg::Vector3D momentum(pT*cos(CombinedSeg[0].globalPosition().phi()),pT*sin(CombinedSeg[0].globalPosition().phi()),pz);
+ //create the error matrix
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(CombinedSeg[0].rError());//delta locR
+ matrix(1,1) = sq(CombinedSeg[0].zError());//delta locz
+ matrix(2,2) = sq(0.00000000001);//delta phi (~0 because we explicitly rotate all tracklets into the middle of the chamber)
+ matrix(3,3) = sq(CombinedSeg[0].alphaError());//delta theta
+ matrix(4,4) = sq(Trk1overPErr);//delta 1/p
+ Tracklet tmpTrk(CombinedSeg[0],momentum,matrix,charge);
+ ATH_MSG_DEBUG( "Track " << tracklets.size() <<" found with p = (" << momentum.x() << ", " << momentum.y() << ", " << momentum.z()
+ << ") and |p| = " << tmpTrk.momentum().mag() << " MeV" );
+ tracklets.push_back(tmpTrk);
+ }
+ }
+ else {
+ //tracklet has a measurable momentum
+ float Trk1overPErr = TrackMomentumError(CleanSegs[st][0][sector].at(i1),CleanSegs[st][1][sector].at(i2));
+ float pT = pTot*sin(CleanSegs[st][0][sector].at(i1).alpha());
+ float pz = pTot*cos(CleanSegs[st][0][sector].at(i1).alpha());
+ Amg::Vector3D momentum(pT*cos(CleanSegs[st][0][sector].at(i1).globalPosition().phi()),
+ pT*sin(CleanSegs[st][0][sector].at(i1).globalPosition().phi()),pz);
+ //create the error matrix
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(CleanSegs[st][0][sector].at(i1).rError());//delta locR
+ matrix(1,1) = sq(CleanSegs[st][0][sector].at(i1).zError());//delta locz
+ matrix(2,2) = sq(0.00000000001);//delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
+ matrix(3,3) = sq(CleanSegs[st][0][sector].at(i1).alphaError());//delta theta
+ matrix(4,4) = sq(Trk1overPErr);//delta 1/p
+ Tracklet tmpTrk(CleanSegs[st][0][sector].at(i1),CleanSegs[st][1][sector].at(i2),momentum,matrix,charge);
+ ATH_MSG_DEBUG( "Track " << tracklets.size() <<" found with p = (" << momentum.x() << ", " << momentum.y() << ", " << momentum.z()
+ << ") and |p| = " << tmpTrk.momentum().mag() << " MeV" );
+ tracklets.push_back(tmpTrk);
+ }
+ }//end barrel chamber selection
+ else if(st >= 3) {//endcap tracklets
+ //always straight tracklets (no momentum measurement possible)
+ std::vector<Muon::MdtPrepData*> mdts = CleanSegs[st][0][sector].at(i1).mdtHitsOnTrack();
+ std::vector<Muon::MdtPrepData*> mdts2 = CleanSegs[st][1][sector].at(i2).mdtHitsOnTrack();
+ for(unsigned int k=0; k<mdts2.size(); ++k) mdts.push_back(mdts2.at(k));
+ std::vector<TrackletSegment> CombinedSeg = TrackletSegmentFitter(mdts);
+ if(CombinedSeg.size() > 0) {
+ float charge = 0;
+ float pT = 100000.0*sin(CombinedSeg[0].alpha());
+ float pz = 100000.0*cos(CombinedSeg[0].alpha());
+ //create the error matrix
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(CombinedSeg[0].rError());//delta locR
+ matrix(1,1) = sq(CombinedSeg[0].zError());//delta locz
+ matrix(2,2) = sq(0.0000001);//delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
+ matrix(3,3) = sq(CombinedSeg[0].alphaError());//delta theta
+ matrix(4,4) = sq(0.00005);//delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
+ Amg::Vector3D momentum(pT*cos(CombinedSeg[0].globalPosition().phi()),
+ pT*sin(CombinedSeg[0].globalPosition().phi()),pz);
+ Tracklet tmpTrk(CombinedSeg[0],momentum,matrix,charge);
+ tracklets.push_back(tmpTrk);
+ }
+ }//end endcap tracklet selection
+
+ }//end tracklet selection (delta alpha & delta b)
+
+ }//end loop on ML2 segments
+ }//end loop on ML1 segments
+ }//end loop on sectors
+ }//end loop on stations
+
+ //Resolve any ambiguous tracklets
+ tracklets = ResolveAmbiguousTracklets(tracklets);
+
+ //convert from tracklets to Trk::Tracks
+ convertToTrackParticles(tracklets);
+
+ return StatusCode::SUCCESS;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ //convert tracklets to Trk::Track and store in a TrackCollection
+ void MSVertexTrackletTool::convertToTrackParticles(std::vector<Tracklet>& tracklets) {
+
+ xAOD::TrackParticleContainer* container = new xAOD::TrackParticleContainer();
+ if(evtStore()->record( container, m_TPContainer ).isFailure()) {
+ ATH_MSG_WARNING("Failed to record the xAOD::TrackParticleContainer with key " << m_TPContainer);
+ return;
+ }
+
+ xAOD::TrackParticleAuxContainer* aux = new xAOD::TrackParticleAuxContainer();
+ if(evtStore()->record( aux, m_TPContainer + "Aux" ).isFailure()) {
+ ATH_MSG_WARNING("Failed to record the xAOD::TrackParticleAuxContainer with key " << m_TPContainer << "Aux");
+ return;
+ }
+ container->setStore(aux);
+
+ for(std::vector<Tracklet>::iterator trkItr=tracklets.begin(); trkItr!=tracklets.end(); ++trkItr) {
+
+ //create the Trk::Perigee for the tracklet
+ AmgSymMatrix(5) * covariance = new AmgSymMatrix(5)(trkItr->errorMatrix());
+ Trk::Perigee * myPerigee = new Trk::Perigee(0.,0.,trkItr->momentum().phi(),trkItr->momentum().theta(),
+ trkItr->charge()/trkItr->momentum().mag(),
+ Trk::PerigeeSurface(trkItr->globalPosition()),covariance);
+
+ //create, store & define the xAOD::TrackParticle
+ xAOD::TrackParticle* trackparticle = new xAOD::TrackParticle();
+
+ container->push_back( trackparticle );
+
+ trackparticle->setFitQuality(1.,(float)trkItr->mdtHitsOnTrack().size());
+ trackparticle->setTrackProperties(xAOD::TrackProperties::LowPtTrack);
+ trackparticle->setDefiningParameters(myPerigee->parameters()[Trk::d0],
+ myPerigee->parameters()[Trk::z0],
+ myPerigee->parameters()[Trk::phi0],
+ myPerigee->parameters()[Trk::theta],
+ myPerigee->parameters()[Trk::qOverP]);
+
+ std::vector<float> covMatrixVec;
+ Amg::compress(*covariance,covMatrixVec);
+ trackparticle->setDefiningParametersCovMatrixVec(covMatrixVec);
+
+ //cleanup memory
+ delete myPerigee;
+ }
+ return;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ int MSVertexTrackletTool::SortMDThits(std::vector<std::vector<Muon::MdtPrepData*> >& SortedMdt) {
+
+ SortedMdt.clear();
+ int nMDT(0);
+ const Muon::MdtPrepDataContainer* mdtTES = 0;
+ if(evtStore()->retrieve(mdtTES,"MDT_DriftCircles").isFailure()) {
+ ATH_MSG_FATAL( "No MdtPrepDataContainer found!!" );
+ return 0;
+ }
+ Muon::MdtPrepDataContainer::const_iterator MDTItr = mdtTES->begin();
+ Muon::MdtPrepDataContainer::const_iterator MDTItrE = mdtTES->end();
+
+ for(; MDTItr != MDTItrE; ++MDTItr) {
+
+ Muon::MdtPrepDataCollection::const_iterator mpdc = (*MDTItr)->begin();
+ Muon::MdtPrepDataCollection::const_iterator mpdcE = (*MDTItr)->end();
+
+ int stName = m_mdtIdHelper->stationName((*mpdc)->identify());
+
+ // Doesn't consider hits belonging to chambers BEE, EEL and EES
+ if(stName == 6 || stName == 14 || stName == 15) continue;
+
+ // Doesn't consider hits belonging to chambers BIS7 and BIS8
+ if(stName == 1 && fabs(m_mdtIdHelper->stationEta((*mpdc)->identify())) >= 7) continue;
+
+ for(; mpdc != mpdcE; ++mpdc) {
+
+ // Removes noisy hits
+ if((*mpdc)->adc() < 50) continue;
+
+ // Removes dead modules or out of time hits
+ if((*mpdc)->status() != 1) continue;
+
+ // Removes tubes out of readout during drift time
+ if((*mpdc)->localPosition()[Trk::locR] == 0.) continue;
+
+ // Should not happen but is seen in the barrel
+ if((*mpdc)->tdc() < 300 || (*mpdc)->tdc() > 2100) continue;
+
+ // Doesn't include hits on tubes that have multiple, good hits
+ Identifier mpdc_id = (*mpdc)->identify();
+ int hitIsUnique = 1;
+ Muon::MdtPrepDataContainer::const_iterator MDTItr2= mdtTES->begin();
+ for(;MDTItr2 != mdtTES->end(); MDTItr2++) {
+ Muon::MdtPrepDataCollection::const_iterator mpdc2 = (*MDTItr)->begin();
+ Muon::MdtPrepDataCollection::const_iterator mpdcE2 = (*MDTItr)->end();
+ for(;mpdc2 != mpdcE2; mpdc2++) {
+ if(mpdc == mpdc2) continue;
+ if((*mpdc2)->adc() < 50) continue;
+ if((*mpdc2)->status() != 1) continue;
+ if((*mpdc2)->localPosition()[Trk::locR] == 0.) continue;
+
+ if((*mpdc2)->identify() == mpdc_id) {
+ hitIsUnique = 0;
+ }
+ }
+ }
+ if(!hitIsUnique) continue;
+
+ nMDT++;
+
+ int foundChamberML = 0;
+ for(unsigned int i = 0; i < SortedMdt.size(); i++) {
+ if(SortedMdt.size() == 0) continue;
+
+ Identifier id1 = SortedMdt.at(i).at(0)->identify();
+ Identifier id2 = (*mpdc)->identify();
+ if(SortMDT(id1, id2)) {
+ foundChamberML = 1;
+ SortedMdt.at(i).push_back(*mpdc);
+ }
+ }
+ if(!foundChamberML) {
+ std::vector<Muon::MdtPrepData*> tempMdt;
+ tempMdt.push_back((*mpdc));
+ SortedMdt.push_back(tempMdt);
+ }
+
+ }//end MdtPrepDataCollection
+ }//end MdtPrepDataContaier
+ //loop over ML, remove any with occupancy > 75%
+ for(std::vector<std::vector<Muon::MdtPrepData*> >::iterator mlIt=SortedMdt.begin(); mlIt!=SortedMdt.end(); ++mlIt) {
+ Identifier id = (*mlIt)[0]->identify();
+ float nTubes = (m_mdtIdHelper->tubeMax(id) - m_mdtIdHelper->tubeMin(id))*(m_mdtIdHelper->tubeLayerMax(id))*(m_mdtIdHelper->multilayerMax(id));
+ float frac = (*mlIt).size()/nTubes;
+ if(frac > 0.75) (*mlIt).clear();
+ else std::sort(mlIt->begin(),mlIt->end(),MSVertexTrackletTool::mdtComp);//sort the MDTs by layer and tube number
+ }
+
+ return nMDT;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ bool MSVertexTrackletTool::SortMDT(Identifier& i1, Identifier& i2) {
+ if(m_mdtIdHelper->stationName(i1) != m_mdtIdHelper->stationName(i2)) return false;
+ if(m_mdtIdHelper->stationEta(i1) != m_mdtIdHelper->stationEta(i2)) return false;
+ if(m_mdtIdHelper->stationPhi(i1) != m_mdtIdHelper->stationPhi(i2)) return false;
+ if(m_mdtIdHelper->multilayer(i1) != m_mdtIdHelper->multilayer(i2)) return false;
+ return true;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ bool MSVertexTrackletTool::mdtComp(const Muon::MdtPrepData* mprd1, const Muon::MdtPrepData* mprd2) {
+ //mdts sorted by layer and tube number
+ if(mdtCompareIdHelper->tubeLayer(mprd1->identify()) > mdtCompareIdHelper->tubeLayer(mprd2->identify())) return false;
+ if(mdtCompareIdHelper->tubeLayer(mprd1->identify()) < mdtCompareIdHelper->tubeLayer(mprd2->identify())) return true;
+ if(mdtCompareIdHelper->tube(mprd1->identify()) < mdtCompareIdHelper->tube(mprd2->identify())) return true;
+ return false;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<TrackletSegment> MSVertexTrackletTool::TrackletSegmentFitter(std::vector<Muon::MdtPrepData*>& mdts) {
+
+ //create the segment seeds
+ std::vector<std::pair<float,float> > SeedParams = SegSeeds(mdts);
+ //fit the segments
+ std::vector<TrackletSegment> segs = TrackletSegmentFitterCore(mdts,SeedParams);
+
+ return segs;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<std::pair<float,float> > MSVertexTrackletTool::SegSeeds(std::vector<Muon::MdtPrepData*>& mdts) {
+
+ std::vector<std::pair<float,float> > SeedParams;
+ //create seeds by drawing the 4 possible lines tangent to the two outermost drift circles
+ //see http://cds.cern.ch/record/620198 (section 4.3) for description of the algorithm
+ //keep all seeds which satisfy the criterion: residual(mdt 2) < m_SeedResidual
+ //NOTE: here there is an assumption that each MDT has a radius of 30mm
+ // -- needs to be revisited when the small tubes in sectors 12 & 14 are installed
+ float x1 = mdts.front()->globalPosition().z();
+ float y1 = mdts.front()->globalPosition().perp();
+ float r1 = fabs(mdts.front()->localPosition()[Trk::locR]);
+
+ float x2 = mdts.back()->globalPosition().z();
+ float y2 = mdts.back()->globalPosition().perp();
+ float r2 = fabs(mdts.back()->localPosition()[Trk::locR]);
+
+ float DeltaX = x2 - x1;
+ float DeltaY = y2 - y1;
+ float DistanceOfCenters = sqrt(DeltaX*DeltaX + DeltaY*DeltaY);
+ if(DistanceOfCenters < 30) return SeedParams;
+ float Alpha0 = acos(DeltaX/DistanceOfCenters);
+
+ //First seed
+ float phi = mdts.front()->globalPosition().phi();
+ float RSum = r1 + r2;
+ if(RSum > DistanceOfCenters) return SeedParams;
+ float Alpha1 = asin(RSum/DistanceOfCenters);
+ float line_theta = Alpha0 + Alpha1;
+ float z_line = x1 + r1*sin(line_theta);
+ float rho_line = y1 - r1*cos(line_theta);
+
+ Amg::Vector3D gPos1(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D gDir(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ Amg::Vector3D globalDir1(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope1 = (globalDir1.perp()/globalDir1.z());
+ float gInter1 = gPos1.perp() - gSlope1*gPos1.z();
+ float resid = SeedResiduals(mdts,gSlope1,gInter1);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope1,gInter1) );
+ //Second seed
+ line_theta = Alpha0 - Alpha1;
+ z_line = x1 - r1*sin(line_theta);
+ rho_line = y1 + r1*cos(line_theta);
+ Amg::Vector3D gPos2(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D globalDir2(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope2= (globalDir2.perp()/globalDir2.z());
+ float gInter2 = gPos2.perp() - gSlope2*gPos2.z();
+ resid = SeedResiduals(mdts,gSlope2,gInter2);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope2,gInter2) );
+
+ float Alpha2 = asin(fabs(r2-r1)/DistanceOfCenters);
+ if(r1 < r2) {
+ //Third seed
+ line_theta = Alpha0 +Alpha2;
+ z_line = x1 - r1*sin(line_theta);
+ rho_line = y1 + r1*cos(line_theta);
+
+ Amg::Vector3D gPos3(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D globalDir3(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope3= (globalDir3.perp()/globalDir3.z());
+ float gInter3 = gPos3.perp() - gSlope3*gPos3.z();
+ resid = SeedResiduals(mdts,gSlope3,gInter3);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope3,gInter3) );
+
+ //Fourth seed
+ line_theta = Alpha0 - Alpha2;
+ z_line = x1 + r1*sin(line_theta);
+ rho_line = y1 - r1*cos(line_theta);
+
+ Amg::Vector3D gPos4(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D globalDir4(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope4= (globalDir4.perp()/globalDir4.z());
+ float gInter4 = gPos4.perp() - gSlope4*gPos4.z();
+ resid = SeedResiduals(mdts,gSlope4,gInter4);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope4,gInter4) );
+ }
+ else {
+ //Third seed
+ line_theta = Alpha0 +Alpha2;
+ z_line = x1 + r1*sin(line_theta);
+ rho_line = y1 - r1*cos(line_theta);
+
+ Amg::Vector3D gPos3(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D globalDir3(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope3= (globalDir3.perp()/globalDir3.z());
+ float gInter3 = gPos3.perp() - gSlope3*gPos3.z();
+ resid = SeedResiduals(mdts,gSlope3,gInter3);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope3,gInter3) );
+
+ //Fourth seed
+ line_theta = Alpha0 - Alpha2;
+ z_line = x1 - r1*sin(line_theta);
+ rho_line = y1 + r1*cos(line_theta);
+
+ Amg::Vector3D gPos4(rho_line*cos(phi), rho_line*sin(phi), z_line);
+ Amg::Vector3D globalDir4(cos(phi)*sin(line_theta), sin(phi)*sin(line_theta), cos(line_theta));
+ float gSlope4= (globalDir4.perp()/globalDir4.z());
+ float gInter4 = gPos4.perp() - gSlope4*gPos4.z();
+ resid = SeedResiduals(mdts,gSlope4,gInter4);
+ if(resid < m_SeedResidual) SeedParams.push_back( std::pair<float,float>(gSlope4,gInter4) );
+ }
+
+ return SeedParams;
+ }
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ float MSVertexTrackletTool::SeedResiduals(std::vector<Muon::MdtPrepData*>& mdts, float slope, float inter) {
+ //calculate the residual of the MDTs not used to create the seed
+ float resid = 0;
+ for(unsigned int i=1; i<(mdts.size()-1); ++i) {
+ float mdtR = mdts.at(i)->globalPosition().perp();
+ float mdtZ = mdts.at(i)->globalPosition().z();
+ float res = fabs( (mdts.at(i)->localPosition()[Trk::locR] - fabs((mdtR-inter-slope*mdtZ)/sqrt(sq(slope)+1)))/( Amg::error(mdts.at(i)->localCovariance(),Trk::locR) ) );
+ if(res > resid) resid = res;
+ }
+ return resid;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<TrackletSegment> MSVertexTrackletTool::TrackletSegmentFitterCore(std::vector<Muon::MdtPrepData*>& mdts,std::vector<std::pair<float,float> >& SeedParams){
+ std::vector<TrackletSegment> segs;
+ int stName = m_mdtIdHelper->stationName(mdts.at(0)->identify());
+ float mlmidpt = 0;
+ if(stName <= 11 || stName == 52) mlmidpt = sqrt(sq(mdts.at(0)->detectorElement()->center().x())+sq(mdts.at(0)->detectorElement()->center().y()));
+ else if(stName <= 21 || stName == 49) mlmidpt = mdts.at(0)->detectorElement()->center().z();
+ else return segs;
+ for(unsigned int i_p=0; i_p<SeedParams.size(); ++i_p) {
+ //Min chi^2 fit from "Precision of the ATLAS Muon Spectrometer" -- M. Woudstra
+ //http://cds.cern.ch/record/620198?ln=en (section 4.3)
+ float chi2(0);
+ float s(0),sz(0),sy(0);
+ //loop on the mdt hits, find the weighted center
+ for(unsigned int i=0; i<mdts.size(); ++i) {
+ float mdt_y = sqrt(sq(mdts.at(i)->globalPosition().x())+sq(mdts.at(i)->globalPosition().y()));
+ float mdt_z = mdts.at(i)->globalPosition().z();
+ float sigma2 = sq(Amg::error(mdts.at(i)->localCovariance(),Trk::locR));
+ s += 1/sigma2;
+ sz += mdt_z/sigma2;
+ sy += mdt_y/sigma2;
+ }
+ float yc = sy/s;
+ float zc = sz/s;
+
+ //Find the initial parameters of the fit
+ float alpha = atan2(SeedParams.at(i_p).first,1.0);
+ if(alpha < 0) alpha += PI;
+ float dalpha = 0;
+ float d = (SeedParams.at(i_p).second - yc + zc*SeedParams.at(i_p).first)*cos(alpha);
+ float dd = 0;
+
+ //require segments to point to the second ML
+ if( fabs(cos(alpha)) > 0.97 && (stName <= 11 || stName == 52) ) continue;
+ if( fabs(cos(alpha)) < 0.03 && ((stName > 11 && stName < 22) || stName == 49) )continue;
+
+ //calculate constants used in the fit
+ float sPz(0),sPy(0),sPyy(0),sPzz(0),sPyz(0),sPyyzz(0);
+ for(unsigned int i=0; i<mdts.size(); ++i) {
+ float mdt_y = sqrt(sq(mdts.at(i)->globalPosition().x())+sq(mdts.at(i)->globalPosition().y()));
+ float mdt_z = mdts.at(i)->globalPosition().z();
+ float sigma2 = sq(Amg::error(mdts.at(i)->localCovariance(),Trk::locR));
+ sPz += (mdt_z-zc)/sigma2;
+ sPy += (mdt_y-yc)/sigma2;
+ sPyy += sq(mdt_y-yc)/sigma2;
+ sPzz += sq(mdt_z-zc)/sigma2;
+ sPyz += (mdt_y-yc)*(mdt_z-zc)/sigma2;
+ sPyyzz += ((mdt_y-yc)-(mdt_z-zc))*((mdt_y-yc)+(mdt_z-zc))/sigma2;
+ }
+
+ //iterative fit
+ int Nitr = 0;
+ float deltaAlpha = 0;
+ float deltad = 0;
+ while(true) {
+ float sumRyi(0),sumRzi(0),sumRi(0);
+ chi2 = 0;
+ Nitr++;
+ for(unsigned int i=0; i<mdts.size(); ++i) {
+ float mdt_y = sqrt(sq(mdts.at(i)->globalPosition().x())+sq(mdts.at(i)->globalPosition().y()));
+ float mdt_z = mdts.at(i)->globalPosition().z();
+ float yPi = -(mdt_z-zc)*sin(alpha) + (mdt_y-yc)*cos(alpha) - d;
+ float signR = -1*yPi/fabs(yPi);
+ float sigma2 = sq(Amg::error(mdts.at(i)->localCovariance(),Trk::locR));
+ float ri = signR*mdts.at(i)->localPosition()[Trk::locR];
+ ////
+ sumRyi += ri*(mdt_y-yc)/sigma2;
+ sumRzi += ri*(mdt_z-zc)/sigma2;
+ sumRi += ri/sigma2;
+ //
+ chi2 += sq(yPi+ri)/sigma2;
+ }
+ float bAlpha = -1*sPyz + cos(alpha)*(sin(alpha)*sPyyzz +2*cos(alpha)*sPyz + sumRzi) + sin(alpha)*sumRyi;
+ float AThTh = sPyy + cos(alpha)*(2*sin(alpha)*sPyz - cos(alpha)*sPyyzz);
+ //the new alpha & d parameters
+ float alphaNew = alpha + bAlpha/AThTh;
+
+ float dNew = sumRi/s;
+ //the errors
+ dalpha = sqrt(1/AThTh);
+ dd = sqrt(1/s);
+ deltaAlpha = fabs(alphaNew-alpha);
+ deltad = fabs(d-dNew);
+ //test if the new segment is different than the previous
+ if(deltaAlpha < 0.0000005 && deltad < 0.000005) break;
+ alpha = alphaNew;
+ d = dNew;
+ //Guard against infinite loops
+ if(Nitr > 10) break;
+ }//end while loop
+
+ //find the chi^2 probability of the segment
+ float chi2Prob = TMath::Prob(chi2,mdts.size()-2);
+ //keep only "good" segments
+ if(chi2Prob > m_minSegFinderChi2) {
+ float z0 = zc - d*sin(alpha);
+ float dz0 = sqrt(sq(dd*sin(alpha))+sq(d*dalpha*cos(alpha)));
+ float y0 = yc + d*cos(alpha);
+ float dy0 = sqrt(sq(dd*cos(alpha))+sq(d*dalpha*sin(alpha)));
+ //find the hit pattern, which side of the wire did the particle pass? (1==Left, 2==Right)
+ /*
+ ( )(/O)( )
+ (./)( )( ) == RRL == 221
+ (O/)( )( )
+ */
+ int pattern(0);
+ if(mdts.size() > 8) pattern = -1;//with more then 8 MDTs the pattern is unique
+ else {
+ for(unsigned int k=0; k<mdts.size(); ++k) {
+ int base = pow(10,k);
+ float mdtR = sqrt(sq(mdts.at(k)->globalPosition().x())+sq(mdts.at(k)->globalPosition().y()));
+ float mdtZ = mdts.at(k)->globalPosition().z();
+ float zTest = (mdtR-y0)/tan(alpha) + z0 - mdtZ;
+ if(zTest > 0) pattern += 2*base;
+ else pattern += base;
+ }
+ }
+
+ //information about the MDT chamber containing the tracklet
+ int chEta = m_mdtIdHelper->stationEta(mdts.at(0)->identify());
+ int chPhi = m_mdtIdHelper->stationPhi(mdts.at(0)->identify());
+
+ //find the position of the tracklet in the global frame
+ float mdtPhi = mdts.at(0)->globalPosition().phi();
+ Amg::Vector3D segpos(y0*cos(mdtPhi),y0*sin(mdtPhi),z0);
+ //create the tracklet
+ TrackletSegment MyTrackletSegment(stName,chEta,chPhi,mlmidpt,alpha,dalpha,segpos,dy0,dz0,mdts,pattern);
+ segs.push_back(MyTrackletSegment);
+ if(pattern == -1) break;//stop if we find a segment with more than 8 hits (guaranteed to be unique!)
+ }
+ }//end loop on segment seeds
+
+ //in case more than 1 segment is reconstructed, check if there are duplicates using the hit patterns
+ if(segs.size() > 1) {
+ std::vector<TrackletSegment> tmpSegs;
+ for(unsigned int i1=0; i1<segs.size(); ++i1) {
+ bool isUnique = true;
+ int pattern1 = segs.at(i1).getHitPattern();
+ for(unsigned int i2=(i1+1); i2<segs.size(); ++i2) {
+ if(pattern1 == -1) break;
+ int pattern2 = segs.at(i2).getHitPattern();
+ if(pattern1 == pattern2) isUnique = false;
+ }
+ if(isUnique) tmpSegs.push_back(segs.at(i1));
+ }
+ segs = tmpSegs;
+ }
+
+ //return the unique segments
+ return segs;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<TrackletSegment> MSVertexTrackletTool::CleanSegments(std::vector<TrackletSegment>& Segs) {
+ std::vector<TrackletSegment> CleanSegs;
+ std::vector<TrackletSegment> segs = Segs;
+ bool keepCleaning(true);
+ int nItr(0);
+ while(keepCleaning) {
+ nItr++;
+ keepCleaning = false;
+ for(std::vector<TrackletSegment>::iterator it=segs.begin(); it!=segs.end(); ++it) {
+ if(it->isCombined()) continue;
+ std::vector<TrackletSegment> segsToCombine;
+ float tanTh1 = tan(it->alpha());
+ float r1 = it->globalPosition().perp();
+ float zi1 = it->globalPosition().z() - r1/tanTh1;
+ //find all segments with similar parameters & attempt to combine
+ for(std::vector<TrackletSegment>::iterator sit=(it+1); sit!=segs.end(); ++sit) {
+ if(sit->isCombined()) continue;
+ if(it->mdtChamber() != sit->mdtChamber()) continue;//require the segments are in the same chamber
+ if( (it->mdtChEta())*(sit->mdtChEta()) < 0 ) continue;//check both segments are on the same side of the detector
+ if( it->mdtChPhi() != sit->mdtChPhi() ) continue;//in the same sector
+ if(fabs(it->alpha() - sit->alpha()) > 0.005) continue;//same trajectory
+ float tanTh2 = tan(sit->alpha());
+ float r2 = sit->globalPosition().perp();
+ float zi2 = sit->globalPosition().z() - r2/tanTh2;
+ //find the distance at the midpoint between the two segments
+ float rmid = (r1+r2)/2.;
+ float z1 = rmid/tanTh1 + zi1;
+ float z2 = rmid/tanTh2 + zi2;
+ float zdist = fabs(z1 - z2);
+ if(zdist < 0.5) {
+ segsToCombine.push_back( *sit );
+ sit->isCombined(true);
+ }
+ }//end sit loop
+
+ //if the segment is unique, keep it
+ if(segsToCombine.size() == 0) {
+ CleanSegs.push_back( *it );
+ }
+ //else, combine all like segments & refit
+ else if(segsToCombine.size() >= 1 ) {
+ //create a vector of all unique MDT hits in the segments
+ std::vector<Muon::MdtPrepData*> mdts = it->mdtHitsOnTrack();
+ for(unsigned int i=0; i<segsToCombine.size(); ++i) {
+ std::vector<Muon::MdtPrepData*> tmpmdts = segsToCombine[i].mdtHitsOnTrack();
+ for(std::vector<Muon::MdtPrepData*>::iterator mit=tmpmdts.begin(); mit!=tmpmdts.end(); ++mit) {
+ bool isNewHit(true);
+ for(std::vector<Muon::MdtPrepData*>::iterator msit=mdts.begin(); msit!=mdts.end(); ++msit) {
+ if((*mit)->identify() == (*msit)->identify()) {
+ isNewHit = false;
+ break;
+ }
+ }
+ if(isNewHit) mdts.push_back( *mit );
+ }
+ }//end segsToCombine loop
+
+ //only need to combine if there are extra hits added to the first segment
+ if(mdts.size() > it->mdtHitsOnTrack().size()) {
+ std::vector<TrackletSegment> refitsegs = TrackletSegmentFitter(mdts);
+ //if the refit fails, what to do?
+ if(refitsegs.size() == 0) {
+ if(segsToCombine.size() == 1) {
+ segsToCombine[0].isCombined(false);
+ CleanSegs.push_back( *it );
+ CleanSegs.push_back(segsToCombine[0]);
+ }
+ else {
+ //loop on the mdts and count the number of segments that share that hit
+ std::vector<int> nSeg;
+ for(unsigned int i=0; i<mdts.size(); ++i) {
+ nSeg.push_back(0);
+ //hit belongs to the first segment
+ for(unsigned int k=0; k<it->mdtHitsOnTrack().size(); ++k) {
+ if(it->mdtHitsOnTrack()[k]->identify() == mdts[i]->identify()) {
+ nSeg[i]++;
+ break;
+ }
+ }
+ //hit belongs to one of the duplicate segments
+ for(unsigned int k=0; k<segsToCombine.size(); ++k) {
+ for(unsigned int m=0; m<segsToCombine[k].mdtHitsOnTrack().size(); ++m) {
+ if(segsToCombine[k].mdtHitsOnTrack()[m]->identify() == mdts[i]->identify()) {
+ nSeg[i]++;
+ break;
+ }
+ }//end loop on mdtHitsOnTrack
+ }//end loop on segsToCombine
+ }//end loop on mdts
+
+ //loop over the duplicates and remove the MDT used by the fewest segments until the fit converges
+ bool keeprefitting(true);
+ int nItr2(0);
+ while(keeprefitting) {
+ nItr2++;
+ int nMinSeg(nSeg[0]);
+ Muon::MdtPrepData* minmdt = mdts[0];
+ std::vector<int> nrfsegs;
+ std::vector<Muon::MdtPrepData*> refitmdts;
+ //loop on MDTs, identify the overlapping set of hits
+ for(unsigned int i=1; i<mdts.size(); ++i) {
+ if(nSeg[i] < nMinSeg) {
+ refitmdts.push_back(minmdt);
+ nrfsegs.push_back(nMinSeg);
+ minmdt = mdts[i];
+ nMinSeg = nSeg[i];
+ }
+ else {
+ refitmdts.push_back(mdts[i]);
+ nrfsegs.push_back(nSeg[i]);
+ }
+ }
+ //reset the list of MDTs & the minimum number of segments an MDT must belong to
+ mdts = refitmdts;
+ nSeg = nrfsegs;
+ //try to fit the new set of MDTs
+ refitsegs = TrackletSegmentFitter(mdts);
+ if(refitsegs.size() > 0) {
+ for(std::vector<TrackletSegment>::iterator rfit=refitsegs.begin(); rfit!=refitsegs.end(); ++rfit) {
+ CleanSegs.push_back( *rfit );
+ }
+ keeprefitting = false;//stop refitting if segments are found
+ }
+ else if(mdts.size() <= 3) {
+ CleanSegs.push_back( *it );
+ keeprefitting = false;
+ }
+ if(nItr2 > 10) break;
+ }//end while
+ }
+ }
+ else {
+ keepCleaning = true;
+ for(std::vector<TrackletSegment>::iterator rfit=refitsegs.begin(); rfit!=refitsegs.end(); ++rfit) {
+ CleanSegs.push_back( *rfit );
+ }
+ }
+ }
+ //if there are no extra MDT hits, keep only the first segment as unique
+ else CleanSegs.push_back( *it );
+ }
+ }//end it loop
+ if(keepCleaning) {
+ segs = CleanSegs;
+ CleanSegs.clear();
+ }
+ if(nItr > 10) break;
+ }//end while
+
+ return CleanSegs;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ bool MSVertexTrackletTool::DeltabCalc(TrackletSegment& ML1seg, TrackletSegment& ML2seg) {
+ float ChMid = (ML1seg.getChMidPoint() + ML2seg.getChMidPoint())/2.0;
+ //Calculate the Delta b (see http://inspirehep.net/record/1266438)
+ float mid1(100),mid2(1000);
+ float deltab(100);
+ if(ML1seg.mdtChamber() <= 11 || ML1seg.mdtChamber() == 52) {
+ //delta b in the barrel
+ mid1 = (ChMid - ML1seg.globalPosition().perp())/tan(ML1seg.alpha()) + ML1seg.globalPosition().z();
+ mid2 = (ChMid - ML2seg.globalPosition().perp())/tan(ML2seg.alpha()) + ML2seg.globalPosition().z();
+ float r01 = ML1seg.globalPosition().perp() - ML1seg.globalPosition().z()*tan(ML1seg.alpha());
+ float r02 = ML2seg.globalPosition().perp() - ML2seg.globalPosition().z()*tan(ML2seg.alpha());
+ deltab = (mid2*tan(ML1seg.alpha()) - ChMid + r01)/(sqrt(1+sq(tan(ML1seg.alpha()))));
+ float deltab2 = (mid1*tan(ML2seg.alpha()) - ChMid + r02)/(sqrt(1+sq(tan(ML2seg.alpha()))));
+ if(fabs(deltab2) < fabs(deltab)) deltab = deltab2;
+ }
+ else {
+ //delta b in the endcap
+ mid1 = ML1seg.globalPosition().perp() + tan(ML1seg.alpha())*(ChMid - ML1seg.globalPosition().z());
+ mid2 = ML2seg.globalPosition().perp() + tan(ML2seg.alpha())*(ChMid - ML2seg.globalPosition().z());
+ float z01 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp()/tan(ML1seg.alpha());
+ float z02 = ML1seg.globalPosition().z() - ML1seg.globalPosition().perp()/tan(ML1seg.alpha());
+ deltab = (mid2/tan(ML1seg.alpha()) - ChMid +z01)/(sqrt(1+sq(1/tan(ML1seg.alpha()))));
+ float deltab2 = (mid1/tan(ML2seg.alpha()) - ChMid +z02)/(sqrt(1+sq(1/tan(ML2seg.alpha()))));
+ if(fabs(deltab2) < fabs(deltab)) deltab = deltab2;
+ }
+
+ //calculate the maximum allowed Delta b based on delta alpha uncertainties and ML spacing
+ double dbmax = 5*fabs(ChMid-ML1seg.getChMidPoint())*sqrt(sq(ML1seg.alphaError()) + sq(ML2seg.alphaError()));
+ if(dbmax > m_maxDeltabCut) dbmax = m_maxDeltabCut;
+ if(fabs(deltab) < dbmax) return true;
+ else return false;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ float MSVertexTrackletTool::TrackMomentum(int chamber,float deltaAlpha) {
+ float pTot = 100000.;
+ //p = k/delta_alpha
+ if(chamber == 0) pTot = k_BIL/fabs(deltaAlpha);
+ else if(chamber == 2) pTot = k_BML/fabs(deltaAlpha);
+ else if(chamber == 3) pTot = k_BMS/fabs(deltaAlpha);
+ else if(chamber == 4) pTot = k_BOL/fabs(deltaAlpha);
+ else if(chamber == 7) pTot = k_BIL/fabs(deltaAlpha);
+ else if(chamber == 8) pTot = k_BML/fabs(deltaAlpha);
+ else if(chamber == 9) pTot = k_BOL/fabs(deltaAlpha);
+ else if(chamber == 10) pTot = k_BOL/fabs(deltaAlpha);
+ else if(chamber == 52) pTot = k_BIL/fabs(deltaAlpha);
+ if(pTot > 10000.) pTot = 100000.;
+
+ return pTot;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ float MSVertexTrackletTool::TrackMomentumError(TrackletSegment& ml1, TrackletSegment& ml2) {
+ //uncertainty on 1/p
+ int ChType = ml1.mdtChamber();
+ float dalpha = sqrt(sq(ml1.alphaError())+sq(ml2.alphaError()));
+ float pErr = dalpha/k_BML;
+ if(ChType == 0) pErr = dalpha/k_BIL;
+ else if(ChType == 2) pErr = dalpha/k_BML;
+ else if(ChType == 3) pErr = dalpha/k_BMS;
+ else if(ChType == 4) pErr = dalpha/k_BOL;
+ else if(ChType == 7) pErr = dalpha/k_BIL;
+ else if(ChType == 8) pErr = dalpha/k_BML;
+ else if(ChType == 9) pErr = dalpha/k_BOL;
+ else if(ChType == 10) pErr = dalpha/k_BOL;
+ else if(ChType == 52) pErr = dalpha/k_BIL;
+
+ return pErr;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ float MSVertexTrackletTool::TrackMomentumError(TrackletSegment& ml1) {
+ //uncertainty in 1/p
+ int ChType = ml1.mdtChamber();
+ float dalpha = fabs(ml1.alphaError());
+ float pErr = dalpha/k_BML;
+ if(ChType == 0) pErr = dalpha/k_BIL;
+ else if(ChType == 2) pErr = dalpha/k_BML;
+ else if(ChType == 3) pErr = dalpha/k_BMS;
+ else if(ChType == 4) pErr = dalpha/k_BOL;
+ else if(ChType == 7) pErr = dalpha/k_BIL;
+ else if(ChType == 8) pErr = dalpha/k_BML;
+ else if(ChType == 9) pErr = dalpha/k_BOL;
+ else if(ChType == 10) pErr = dalpha/k_BOL;
+ else if(ChType == 52) pErr = dalpha/k_BIL;
+
+ return pErr;
+ }
+
+
+//** ----------------------------------------------------------------------------------------------------------------- **//
+
+
+ std::vector<Tracklet> MSVertexTrackletTool::ResolveAmbiguousTracklets(std::vector<Tracklet>& tracks) {
+
+ ATH_MSG_DEBUG( "In ResolveAmbiguousTracks" );
+
+ //////////////////////////////
+ // considering only tracklets with the number of associated hits
+ // being more than 3/4 the number of layers in the MS chamber
+
+ if (m_tightTrackletRequirement) {
+
+ std::vector<Tracklet> myTracks = tracks;
+ tracks.clear();
+
+ for(unsigned int tk1=0; tk1<myTracks.size(); ++tk1) {
+ Identifier id1 = myTracks.at(tk1).getML1seg().mdtHitsOnTrack().at(0)->identify();
+ Identifier id2 = myTracks.at(tk1).getML2seg().mdtHitsOnTrack().at(0)->identify();
+ int nLayerML1 = m_mdtIdHelper->tubeLayerMax(id1);
+ int nLayerML2 = m_mdtIdHelper->tubeLayerMax(id2);
+ float ratio = (float)(myTracks.at(tk1).mdtHitsOnTrack().size())/(nLayerML1+nLayerML2);
+ if (ratio>0.75) tracks.push_back(myTracks.at(tk1));
+ }
+ }
+
+ std::vector<Tracklet> UniqueTracks;
+ std::vector<unsigned int> AmbigTrks;
+ int nBarrelAmbiguousTrks = 0;
+ int nEndcapAmbiguousTrks = 0;
+ for(unsigned int tk1=0; tk1<tracks.size(); ++tk1) {
+ int nShared = 0;
+ //check if any Ambiguity has been broken
+ bool isResolved = false;
+ for(unsigned int ck=0; ck<AmbigTrks.size(); ++ck) {
+ if(tk1 == AmbigTrks.at(ck)) {
+ isResolved = true;
+ break;
+ }
+ }
+ if(isResolved) continue;
+ std::vector<Tracklet> AmbigTracks;
+ AmbigTracks.push_back(tracks.at(tk1));
+ //get a point on the track
+ float Trk1ML1R = tracks.at(tk1).getML1seg().globalPosition().perp();
+ float Trk1ML1Z = tracks.at(tk1).getML1seg().globalPosition().z();
+ float Trk1ML2R = tracks.at(tk1).getML2seg().globalPosition().perp();
+ float Trk1ML2Z = tracks.at(tk1).getML2seg().globalPosition().z();
+
+ //loop over the rest of the tracks and find any amibuities
+ for(unsigned int tk2=(tk1+1); tk2<tracks.size(); ++tk2) {
+ if(tracks.at(tk1).mdtChamber() == tracks.at(tk2).mdtChamber() &&
+ tracks.at(tk1).mdtChPhi() == tracks.at(tk2).mdtChPhi() &&
+ (tracks.at(tk1).mdtChEta())*(tracks.at(tk2).mdtChEta()) > 0 ) {
+ //check if any Ambiguity has been broken
+ for(unsigned int ck=0; ck<AmbigTrks.size(); ++ck) {
+ if(tk2 == AmbigTrks.at(ck)) {
+ isResolved = true;
+ break;
+ }
+ }
+ if(isResolved) continue;
+ //get a point on the track
+ float Trk2ML1R = tracks.at(tk2).getML1seg().globalPosition().perp();
+ float Trk2ML1Z = tracks.at(tk2).getML1seg().globalPosition().z();
+ float Trk2ML2R = tracks.at(tk2).getML2seg().globalPosition().perp();
+ float Trk2ML2Z = tracks.at(tk2).getML2seg().globalPosition().z();
+
+ //find the distance between the tracks
+ float DistML1(1000),DistML2(1000);
+ if(tracks.at(tk1).mdtChamber() <= 11 || tracks.at(tk1).mdtChamber() == 52) {
+ DistML1 = fabs(Trk1ML1Z - Trk2ML1Z);
+ DistML2 = fabs(Trk1ML2Z - Trk2ML2Z);
+ }
+ else if(tracks.at(tk1).mdtChamber() <= 21 || tracks.at(tk1).mdtChamber() == 49) {
+ DistML1 = fabs(Trk1ML1R - Trk2ML1R);
+ DistML2 = fabs(Trk1ML2R - Trk2ML2R);
+ }
+ if(DistML1 < 40 || DistML2 < 40) {
+ //find how many MDTs the tracks share
+ std::vector<Muon::MdtPrepData*> mdt1 = tracks.at(tk1).mdtHitsOnTrack();
+ std::vector<Muon::MdtPrepData*> mdt2 = tracks.at(tk2).mdtHitsOnTrack();
+ nShared = 0;
+ for(unsigned int m1=0; m1<mdt1.size(); ++m1) {
+ for(unsigned int m2=0; m2<mdt2.size(); ++m2) {
+ if(mdt1.at(m1)->identify() == mdt2.at(m2)->identify()) {
+ nShared++;
+ break;
+ }
+ }
+ }
+
+ if(nShared <= 1) continue;//if the tracks share only 1 hits move to next track
+ //store the track as ambiguous
+ AmbigTracks.push_back(tracks.at(tk2));
+ AmbigTrks.push_back(tk2);
+ }
+ }//end chamber match
+ }//end tk2 loop
+
+ if(AmbigTracks.size() == 1) {
+ UniqueTracks.push_back(tracks.at(tk1));
+ continue;
+ }
+ if(tracks.at(tk1).mdtChamber() <= 11 || tracks.at(tk1).mdtChamber() == 52) nBarrelAmbiguousTrks += (AmbigTracks.size() - 1) ;
+ else nEndcapAmbiguousTrks += (AmbigTracks.size() - 1) ;
+ //Deal with any ambiguities
+ //Barrel tracks
+ if(tracks.at(tk1).mdtChamber() <= 11 || tracks.at(tk1).mdtChamber() == 52) {
+ bool hasMomentum = true;
+ if(tracks.at(tk1).charge() == 0) hasMomentum = false;
+ float aveX(0),aveY(0),aveZ(0),aveAlpha(0);
+ float aveP(0),nAmbigP(0),TrkCharge(tracks.at(tk1).charge());
+ bool allSameSign(true);
+ for(unsigned int i=0; i<AmbigTracks.size(); ++i) {
+ if(!hasMomentum) {
+ aveX += AmbigTracks.at(i).globalPosition().x();
+ aveY += AmbigTracks.at(i).globalPosition().y();
+ aveZ += AmbigTracks.at(i).globalPosition().z();
+ aveAlpha += AmbigTracks.at(i).getML1seg().alpha();
+ }
+ else {
+ //check the charge is the same
+ if(fabs(AmbigTracks.at(i).charge()-TrkCharge) > 0.1) allSameSign = false;
+ //find the average momentum
+ aveP += AmbigTracks.at(i).momentum().mag();
+ nAmbigP++;
+ aveAlpha += AmbigTracks.at(i).alpha();
+ aveX += AmbigTracks.at(i).globalPosition().x();
+ aveY += AmbigTracks.at(i).globalPosition().y();
+ aveZ += AmbigTracks.at(i).globalPosition().z();
+ }
+ }//end loop on ambiguous tracks
+ if(!hasMomentum) {
+ aveX = aveX/(float)AmbigTracks.size();
+ aveY = aveY/(float)AmbigTracks.size();
+ aveZ = aveZ/(float)AmbigTracks.size();
+ Amg::Vector3D gpos(aveX,aveY,aveZ);
+ aveAlpha = aveAlpha/(float)AmbigTracks.size();
+ float alphaErr = tracks.at(tk1).getML1seg().alphaError();
+ float rErr = tracks.at(tk1).getML1seg().rError();
+ float zErr = tracks.at(tk1).getML1seg().zError();
+ int Chamber = tracks.at(tk1).mdtChamber();
+ int ChEta = tracks.at(tk1).mdtChEta();
+ int ChPhi = tracks.at(tk1).mdtChPhi();
+ //
+ TrackletSegment aveSegML1(Chamber,ChEta,ChPhi,tracks.at(tk1).getML1seg().getChMidPoint(),aveAlpha,alphaErr,gpos,rErr,zErr,tracks.at(tk1).getML1seg().mdtHitsOnTrack(),0);
+ float pT = 10000.0*sin(aveSegML1.alpha());
+ float pz = 10000.0*cos(aveSegML1.alpha());
+ Amg::Vector3D momentum(pT*cos(aveSegML1.globalPosition().phi()),pT*sin(aveSegML1.globalPosition().phi()),pz);
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(tracks.at(tk1).getML1seg().rError());//delta R
+ matrix(1,1) = sq(tracks.at(tk1).getML1seg().zError());//delta z
+ matrix(2,2) = sq(0.0000001);//delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
+ matrix(3,3) = sq(tracks.at(tk1).getML1seg().alphaError());//delta theta
+ matrix(4,4) = sq(0.00005);//delta 1/p
+ Tracklet aveTrack(aveSegML1,momentum,matrix,0);
+ UniqueTracks.push_back(aveTrack);
+ }
+ else if(allSameSign) {
+ aveP = aveP/nAmbigP;
+ float pT = aveP*sin(tracks.at(tk1).getML1seg().alpha());
+ float pz = aveP*cos(tracks.at(tk1).getML1seg().alpha());
+ Amg::Vector3D momentum(pT*cos(tracks.at(tk1).globalPosition().phi()),pT*sin(tracks.at(tk1).globalPosition().phi()),pz);
+ Tracklet MyTrack = tracks.at(tk1);
+ MyTrack.momentum(momentum);
+ MyTrack.charge(tracks.at(tk1).charge());
+ UniqueTracks.push_back(MyTrack);
+ }
+ else {
+ aveX = aveX/(float)AmbigTracks.size();
+ aveY = aveY/(float)AmbigTracks.size();
+ aveZ = aveZ/(float)AmbigTracks.size();
+ Amg::Vector3D gpos(aveX,aveY,aveZ);
+ aveAlpha = aveAlpha/(float)AmbigTracks.size();
+ float alphaErr = tracks.at(tk1).getML1seg().alphaError();
+ float rErr = tracks.at(tk1).getML1seg().rError();
+ float zErr = tracks.at(tk1).getML1seg().zError();
+ int Chamber = tracks.at(tk1).mdtChamber();
+ int ChEta = tracks.at(tk1).mdtChEta();
+ int ChPhi = tracks.at(tk1).mdtChPhi();
+ TrackletSegment aveSegML1(Chamber,ChEta,ChPhi,tracks.at(tk1).getML1seg().getChMidPoint(),aveAlpha,alphaErr,gpos,rErr,zErr,tracks.at(tk1).getML1seg().mdtHitsOnTrack(),0);
+ float pT = 10000.0*sin(aveSegML1.alpha());
+ float pz = 10000.0*cos(aveSegML1.alpha());
+ Amg::Vector3D momentum(pT*cos(aveSegML1.globalPosition().phi()),pT*sin(aveSegML1.globalPosition().phi()),pz);
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(tracks.at(tk1).getML1seg().rError());//delta R
+ matrix(1,1) = sq(tracks.at(tk1).getML1seg().zError());//delta z
+ matrix(2,2) = sq(0.0000001);//delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
+ matrix(3,3) = sq(tracks.at(tk1).getML1seg().alphaError());//delta theta
+ matrix(4,4) = sq(0.00005);//delta 1/p
+ Tracklet aveTrack(aveSegML1,momentum,matrix,0);
+ UniqueTracks.push_back(aveTrack);
+ }
+ }//end barrel Tracks
+
+ //Endcap tracks
+ else if( (tracks.at(tk1).mdtChamber() > 11 && tracks.at(tk1).mdtChamber() <= 21) || tracks.at(tk1).mdtChamber() == 49 ) {
+ std::vector<Muon::MdtPrepData*> AllMdts;
+ for(unsigned int i=0; i<AmbigTracks.size(); ++i) {
+ std::vector<Muon::MdtPrepData*> mdts = AmbigTracks.at(i).mdtHitsOnTrack();
+ std::vector<Muon::MdtPrepData*> tmpAllMdt = AllMdts;
+ for(unsigned int m1=0; m1<mdts.size(); ++m1) {
+ bool isNewHit = true;
+ for(unsigned int m2=0; m2<tmpAllMdt.size(); ++m2) {
+ if(mdts.at(m1)->identify() == tmpAllMdt.at(m2)->identify()) {
+ isNewHit = false;
+ break;
+ }
+ }
+ if(isNewHit) AllMdts.push_back(mdts.at(m1));
+ }//end loop on mdts
+ }//end loop on ambiguous tracks
+
+ std::vector<TrackletSegment> MyECsegs = TrackletSegmentFitter(AllMdts);
+ if(MyECsegs.size() > 0) {
+ TrackletSegment ECseg = MyECsegs.at(0);
+ ECseg.clearMdt();
+ float pT = 10000.0*sin(ECseg.alpha());
+ float pz = 10000.0*cos(ECseg.alpha());
+ Amg::Vector3D momentum(pT*cos(ECseg.globalPosition().phi()),pT*sin(ECseg.globalPosition().phi()),pz);
+ AmgSymMatrix(5) matrix;
+ matrix.setIdentity();
+ matrix(0,0) = sq(ECseg.rError());//delta R
+ matrix(1,1) = sq(ECseg.zError());//delta z
+ matrix(2,2) = sq(0.0000001);//delta phi (~0 because we explicitly rotate all tracks into the middle of the chamber)
+ matrix(3,3) = sq(ECseg.alphaError());//delta theta
+ matrix(4,4) = sq(0.00005);//delta 1/p (endcap tracks are straight lines with no momentum that we can measure ...)
+ Tracklet MyCombTrack(MyECsegs.at(0),ECseg,momentum,matrix,0);
+ UniqueTracks.push_back(MyCombTrack);
+ }
+ else UniqueTracks.push_back(tracks.at(tk1));
+ }//end endcap tracks
+
+ }//end loop on tracks -- tk1
+
+ return UniqueTracks;
+ }
+
+}
+
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.h b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.h
new file mode 100755
index 0000000000000..b072ce1e4de07
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/MSVertexTrackletTool.h
@@ -0,0 +1,88 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MSVERTEXTRACKLETTOOL_H
+#define MSVERTEXTRACKLETTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "MSVertexToolInterfaces/IMSVertexTrackletTool.h"
+#include "MuonIdHelpers/MuonIdHelperTool.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "GeoAdaptors/GeoMuonHits.h"
+#include "Identifier/Identifier.h"
+#include "MSVertexUtils/TrackletSegment.h"
+#include "MSVertexUtils/Tracklet.h"
+#include <utility>
+#include <vector>
+
+
+namespace Muon {
+
+ class MSVertexTrackletTool : virtual public IMSVertexTrackletTool, public AthAlgTool
+ {
+
+ public :
+ /** default constructor */
+ MSVertexTrackletTool (const std::string& type, const std::string& name, const IInterface* parent);
+ /** destructor */
+ virtual ~MSVertexTrackletTool();
+
+ static const InterfaceID& interfaceID();
+
+ virtual StatusCode initialize(void);
+ virtual StatusCode finalize(void);
+
+ private:
+ //tool handles & private data members
+
+ const MdtIdHelper* m_mdtIdHelper;
+ static const MdtIdHelper* mdtCompareIdHelper;
+
+ std::string m_TPContainer;
+ float m_SeedResidual;
+ float m_minSegFinderChi2;
+ float m_BarrelDeltaAlphaCut;
+ float m_maxDeltabCut;
+ float m_EndcapDeltaAlphaCut;
+ float m_DeltaAlphaCut;
+ float m_DeltabCut;
+ float m_TrackPhiAngle;
+
+ bool m_tightTrackletRequirement;
+
+ int nMDT;
+ float PI;
+ float k_BIL;
+ float k_BML;
+ float k_BMS;
+ float k_BOL;
+
+
+ public:
+ StatusCode findTracklets(std::vector<Tracklet>& traklets);
+
+ private:
+ //private functions
+ int SortMDThits(std::vector<std::vector<Muon::MdtPrepData*> >& SortedMdt);
+ bool SortMDT(Identifier& i1, Identifier& i2);
+ std::vector<TrackletSegment> TrackletSegmentFitter(std::vector<Muon::MdtPrepData*>& mdts);
+ std::vector<TrackletSegment> TrackletSegmentFitterCore(std::vector<Muon::MdtPrepData*>& mdts,std::vector<std::pair<float,float> >& SeedParams);
+ std::vector<std::pair<float,float> > SegSeeds(std::vector<Muon::MdtPrepData*>& mdts);
+ float SeedResiduals(std::vector<Muon::MdtPrepData*>& mdts, float slope, float inter);
+ std::vector<TrackletSegment> CleanSegments(std::vector<TrackletSegment>& segs);
+ bool DeltabCalc(TrackletSegment& ML1seg, TrackletSegment& ML2seg);
+ float TrackMomentum(int chamber,float deltaAlpha);
+ float TrackMomentumError(TrackletSegment& ml1, TrackletSegment& ml2);
+ float TrackMomentumError(TrackletSegment& ml1);
+ std::vector<Tracklet> ResolveAmbiguousTracklets(std::vector<Tracklet>& tracks);
+ void convertToTrackParticles(std::vector<Tracklet>& tracklets);
+ float sq(float x) { return (x)*(x); }
+ static bool mdtComp(const Muon::MdtPrepData* mprd1, const Muon::MdtPrepData* mprd2);
+
+ };
+
+
+}
+#endif
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_entries.cxx b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_entries.cxx
new file mode 100644
index 0000000000000..534bae88dc449
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_entries.cxx
@@ -0,0 +1,19 @@
+
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "../MSVertexTrackletTool.h"
+#include "../MSVertexRecoTool.h"
+
+using namespace Muon;
+
+DECLARE_TOOL_FACTORY( MSVertexTrackletTool )
+DECLARE_TOOL_FACTORY( MSVertexRecoTool )
+
+DECLARE_FACTORY_ENTRIES( MSVertexTools )
+{
+ DECLARE_TOOL( MSVertexTrackletTool )
+ DECLARE_TOOL( MSVertexRecoTool )
+}
+
+
+
+
diff --git a/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_load.cxx b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_load.cxx
new file mode 100644
index 0000000000000..2d776463a2c5d
--- /dev/null
+++ b/MuonSpectrometer/MSVertexReconstruction/MSVertexTools/src/components/MSVertexTools_load.cxx
@@ -0,0 +1,5 @@
+
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( MSVertexTools )
+
--
GitLab