diff --git a/Trigger/TrigFTK/FTK_DataProviderSvc/FTK_DataProviderSvc/FTK_DataProviderSvc.h b/Trigger/TrigFTK/FTK_DataProviderSvc/FTK_DataProviderSvc/FTK_DataProviderSvc.h
index 8bc315c3c9d040c64ec922ff68fd084e0ec9a0f4..c420638fb8c5a4e82fad9c679d2a722ab21f0e38 100644
--- a/Trigger/TrigFTK/FTK_DataProviderSvc/FTK_DataProviderSvc/FTK_DataProviderSvc.h
+++ b/Trigger/TrigFTK/FTK_DataProviderSvc/FTK_DataProviderSvc/FTK_DataProviderSvc.h
@@ -32,6 +32,7 @@
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/TrackParticleContainer.h"
#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
+#include "FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h"
/// Forward Declarations ///
class AtlasDetectorID;
@@ -157,14 +158,16 @@ class FTK_DataProviderSvc : public virtual IFTK_DataProviderSvc, virtual public
ToolHandle< InDet::IVertexFinder > m_VertexFinderTool;
ToolHandle< IFTK_VertexFinderTool > m_RawVertexFinderTool;
ToolHandle< Trk::IRIO_OnTrackCreator > m_ROTcreator;
+ ToolHandle< IFTK_DuplicateTrackRemovalTool > m_DuplicateTrackRemovalTool;
double m_trainingBeamspotX;
double m_trainingBeamspotY;
double m_trainingBeamspotZ;
double m_trainingBeamspotTiltX;
double m_trainingBeamspotTiltY;
-
-
+
+ bool m_remove_duplicates;
+
const FTK_RawTrackContainer* m_ftk_tracks;
// Track Cache
diff --git a/Trigger/TrigFTK/FTK_DataProviderSvc/src/FTK_DataProviderSvc.cxx b/Trigger/TrigFTK/FTK_DataProviderSvc/src/FTK_DataProviderSvc.cxx
index 6d0225b4adc6b6df150918f077d0398c60e52e48..e0b3bae6c0f7304052c016ef5393deae4abfbed8 100644
--- a/Trigger/TrigFTK/FTK_DataProviderSvc/src/FTK_DataProviderSvc.cxx
+++ b/Trigger/TrigFTK/FTK_DataProviderSvc/src/FTK_DataProviderSvc.cxx
@@ -53,6 +53,7 @@
#include "PixelConditionsServices/IPixelOfflineCalibSvc.h"
#include "TrkToolInterfaces/ITrackSummaryTool.h"
#include "TrkTrackSummary/TrackSummary.h"
+#include "FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h"
#include "FTK_RecToolInterfaces/IFTK_VertexFinderTool.h"
#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
@@ -102,11 +103,13 @@ FTK_DataProviderSvc::FTK_DataProviderSvc(const std::string& name, ISvcLocator* s
m_VertexFinderTool("InDet::InDetIterativePriVxFinderTool"),
m_RawVertexFinderTool("FTK_VertexFinderTool"),
m_ROTcreator("Trk::IRIO_OnTrackCreator/FTK_ROTcreatorTool"),
+ m_DuplicateTrackRemovalTool("FTK_DuplicateTrackRemovalTool"),
m_trainingBeamspotX(0.),
m_trainingBeamspotY(0.),
m_trainingBeamspotZ(0.),
m_trainingBeamspotTiltX(0.),
m_trainingBeamspotTiltY(0.),
+ m_remove_duplicates(false),
m_ftk_tracks(nullptr),
m_conv_tracks(nullptr),
m_refit_tracks(nullptr),
@@ -167,6 +170,7 @@ FTK_DataProviderSvc::FTK_DataProviderSvc(const std::string& name, ISvcLocator* s
declareProperty("TrackFitter", m_trackFitter);
declareProperty("UncertaintyTool",m_uncertaintyTool);
declareProperty("TrackSummaryTool", m_trackSumTool);
+ declareProperty("DuplicateTrackRemovalTool",m_DuplicateTrackRemovalTool);
declareProperty("TrackParticleCreatorTool", m_particleCreatorTool);
declareProperty("VertexFinderTool",m_VertexFinderTool);
declareProperty("ROTcreatorTool",m_ROTcreator);
@@ -194,7 +198,7 @@ FTK_DataProviderSvc::FTK_DataProviderSvc(const std::string& name, ISvcLocator* s
declareProperty("CorrectSCTClusters",m_correctSCTClusters);
declareProperty("setBroadPixelClusterOnTrackErrors",m_broadPixelErrors);
declareProperty("setBroadSCT_ClusterOnTrackErrors",m_broadSCT_Errors);
-
+ declareProperty("RemoveDuplicates",m_remove_duplicates);
}
@@ -240,6 +244,8 @@ StatusCode FTK_DataProviderSvc::initialize() {
ATH_CHECK(m_particleCreatorTool.retrieve());
ATH_MSG_INFO( " getting vertexFinderTool tool with name " << m_VertexFinderTool.name());
ATH_CHECK(m_VertexFinderTool.retrieve());
+ ATH_MSG_INFO( " getting DuplicateTrackRemovalTool tool with name " << m_DuplicateTrackRemovalTool.name());
+ ATH_CHECK(m_DuplicateTrackRemovalTool.retrieve());
ATH_MSG_INFO( " getting FTK_RawTrackVertexFinderTool tool with name " << m_RawVertexFinderTool.name());
ATH_CHECK(m_RawVertexFinderTool.retrieve());
ATH_MSG_INFO( " getting ROTcreator tool with name " << m_ROTcreator.name());
@@ -1002,23 +1008,33 @@ void FTK_DataProviderSvc::getFTK_RawTracksFromSG(){
}
// new event - get the tracks from StoreGate
-
if (!m_storeGate->contains<FTK_RawTrackContainer>(m_RDO_key)) {
ATH_MSG_DEBUG( "getFTK_RawTracksFromSG: FTK tracks "<< m_RDO_key <<" not found in StoreGate !");
return;
- } else {
- ATH_MSG_VERBOSE( "getFTK_RawTracksFromSG: Doing storegate retreive");
- StatusCode sc = m_storeGate->retrieve(m_ftk_tracks, m_RDO_key);
- if (sc.isFailure()) {
- ATH_MSG_VERBOSE( "getFTK_RawTracksFromSG: Failed to get FTK Tracks Container");
- return;
+ } else {
+ if (m_remove_duplicates){//get all tracks, and then call duplicate rmoval tool
+ ATH_MSG_INFO( "getFTK_RawTracksFromSG: Doing storegate retreive and then duplicate removal");
+ const FTK_RawTrackContainer* temporaryTracks=nullptr;
+ StatusCode sc = m_storeGate->retrieve(temporaryTracks, m_RDO_key);
+ m_ftk_tracks = m_DuplicateTrackRemovalTool->removeDuplicates(temporaryTracks);
+ if (sc.isFailure()) {
+ ATH_MSG_WARNING( "getFTK_RawTracksFromSG: Failed to get FTK Tracks Container when using removeDumplicates ");
+ return;
+ }
}
- ATH_MSG_DEBUG( "getFTK_RawTracksFromSG: Got " << m_ftk_tracks->size() << " raw FTK tracks (RDO) from StoreGate ");
+ else{//the original way
+ ATH_MSG_VERBOSE( "getFTK_RawTracksFromSG: Doing storegate retreive");
+ StatusCode sc = m_storeGate->retrieve(m_ftk_tracks, m_RDO_key);
+ if (sc.isFailure()) {
+ ATH_MSG_WARNING( "getFTK_RawTracksFromSG: Failed to get FTK Tracks Container");
+ return;
+ }
+ }
+ ATH_MSG_INFO( "getFTK_RawTracksFromSG: Got " << m_ftk_tracks->size() << " raw FTK tracks (RDO) from StoreGate ");
m_gotRawTracks = true;
-
}
+
// Creating collection for pixel clusters
-
m_PixelClusterContainer = new InDet::PixelClusterContainer(m_pixelId->wafer_hash_max());
m_PixelClusterContainer->addRef();
StatusCode sc = m_storeGate->record(m_PixelClusterContainer,m_PixelClusterContainerName);
diff --git a/Trigger/TrigFTK/FTK_RecToolInterfaces/FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h b/Trigger/TrigFTK/FTK_RecToolInterfaces/FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..a182148a6e478ebb35937dfe74ae966106f06e8c
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecToolInterfaces/FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h
@@ -0,0 +1,35 @@
+
+//abstract interface
+
+#ifndef __IFTK_DUPLICATETRACKREMOVAL_TOOL_H__
+#define __IFTK_DUPLICATETRACKREMOVAL_TOOL_H__
+
+#include "GaudiKernel/IAlgTool.h"
+//#include "AthenaKernel/IOVSvcDefs.h"
+class VxContainer;
+#include "TrigFTK_RawData/FTK_RawTrack.h"
+#include "TrigFTK_RawData/FTK_RawTrackContainer.h"
+#include "xAODTracking/VertexFwd.h"
+//#include "xAODTracking/TrackParticleFwd.h"
+#include "xAODTracking/VertexContainerFwd.h"
+#include "xAODTracking/VertexAuxContainer.h"
+#include "TrkTrack/TrackCollection.h"
+//#include "xAODTracking/TrackParticleContainerFwd.h"
+//namespace Trk {
+//// class Track;
+//}
+static const InterfaceID IID_IFTK_DuplicateTrackRemovalTool("IFTK_DuplicateTrackRemovalTool",1,0);
+
+class IFTK_DuplicateTrackRemovalTool : virtual public IAlgTool {
+
+ public:
+ /** other standard AlgTool methods */
+ static const InterfaceID& interfaceID () //!< the Tool's interface
+ { return IID_IFTK_DuplicateTrackRemovalTool; }
+
+ virtual FTK_RawTrackContainer* removeDuplicates(const FTK_RawTrackContainer* trks) = 0;
+
+ private:
+};
+
+#endif
diff --git a/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_DuplicateTrackRemovalTool.h b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_DuplicateTrackRemovalTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..4f8c5994152d6df061abd6d46d396743aad61cfc
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_DuplicateTrackRemovalTool.h
@@ -0,0 +1,122 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+////////////////////////////////////////////////////////////////////////////////
+// FTK_DuplicateTrackRemovalTool tool
+// -------------------------------
+// ATLAS Collaboration
+//
+// Remove duplicate (overlapping) tracks
+//
+// June 2017: Tool created
+//
+// Author: Andy Haas, NYU
+// e-mail: ahaas@cern.ch
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef __TRIG_FTK_DUPLICATETRACKREMOVAL_TOOL_H__
+#define __TRIG_FTK_DUPLICATETRACKREMOVAL_TOOL_H__
+
+#include "GaudiKernel/ToolHandle.h"
+
+#include "AthenaBaseComps/AthAlgTool.h"
+
+#include "FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h"
+#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
+#include "TrkVxEdmCnv/IVxCandidateXAODVertex.h"
+
+#include "TrigFTK_RawData/FTK_RawTrack.h"
+#include "TrigFTK_RawData/FTK_RawTrackContainer.h"
+#include "TrkTrack/TrackCollection.h"
+#include "VxVertex/VxContainer.h"
+//#include "xAODTracking/Vertex.h"
+//#include "xAODTracking/TrackParticle.h"
+//#include "xAODTracking/VertexContainer.h"
+//#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/VertexFwd.h"
+//#include "xAODTracking/TrackParticleFwd.h"
+#include "xAODTracking/VertexContainerFwd.h"
+#include "xAODTracking/VertexAuxContainer.h"
+//#include "xAODTracking/TrackParticleContainerFwd.h"
+//#include "Tracking/TrkVertexFitter/TrkVxEdmCnv/TrkVxEdmCnv/IVxCandidateXAODVertex.h"
+class VxContainer;
+using std::vector;
+namespace Trk {
+ class Track;
+ //class Trk::IVxCandidateXAODVertex;
+}
+class FTK_DuplicateTrackRemovalTool : public AthAlgTool, virtual public IFTK_DuplicateTrackRemovalTool
+{
+ //struct to hold track parameters
+ struct MyTrack{
+ MyTrack(int idx,float pt,float theta,float phi,float d0,float z0,float pterr,float thetaerr,float phierr,float d0err,float z0err){
+ m_idx=idx;
+ m_pt=pt;
+ m_theta=theta;
+ m_phi=phi;
+ m_d0=d0;
+ m_z0=z0;
+
+ m_pterr=pterr;
+ // m_qoverperr=qoverperr;
+ m_thetaerr=thetaerr;
+ m_phierr=phierr;
+ m_d0err=d0err;
+ m_z0err=z0err;
+ }
+ int m_idx;
+ double m_pt;
+ double m_theta;
+ double m_phi;
+ double m_d0;
+ double m_z0;
+
+ double m_pterr;
+ //double m_qoverperr=qoverperr;
+ double m_thetaerr;
+ double m_phierr;
+ double m_d0err;
+ double m_z0err;
+
+ //Track pT sort
+ bool operator<(const MyTrack &a)const{
+ return fabs(m_pt) > fabs(a.m_pt);
+ }
+ };
+
+ public:
+
+ FTK_DuplicateTrackRemovalTool( const std::string&, const std::string&, const IInterface* );
+ virtual ~FTK_DuplicateTrackRemovalTool(){};
+
+ virtual StatusCode initialize();
+ virtual StatusCode finalize ();
+
+ FTK_RawTrackContainer* removeDuplicates(const FTK_RawTrackContainer* trks);
+
+ private:
+
+ bool m_barrelOnly;
+ bool m_hasIBL;
+ double m_cluster_size;
+ double m_chi2cut;
+ double m_constTrkPt;
+ double m_constTrkEta;
+ double m_z0errfactor;
+
+ FTK_RawTrackContainer* m_trks_nodups;
+
+ std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> findVertex(vector<MyTrack> trks);
+ double ctheta2eta(double cot);
+ vector<MyTrack> getTracks(const FTK_RawTrackContainer* trks);
+ vector<MyTrack> getTracks(const TrackCollection* trks);
+
+ ToolHandle<IFTK_UncertaintyTool> m_uncertaintyTool;
+};
+
+
+#endif
+
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/FTK_DuplicateTrackRemovalTool.cxx b/Trigger/TrigFTK/FTK_RecTools/src/FTK_DuplicateTrackRemovalTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a6e6163215a10026d700605cb549e2d911c8035f
--- /dev/null
+++ b/Trigger/TrigFTK/FTK_RecTools/src/FTK_DuplicateTrackRemovalTool.cxx
@@ -0,0 +1,419 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <cmath>
+#include <iostream>
+#include <vector>
+#include "TVector3.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "FTK_RecTools/FTK_DuplicateTrackRemovalTool.h"
+#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
+#include "TrigFTK_RawData/FTK_RawTrack.h"
+#include "TrkTrack/TrackCollection.h"
+#include "VxVertex/VxContainer.h"
+#include "VxVertex/ExtendedVxCandidate.h"
+#include "VxVertex/RecVertex.h"
+#include "VxVertex/Vertex.h"
+#include "VxVertex/VxTrackAtVertex.h"
+#include "TrigFTK_RawData/FTK_RawTrackContainer.h"
+#include "TrkParameters/TrackParameters.h"
+#include "TrkTrack/Track.h"
+
+#include "xAODTracking/TrackParticleAuxContainer.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/VertexAuxContainer.h"
+#include "VxVertex/VxCandidate.h"
+#include "TrkVxEdmCnv/IVxCandidateXAODVertex.h"
+#include <map>
+#include <vector>
+#include <utility>
+//#include "TrkSurfaces/PlaneSurface.h"
+
+using std::map; using std::string;
+using std::cout; using std::endl;
+using std::vector; using std::iterator;
+using Trk::RecVertex;using Trk::VxTrackAtVertex;
+
+FTK_DuplicateTrackRemovalTool::FTK_DuplicateTrackRemovalTool(const std::string& t,
+ const std::string& n,
+ const IInterface* p ):
+ AthAlgTool(t,n,p),
+ m_barrelOnly(false),
+ m_hasIBL(false),
+ m_cluster_size(1.87),
+ m_chi2cut(5.),
+ m_constTrkPt(1.),
+ m_constTrkEta(1.1),
+ m_z0errfactor(1.0),
+ m_uncertaintyTool("FTK_UncertaintyTool",this)
+{
+ declareInterface< IFTK_DuplicateTrackRemovalTool >( this );
+ declareProperty( "HasIBL", m_hasIBL);
+ declareProperty( "BarrelOnly", m_barrelOnly);
+ declareProperty( "Cluster_size", m_cluster_size);
+ declareProperty( "Chi2cut", m_chi2cut);
+ declareProperty( "ConstTrkPt", m_constTrkPt);
+ declareProperty( "ConstTrkEta", m_constTrkEta);
+ declareProperty( "UncertaintyTool", m_uncertaintyTool);
+ declareProperty( "z0ErrFactor", m_z0errfactor);
+}
+
+StatusCode FTK_DuplicateTrackRemovalTool::initialize() {
+
+ StatusCode sc = AlgTool::initialize();
+ MsgStream athenaLog(msgSvc(), name());
+
+ m_trks_nodups = new FTK_RawTrackContainer; //we DO own the tracks
+ //(SG::VIEW_ELEMENTS);//we don't own the tracks, we're just going to hold them
+
+ athenaLog << MSG::INFO << "FTK_DuplicateTrackRemovalTool initialized "<< endmsg;
+ return sc;
+}
+
+StatusCode FTK_DuplicateTrackRemovalTool::finalize() {
+ StatusCode sc = AlgTool::finalize();
+ return sc;
+}
+
+FTK_RawTrackContainer* FTK_DuplicateTrackRemovalTool::removeDuplicates(const FTK_RawTrackContainer* trks){
+ ATH_MSG_INFO("ACH99 - I'm in removeDuplicates!");
+ m_trks_nodups->clear();
+ m_trks_nodups->reserve(trks->size());
+ for (unsigned int i = 0; i!=trks->size(); i++) m_trks_nodups->push_back(new FTK_RawTrack(*(trks->at(i))));
+ return m_trks_nodups;
+}
+
+std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> FTK_DuplicateTrackRemovalTool::findVertex(vector<MyTrack> mytrk)
+{
+
+ MsgStream athenaLog(msgSvc(), name());
+ double mGkx=0.0001;
+ double mGky=0.0001;
+ double mGkz=16.;
+ athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: barrel " << m_barrelOnly<< endmsg;
+ xAOD::VertexContainer* theVertexContainer = new xAOD::VertexContainer;
+ xAOD::VertexAuxContainer* theVertexAuxContainer = new xAOD::VertexAuxContainer;
+ theVertexContainer->setStore( theVertexAuxContainer );
+ vector<MyTrack>::iterator trkbegin;
+ vector<MyTrack>::iterator trkend;
+ vector<MyTrack> trkatvtx;
+
+ athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: total " << mytrk.size()<< " tracks "<< endmsg;
+ sort(mytrk.begin(),mytrk.end());//sorting mytrk by pt(Track7.h)
+
+ int numdo=0;
+ do{
+ athenaLog << MSG::DEBUG << "findVertex: finding vertex "<< endmsg;
+ //////////////////////
+ //clustering
+ //////////////////////
+ numdo++;
+ vector<MyTrack> vxtrk;
+ //high pt track selecting
+ trkbegin=mytrk.begin();
+ trkend=mytrk.end();
+ float seed_z=(*trkbegin).m_z0;
+
+ //selecting tracks near high pt track
+ for(vector<MyTrack>::iterator i=trkbegin ; i<trkend;){
+ if(fabs((*i).m_z0 - seed_z) < m_cluster_size){
+ vxtrk.push_back(*i);
+ i=mytrk.erase(i);
+ }
+ else{
+ i++;
+ }
+ trkbegin=mytrk.begin();
+ trkend=mytrk.end();
+ }
+
+ //calculating seed z position(z average weighted by z error)
+ double oldz=0;
+ double z_weight=0;
+ trkbegin=vxtrk.begin();
+ trkend=vxtrk.end();
+ for(vector<MyTrack>::iterator i=trkbegin ; i<trkend;i++){
+ oldz += (*i).m_z0*(*i).m_z0err;
+ z_weight += (*i).m_z0err;
+ }
+ oldz /= z_weight;
+
+ double oldx=0;
+ double oldy=0;
+ double m_Gk[3][3]={{mGkx,0.,0.},{0.,mGky,0.},{0.,0.,mGkz}};//error
+ Amg::MatrixX C22_mat(3,3);
+ double chi2=0;
+
+ trkbegin=vxtrk.begin();
+ trkend=vxtrk.end();
+ int trknumber=0;
+ athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: fit VTX, using "<<vxtrk.size()<<" trks."<< endmsg;
+ for( vector<MyTrack>::iterator i=trkbegin ; i<trkend ;){
+ trknumber++;
+ double tmpchi2 = chi2;
+ athenaLog << MSG::VERBOSE << "getTracks: "<<trknumber<<", pt "<<(*i).m_pt<<", eta "<<(*i).m_theta<<", phi "<<(*i).m_phi<<", d0 "<<(*i).m_d0<<", z0 "<<(*i).m_z0<<", pt err "<<(*i).m_pterr<<", theta err "<<(*i).m_thetaerr<<", phi err "<< (*i).m_phierr<<", d0 err"<< (*i).m_d0err<<", z0 err "<<(*i).m_z0err<< endmsg;
+ //track parameter reading
+ double xv,yv,zv,P0,phi0,theta0;
+ TVector3 trk;
+ xv=oldx;
+ yv=oldy;
+ zv=oldz;
+ P0=(*i).m_pt;
+ phi0=(*i).m_phi;
+ theta0=(*i).m_theta;
+
+ //parameter of angle
+
+ double cosPhi0,sinPhi0,sinPsi,cosPsi;
+ double psi,ctt,sint;
+ double alpha=0.02997*20.84;///1000.0;
+
+ cosPhi0=cos(phi0);sinPhi0=sin(phi0);
+ sinPsi=-alpha*(xv*cosPhi0+yv*sinPhi0)/P0;
+ cosPsi=sqrt(1.0-sinPsi*sinPsi);
+ psi=asin(sinPsi);
+ sint=sin(theta0);
+ ctt=cos(theta0)/sint;
+
+ //difference of values expected and observed
+
+ double m_A[2][3],m_B[2][3],m_u[2],m_h[2],m_resid[2];
+
+ m_A[0][0]=-sin(phi0+psi)/cosPsi;
+ m_A[0][1]= cos(phi0+psi)/cosPsi;
+ m_A[0][2]=0.0;
+
+ m_A[1][0]=-ctt*cosPhi0/cosPsi;
+ m_A[1][1]=-ctt*sinPhi0/cosPsi;
+ m_A[1][2]=1.0;
+
+ m_B[0][0]=-xv*m_A[0][1]+yv*m_A[0][0];
+ m_B[0][1]=0.0;
+ m_B[0][2]=(1.0-1.0/cosPsi)/alpha;
+
+ m_B[1][0]=-xv*m_A[1][1]+yv*m_A[1][0];
+ m_B[1][1]=-P0*psi/(alpha*sint*sint);
+ m_B[1][2]=ctt*(psi-sinPsi/cosPsi)/alpha;
+
+ m_u[0]=(*i).m_d0;
+ m_u[1]=(*i).m_z0;
+
+ m_h[0]=yv*cosPhi0-xv*sinPhi0+P0*(1-cosPsi)/alpha;
+ m_h[1]=zv+P0*ctt*psi/alpha;
+
+ m_resid[0]=m_u[0]-m_h[0];
+ m_resid[1]=m_u[1]-m_h[1];
+
+ //error of difference
+
+ double m_Vqq[3][3]={{(*i).m_phierr*(*i).m_phierr,0.,0.},{0.,(*i).m_thetaerr*(*i).m_thetaerr,0.},{0.,0.,(*i).m_pterr*(*i).m_pterr}};
+ double m_Vuu[2][2]={{(*i).m_d0err*(*i).m_d0err,0.},{0.,(*i).m_z0err*(*i).m_z0err}};
+ double AC[2][3],BV[2][3],Sk[2][2];
+
+ for(int i1=0;i1<2;i1++) for(int j=0;j<2;j++) Sk[i1][j]=m_Vuu[i1][j];
+ for(int i1=0;i1<2;i1++) for(int j=0;j<3;j++)
+ {
+ AC[i1][j]=0.0;
+ for(int k=0;k<3;k++) AC[i1][j]+=m_A[i1][k]*m_Gk[j][k];
+ }
+ for(int i1=0;i1<2;i1++) for(int j=0;j<2;j++)
+ {
+ for(int k=0;k<3;k++) Sk[i1][j]+=AC[i1][k]*m_A[j][k];
+ }
+ for(int i1=0;i1<2;i1++)
+ for(int j=0;j<3;j++)
+ {
+ BV[i1][j]=0.0;
+ for(int k=0;k<3;k++) BV[i1][j]+=m_B[i1][k]*m_Vqq[k][j];
+ }
+ for(int i1=0;i1<2;i1++)
+ for(int j=0;j<2;j++)
+ {
+ for(int k=0;k<3;k++) Sk[i1][j]+=BV[i1][k]*m_B[j][k];
+ }
+
+ //error determinant
+ double detr,m_V[2][2];
+
+ detr=1.0/(Sk[0][0]*Sk[1][1]-Sk[0][1]*Sk[1][0]);
+ m_V[0][0]=Sk[1][1]*detr;
+ m_V[1][1]=Sk[0][0]*detr;
+ m_V[0][1]=m_V[1][0]=-Sk[0][1]*detr;
+
+ //chi2
+ chi2=m_V[0][0]*m_resid[0]*m_resid[0]+m_V[1][1]*m_resid[1]*m_resid[1]+2.0*m_V[0][1]*m_resid[1]*m_resid[0];
+ if(chi2>m_chi2cut || chi2<0){
+ i=vxtrk.erase(i);
+ chi2 = tmpchi2;
+ }
+ else{
+ trkatvtx.push_back(*i);
+ i++;
+ }
+ trkbegin=vxtrk.begin();
+ trkend=vxtrk.end();
+ //vertex position & covariance update
+ double K[2][3],dV[3];
+
+ for(int i2=0;i2<2;i2++){
+ for(int j=0;j<3;j++){
+ for(int k=0;k<2;k++)K[i2][j]=AC[k][j]*m_V[k][i2];
+ }
+ }
+
+ for(int i2=0;i2<3;i2++){
+ dV[i2]=K[0][i2]*m_resid[0]+K[1][i2]*m_resid[1];
+ for(int j2=i2;j2<3;j2++){
+ m_Gk[i2][j2]-=K[0][i2]*AC[0][j2]+K[1][i2]*AC[1][j2];
+ m_Gk[j2][i2]=m_Gk[i2][j2];
+ }
+ }
+ oldx+=dV[0];
+ oldy+=dV[1];
+ oldz+=dV[2];
+ C22_mat(0,0)=m_Gk[0][0];C22_mat(0,1)=m_Gk[0][1];C22_mat(0,2)=m_Gk[0][0];
+ C22_mat(1,0)=m_Gk[1][0];C22_mat(1,1)=m_Gk[1][1];C22_mat(1,2)=m_Gk[1][2];
+ C22_mat(2,0)=m_Gk[2][0];C22_mat(2,1)=m_Gk[2][1];C22_mat(2,2)=m_Gk[2][2];
+ }//track loop end
+ if (vxtrk.size()==0)continue;
+ athenaLog << MSG::DEBUG << "findVertex: find vertex at "<< oldx<<";"<<oldy<<";"<<oldz<< endmsg;
+ Amg::Vector3D frameOrigin(oldx,oldy,oldz);
+ int ndf = 2*vxtrk.size()-3;
+ xAOD::Vertex* vx = new xAOD::Vertex;
+ vx->makePrivateStore();
+ vx->setPosition (frameOrigin);
+ vx->setCovariancePosition (C22_mat);
+ vx->setFitQuality(chi2,static_cast<float>(ndf));
+ vx->setVertexType(xAOD::VxType::PriVtx);
+
+ std::vector<VxTrackAtVertex> & tracksAtVertex = vx->vxTrackAtVertex(); tracksAtVertex.clear();
+ // Store the tracks at vertex
+ Amg::Vector3D Vertex(frameOrigin[0],frameOrigin[1],frameOrigin[2]);
+ const Trk::PerigeeSurface Surface(Vertex);
+ Trk::Perigee * refittedPerigee(0);
+ vector<MyTrack>::iterator i;
+ for (i = trkatvtx.begin(); i != trkatvtx.end() ; ++i)
+ {
+ AmgSymMatrix(5) * CovMtxP = new AmgSymMatrix(5);
+ double eta = ctheta2eta((*i).m_theta);
+ double etaerr=(*i).m_thetaerr*cosh(eta);
+ double qoverperr=cosh(eta)? 2./cosh(eta) * sqrt(4*(*i).m_pterr*(*i).m_pterr + 1/4/(*i).m_pt/(*i).m_pt*tanh(eta)*tanh(eta)*etaerr*etaerr) : 10;
+ (*CovMtxP)(0,0)=(*i).m_d0err*(*i).m_d0err;
+ (*CovMtxP)(0,0)=(*i).m_z0err*(*i).m_z0err;
+ (*CovMtxP)(0,0)=(*i).m_phierr*(*i).m_phierr;
+ (*CovMtxP)(0,0)=(*i).m_thetaerr*(*i).m_thetaerr;
+ (*CovMtxP)(0,0)=qoverperr*qoverperr;
+ double qoverp = 2/(*i).m_pt/cosh(eta);
+ refittedPerigee = new Trk::Perigee ((*i).m_d0,(*i).m_z0,(*i).m_phi,(*i).m_theta,qoverp, Surface, CovMtxP);
+ tracksAtVertex.emplace_back(-1, refittedPerigee);
+ }
+
+ theVertexContainer->push_back(vx);
+
+ trkatvtx.clear();
+ //delete vx;
+ }while(mytrk.size()>0);//vertex loop end
+
+ for (auto pv = theVertexContainer->begin(); pv != theVertexContainer->end(); ++pv) {
+ athenaLog << MSG::DEBUG << "vertex x = " <<(*pv)->position().x()<< endmsg;
+ }
+ return std::make_pair(theVertexContainer, theVertexAuxContainer);
+ //return theVertexContainer;
+
+}
+
+
+vector<FTK_DuplicateTrackRemovalTool::MyTrack> FTK_DuplicateTrackRemovalTool::getTracks(const FTK_RawTrackContainer* trks){
+
+ MsgStream athenaLog(msgSvc(), name());
+ vector<MyTrack> mytrk;
+ athenaLog << MSG::DEBUG << "getRawTracks: find "<< trks->size()<< " tracks "<< endmsg;
+ for (unsigned int ftk_track_index = 0; ftk_track_index != trks->size(); ++ftk_track_index){
+ const FTK_RawTrack* ftk_track= trks->at(ftk_track_index);
+ float trk_theta = std::atan2(1.0,ftk_track->getCotTh());
+ double invpt= ftk_track->getInvPt();
+ float trk_eta = -std::log(std::tan(trk_theta/2));
+ float trk_phi = ftk_track->getPhi();
+ float trk_d0 = ftk_track->getD0();
+ float trk_z0 = ftk_track->getZ0();
+ float trk_pt = 1/invpt;
+/* double trk_phierr=sqrt(3.446e-7+29.24*invpt*invpt);
+ double trk_thetaerr=sqrt(7.093e-6+38.4*invpt*invpt);
+ double trk_pterr=1;
+ double trk_d0err=sqrt(1.662e-3+6.947e4*invpt*invpt);
+ double trk_z0err=sqrt(6.472e-2+1.587e5*invpt*invpt);
+*/
+ double trk_d0err= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::d0, FTKTrackParam::d0));
+ double trk_z0err= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::z0, FTKTrackParam::z0))*m_z0errfactor;
+ double trk_phierr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::phi, FTKTrackParam::phi));
+ double trk_thetaerr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::theta, FTKTrackParam::theta));
+ double trk_pterr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::pt, FTKTrackParam::pt));
+
+ if(fabs(trk_pt)<m_constTrkPt)continue;
+ if(m_barrelOnly&&fabs(trk_eta)>m_constTrkEta)continue;
+ if(ftk_track_index<10)athenaLog << MSG::DEBUG << "getRawTracks: get pt "<< trk_pt<<"+/-"<<trk_pterr << ", eta "<<trk_eta<<"+/-"<<trk_thetaerr<<", phi"<<trk_phi<<"+/-"<<trk_phierr<<", d0"<<trk_d0<<"+/-"<<trk_d0err<<", z0"<<trk_z0<<"+/-"<<trk_z0err <<endmsg;
+ MyTrack Trk(ftk_track_index,trk_pt,trk_theta,trk_phi,trk_d0,trk_z0,trk_pterr,trk_thetaerr,trk_phierr,trk_d0err,trk_z0err);
+ mytrk.push_back(Trk);
+ }
+ athenaLog << MSG::DEBUG << "getRawTracks: total "<< mytrk.size() << " track "<< endmsg;
+ return mytrk;
+}
+
+vector<FTK_DuplicateTrackRemovalTool::MyTrack> FTK_DuplicateTrackRemovalTool::getTracks(const TrackCollection* trks){
+ MsgStream athenaLog(msgSvc(), name());
+ vector<MyTrack> mytrk;
+ athenaLog << MSG::DEBUG << "getRefitTracks: find "<< trks->size()<< " tracks "<< endmsg;
+ TrackCollection::const_iterator track_it = trks->begin();
+ TrackCollection::const_iterator last_track = trks->end();
+ for (int iTrack = 0 ; track_it!= last_track; track_it++, iTrack++) {
+// double invpt= sin(trk_Theta)*ftk_track->getCurv()/2.;
+ double trk_qOverp= (*track_it)->perigeeParameters()->parameters()[Trk::qOverP];
+ float trk_theta = (*track_it)->perigeeParameters()->parameters()[Trk::theta];
+ float trk_eta = (*track_it)->perigeeParameters()->eta();//ctheta2eta(trk_theta);
+ float trk_phi = (*track_it)->perigeeParameters()->parameters()[Trk::phi0];
+ float trk_d0 = (*track_it)->perigeeParameters()->parameters()[Trk::d0];
+ float trk_z0 = (*track_it)->perigeeParameters()->parameters()[Trk::z0];
+ float trk_pt = 1/(trk_qOverp/sin(trk_theta));
+/*
+ float invpt = 1/(trk_pt);
+ double trk_phierr=sqrt(3.446e-7+29.24*invpt*invpt);
+ double trk_thetaerr=sqrt(7.093e-6+38.4*invpt*invpt);
+ double trk_pterr=1;
+ double trk_d0err=sqrt(1.662e-3+6.947e4*invpt*invpt);
+ double trk_z0err=sqrt(6.472e-2+1.587e5*invpt*invpt)*m_z0errfactor;
+*/
+ const AmgSymMatrix(5) *cov=(*track_it)->perigeeParameters()->covariance();
+ double trk_d0err=sqrt((*cov)(0,0));
+ double trk_z0err=sqrt((*cov)(1,1));
+ double trk_phierr=sqrt((*cov)(2,2));
+ double trk_thetaerr=sqrt((*cov)(3,3));
+ double trk_qOverperr=sqrt((*cov)(4,4));
+/* double trk_d0err=(*track_it)->measurementsOnTrack()->localCovariance().error(Trk::d0);
+ double trk_z0err=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::z0);
+ double trk_phi0err=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::phi0);
+ double trk_thetaerr=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::theta);
+ double trk_qOverperr=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::qOverP);
+*/ double trk_pterr=-sin(trk_theta)*trk_qOverperr/trk_qOverp/trk_qOverp+cos(trk_theta)*trk_thetaerr/trk_qOverp;
+ if(fabs(trk_pt)<m_constTrkPt)continue;
+ if(m_barrelOnly&&fabs(trk_eta)>m_constTrkEta)continue;
+ if(iTrack<10)athenaLog << MSG::DEBUG << "getRefitTracks: get pt "<< trk_pt<<"+/-"<<trk_qOverperr << ", eta "<<trk_eta<<"+/-"<<trk_thetaerr<<", phi"<<trk_phi<<"+/-"<<trk_phierr<<", d0"<<trk_d0<<"+/-"<<trk_d0err<<", z0"<<trk_z0<<"+/-"<<trk_z0err <<endmsg;
+ mytrk.push_back(MyTrack(iTrack,trk_pt,trk_theta,trk_phi,trk_d0,trk_z0,trk_pterr,trk_thetaerr,trk_phierr,trk_d0err,trk_z0err));
+ }
+ return mytrk;
+}
+
+double FTK_DuplicateTrackRemovalTool::ctheta2eta(double cot){
+ double eta=0.0;
+ double tempX=0.0;
+
+ if(atan(1.0/cot) >= 0) tempX = atan(1.0/cot);
+ if(atan(1.0/cot) < 0) tempX = atan(1.0/cot)+M_PI;
+
+ eta = -log(tan(tempX/2.0));
+
+ return eta;
+}
+
diff --git a/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx b/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
index f376046dab2dd39d49f034fe50f97b6bb725c999..ee4c8bcdf7c0fc3635a6597d65705c74b5435457 100644
--- a/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
+++ b/Trigger/TrigFTK/FTK_RecTools/src/components/FTK_RecTools_entries.cxx
@@ -3,13 +3,16 @@
#include "FTK_RecTools/FTK_VertexFinderTool.h"
#include "FTK_RecTools/FTK_PixelClusterOnTrackTool.h"
#include "FTK_RecTools/FTK_SCTClusterOnTrackTool.h"
+#include "FTK_RecTools/FTK_DuplicateTrackRemovalTool.h"
DECLARE_TOOL_FACTORY(FTK_VertexFinderTool)
DECLARE_TOOL_FACTORY(FTK_PixelClusterOnTrackTool)
DECLARE_TOOL_FACTORY(FTK_SCTClusterOnTrackTool)
+DECLARE_TOOL_FACTORY(FTK_DuplicateTrackRemovalTool)
DECLARE_FACTORY_ENTRIES(FTK_RecTools) {
DECLARE_TOOL(FTK_VertexFinderTool)
DECLARE_TOOL(FTK_PixelClusterOnTrackTool)
DECLARE_TOOL(FTK_SCTClusterOnTrackTool)
+ DECLARE_TOOL(FTK_DuplicateTrackRemovalTool)
}