From 2840188565fd8b8166d5ce78ca5c1d7b82a36497 Mon Sep 17 00:00:00 2001
From: Andy Salnikov <salnikov@slac.stanford.edu>
Date: Sun, 22 Nov 2020 03:57:39 +0100
Subject: [PATCH] Improve vertex finding / clustering in beamspot calibration.
Clusterer now has an option to do clustering at beampot position
instead of (0,0) position. This improves clustering and makes resulting
vertex fit less biased, and matches offline beam position.
---
.../TrigT2BeamSpot/CMakeLists.txt | 5 +-
.../TrigT2BeamSpot/T2VertexBeamSpotTool.h | 8 +--
.../python/T2VertexBeamSpotConfig.py | 1 +
.../python/T2VertexBeamSpotMonitoring.py | 2 +
.../TrigT2BeamSpot/src/T2TrackClusterer.cxx | 63 ++++++++++++++++---
.../TrigT2BeamSpot/src/T2TrackClusterer.h | 26 +++++++-
.../src/T2VertexBeamSpotTool.cxx | 34 ++++++----
7 files changed, 112 insertions(+), 27 deletions(-)
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/CMakeLists.txt b/Trigger/TrigAlgorithms/TrigT2BeamSpot/CMakeLists.txt
index ce22c00a73fe..dcaab83db172 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/CMakeLists.txt
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/CMakeLists.txt
@@ -11,7 +11,10 @@ atlas_add_component( TrigT2BeamSpot
src/*.cxx
src/components/TrigT2*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
- LINK_LIBRARIES ${ROOT_LIBRARIES} AthContainers AthenaBaseComps AthenaMonitoringKernelLib BeamSpotConditionsData EventPrimitives GaudiKernel StoreGateLib TrigInDetEvent TrigInDetToolInterfacesLib TrigInterfacesLib TrigNavigationLib TrigSteeringEvent TrigTimeAlgsLib TrkParameters TrkTrack TrkTrackSummary xAODEventInfo )
+ LINK_LIBRARIES ${ROOT_LIBRARIES} AthContainers AthenaBaseComps AthenaMonitoringKernelLib
+ BeamSpotConditionsData EventPrimitives GaudiKernel StoreGateLib TrigInDetEvent
+ TrigInDetToolInterfacesLib TrigInterfacesLib TrigNavigationLib TrigSteeringEvent
+ TrigTimeAlgsLib TrkParameters TrkTrack TrkTrackSummary xAODEventInfo )
# Install files from the package:
atlas_install_python_modules( python/*.py POST_BUILD_CMD ${ATLAS_FLAKE8} )
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/TrigT2BeamSpot/T2VertexBeamSpotTool.h b/Trigger/TrigAlgorithms/TrigT2BeamSpot/TrigT2BeamSpot/T2VertexBeamSpotTool.h
index ba80ec2c7b4a..bbff59fcf926 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/TrigT2BeamSpot/T2VertexBeamSpotTool.h
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/TrigT2BeamSpot/T2VertexBeamSpotTool.h
@@ -29,6 +29,7 @@
#include "../src/T2Track.h"
#include "../src/T2BeamSpot.h"
#include "../src/T2SplitVertex.h"
+#include "../src/T2TrackClusterer.h"
//Athena tools
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
@@ -54,7 +55,6 @@ namespace HLT {
}
namespace PESA {
class T2SplitVertex;
- class T2TrackClusterer;
/**
* This class uses primary vertex reconstruction to measure
* and monitor the LHC beam as seen by the ATLAS detector.
@@ -116,6 +116,8 @@ namespace PESA {
bool m_splitWholeCluster;
bool m_reclusterSplit;
double m_trackSeedPt;
+ std::string m_clusterPerigee = "original"; // one of "original", "beamspot", "beamline"
+ T2TrackClusterer::TrackPerigee m_clusterTrackPerigee = T2TrackClusterer::perigee_original;
/* Track selection criteria */
unsigned m_totalNTrkMin;
@@ -177,10 +179,8 @@ namespace PESA {
std::string m_vertexCollName;
private:
- ToolHandle<GenericMonitoringTool> m_monTool{this,"MonTool","","Monitoring tool"};
-
-
+ ToolHandle<GenericMonitoringTool> m_monTool{this,"MonTool","","Monitoring tool"};
};
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotConfig.py b/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotConfig.py
index 9a82a28453f4..7b7a431bcbcf 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotConfig.py
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotConfig.py
@@ -21,6 +21,7 @@ InDetTrigMTBeamSpotTool = PESA__T2VertexBeamSpotTool(
MonTool = bsToolMonitoring,
WeightClusterZ = True, # Use the track Z0 weighted cluster Z position as seed
ReclusterSplit = False, # Recluster split track collections before vertex fitting
+ ClusterPerigee = "beamspot",
nSplitVertices = 2, # Turn on (>1) or off vertex splitting
TotalNTrackMin = 4, # Minimum number of tracks required in an event
TrackMinPt = 0.5, # Minimum track pT to be considered for vertexing
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotMonitoring.py b/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotMonitoring.py
index 94b7896b2184..87c857417a5a 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotMonitoring.py
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/python/T2VertexBeamSpotMonitoring.py
@@ -237,6 +237,8 @@ class T2VertexBeamSpotToolMonitoring(BaseMonitoringTool):
title="ClusterDeltaZ0; Delta of track Z0 and cluster Z position [mm]; Number of tracks")
self.makeHisto1D('ClusterZ0Pull', 'TH1I', 100, -10.0, 10.0,
title="ClusterZ0Pull; Pull of track Z0 with respect to cluster Z position; Number of tracks")
+ self.makeHisto1D('ClusterClusterDeltaZ0', 'TH1I', 200, -100.0, 100.0,
+ title="ClusterClusterDeltaZ0; Delta of cluster-cluster Z position [mm]; Number of clusters")
def defineVertexHistos(self, detail):
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.cxx b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.cxx
index 4454071b7e0c..fdce5e5daa8c 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.cxx
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.cxx
@@ -10,14 +10,34 @@
#include "TrkParameters/TrackParameters.h"
#include "EventPrimitives/EventPrimitivesHelpers.h"
+#include "BeamSpotConditionsData/BeamSpotData.h"
#include <cmath>
+#include <stdexcept>
+
using std::abs;
using std::sqrt;
+using std::sin;
+using std::cos;
+using std::tan;
using namespace PESA;
+T2TrackClusterer::TrackPerigee
+T2TrackClusterer::trackPerigeeFromString(const std::string& perigeeStr)
+{
+ if (perigeeStr == "original") {
+ return perigee_original;
+ } else if (perigeeStr == "beamspot") {
+ return perigee_beamspot;
+ } else if (perigeeStr == "beamline") {
+ return perigee_beamline;
+ } else {
+ throw std::runtime_error("Invalid value of ClusterPerigee parameter: " + perigeeStr);
+ }
+}
+
double
T2TrackClusterer::trackWeight( const Trk::Track& track ) const
{
@@ -27,7 +47,7 @@ T2TrackClusterer::trackWeight( const Trk::Track& track ) const
}
const TrackCollection&
-T2TrackClusterer::cluster( const TrackCollection& tracks )
+T2TrackClusterer::cluster( const TrackCollection& tracks, const InDet::BeamSpotData* beamspot )
{
m_seedZ0 = 0.;
m_totalZ0Err = 0.;
@@ -41,9 +61,9 @@ T2TrackClusterer::cluster( const TrackCollection& tracks )
}
const Trk::Track* seedTrack = *tracks.begin();
-
- const Trk::TrackParameters* seedTrackPars = seedTrack->perigeeParameters();
- const double seedPT = std::abs(sin(seedTrackPars->parameters()[Trk::theta])/seedTrackPars->parameters()[Trk::qOverP]);
+
+ auto& params = seedTrack->perigeeParameters()->parameters();
+ const double seedPT = std::abs(sin(params[Trk::theta])/params[Trk::qOverP]);
if ( seedPT < m_minPT )
{
@@ -51,16 +71,16 @@ T2TrackClusterer::cluster( const TrackCollection& tracks )
return *m_cluster.asDataVector();
}
+
double sumWeight = trackWeight( *seedTrack );
- m_seedZ0 = seedTrackPars->parameters()[Trk::z0]; // becomes the weighted-average z0 of the cluster
+ m_seedZ0 = trackPerigeeZ0(*seedTrack, beamspot);
m_cluster.push_back( seedTrack );
for ( TrackCollection::const_iterator track = tracks.begin() + 1;
track != tracks.end(); ++track )
{
- const Trk::TrackParameters* trackPars = (*track)->perigeeParameters();
- const double trackZ0 = trackPars->parameters()[Trk::z0];
+ const double trackZ0 = trackPerigeeZ0(**track, beamspot);
const double deltaZ = trackZ0 - m_seedZ0;
const double weight = trackWeight( **track );
@@ -81,3 +101,32 @@ T2TrackClusterer::cluster( const TrackCollection& tracks )
return *m_cluster.asDataVector();
}
+
+double
+T2TrackClusterer::trackPerigeeZ0(const Trk::Track& track, const InDet::BeamSpotData* beamspot) const
+{
+ auto& params0 = track.perigeeParameters()->parameters();
+ if (m_trackPerigee == perigee_original or beamspot == nullptr) {
+ return params0[Trk::z0];
+ } else if (m_trackPerigee == perigee_beamspot or m_trackPerigee == perigee_beamline) {
+ // TODO: beamline is doing the same as beamspot for now, I think it should be OK
+ // with our track resolution, somethig to check in the future (when I retire)
+
+ // This assumes that track perigee is defined at 0,0, transform is:
+ // z = z_0 + (B_x*cos(phi) + B_y*sin(phi)) / tan(theta)
+
+ double z_0 = params0[Trk::z0];
+ double phi = params0[Trk::phi];
+ double theta = params0[Trk::theta];
+ auto& beamPos = beamspot->beamPos();
+ double B_x = beamPos[0];
+ double B_y = beamPos[1];
+
+ double z = z_0 + (B_x*cos(phi) + B_y*sin(phi)) / tan(theta);
+ return z;
+
+ } else {
+ // fallback to "original"
+ return params0[Trk::z0];
+ }
+}
\ No newline at end of file
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.h b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.h
index 130345b740f8..d58500686d0e 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.h
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2TrackClusterer.h
@@ -33,6 +33,8 @@
#include <vector>
#include <cmath>
+namespace InDet { class BeamSpotData; }
+
namespace PESA
{
@@ -40,14 +42,28 @@ namespace PESA
{
public:
+ // option for calculating track perigee parameters
+ enum TrackPerigee {
+ perigee_original, // Use original track perigee parameters, usually means
+ // perigee w.r.t. coordinate system origin x,y=0,0
+ perigee_beamspot, // Perigee to a beamspot position, ignoring tilt
+ perigee_beamline, // Perigee to a beam line, including tilt
+ };
+
+ // convert string to above enum value
+ static TrackPerigee trackPerigeeFromString(const std::string& perigeeStr);
+
// Constructor
- T2TrackClusterer( double deltaZ = 10.*Gaudi::Units::mm, double minPT = 1.*Gaudi::Units::GeV, bool weightedZ = true, unsigned maxSize = 10000. )
+ // clusterPerigee is one of '0'=origin, 's'=beamspot, 'l'=beamline
+ T2TrackClusterer( double deltaZ = 10.*Gaudi::Units::mm, double minPT = 1.*Gaudi::Units::GeV, bool weightedZ = true, unsigned maxSize = 10000.,
+ TrackPerigee trackPerigee = perigee_original)
: m_deltaZ ( deltaZ )
, m_minPT ( minPT )
, m_weightedZ ( weightedZ )
, m_maxSize ( maxSize )
, m_seedZ0 ( 0. )
, m_totalZ0Err( 0. )
+ , m_trackPerigee( trackPerigee )
{}
// Accessors
@@ -60,9 +76,13 @@ namespace PESA
// Methods
double trackWeight( const Trk::Track& track ) const;
- const TrackCollection& cluster( const TrackCollection& tracks );
+ const TrackCollection& cluster( const TrackCollection& tracks, const InDet::BeamSpotData* beamspot = nullptr );
private:
+
+ // return perigee z0 for a track
+ double trackPerigeeZ0(const Trk::Track& track, const InDet::BeamSpotData* beamspot) const;
+
// Data members
const double m_deltaZ;
const double m_minPT;
@@ -72,6 +92,8 @@ namespace PESA
double m_seedZ0;
double m_totalZ0Err;
+ const TrackPerigee m_trackPerigee = perigee_original;
+
ConstDataVector<TrackCollection> m_cluster;
ConstDataVector<TrackCollection> m_unusedTracks;
};
diff --git a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2VertexBeamSpotTool.cxx b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2VertexBeamSpotTool.cxx
index e9b9fa304a8c..c6a813924cc6 100644
--- a/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2VertexBeamSpotTool.cxx
+++ b/Trigger/TrigAlgorithms/TrigT2BeamSpot/src/T2VertexBeamSpotTool.cxx
@@ -14,7 +14,6 @@
//============================================================
// This algorithm
#include "TrigT2BeamSpot/T2VertexBeamSpotTool.h"
-#include "T2TrackClusterer.h"
#include "T2TrackManager.h"
#include "T2Timer.h"
// Specific to this algorithm
@@ -85,6 +84,7 @@ PESA::T2VertexBeamSpotTool::T2VertexBeamSpotTool( const std::string& type, const
declareProperty("WeightClusterZ", m_weightSeed = true );
declareProperty("ReclusterSplit", m_reclusterSplit = true );
declareProperty("SplitWholeCluster", m_splitWholeCluster = false );
+ declareProperty("ClusterPerigee", m_clusterPerigee = "beamspot" );
// Track Selection criteria
declareProperty("TrackMinPt", m_minTrackPt = 0.7*GeV );
@@ -125,9 +125,6 @@ PESA::T2VertexBeamSpotTool::T2VertexBeamSpotTool( const std::string& type, const
// Single interaction trigger
declareProperty("minNpvTrigger", m_minNpvTrigger = 0);
declareProperty("maxNpvTrigger", m_maxNpvTrigger = 2);
-
-
-
}
@@ -145,6 +142,8 @@ StatusCode T2VertexBeamSpotTool::initialize(){
//Retrieve primary vertex fitter tool
ATH_CHECK( m_primaryVertexFitterTool.retrieve() );
+ m_clusterTrackPerigee = T2TrackClusterer::trackPerigeeFromString(m_clusterPerigee);
+
return StatusCode::SUCCESS;
}
@@ -253,9 +252,18 @@ unsigned int T2VertexBeamSpotTool::reconstructVertices( ConstDataVector<TrackCol
auto mon = Monitored::Group(m_monTool, timeToSortTracks );
}
+ // Extract beam spot parameters
+ SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
+ const InDet::BeamSpotData* indetBeamSpot(*beamSpotHandle);
+ const T2BeamSpot beamSpot(indetBeamSpot);
+ ATH_MSG_DEBUG( "Beamspot from BeamCondSvc: " << beamSpot);
+
// Prepare a track clustering algorithm with the given parameters
// Should we leave this as a standalone class or merge with the tool?
- T2TrackClusterer trackClusterer( m_trackClusDZ, m_trackSeedPt, m_weightSeed, m_vtxNTrkMax );
+ T2TrackClusterer trackClusterer( m_trackClusDZ, m_trackSeedPt, m_weightSeed, m_vtxNTrkMax,
+ m_clusterTrackPerigee );
+
+ std::vector<double> clusterZ0; // Z0 for each cluster, for monitoring only
// Create clusters from the track collection until all its tracks are used
while ( ! mySelectedTrackCollection.empty() ) {
@@ -270,7 +278,7 @@ unsigned int T2VertexBeamSpotTool::reconstructVertices( ConstDataVector<TrackCol
// as the seed.
{
auto timeToZCluster = Monitored::Timer("TIME_toZCluster");
- trackClusterer.cluster( *mySelectedTrackCollection.asDataVector() );
+ trackClusterer.cluster( *mySelectedTrackCollection.asDataVector(), indetBeamSpot );
auto mon = Monitored::Group(m_monTool, timeToZCluster );
}
@@ -330,13 +338,6 @@ unsigned int T2VertexBeamSpotTool::reconstructVertices( ConstDataVector<TrackCol
// Update vertex counter
nVtx++;
- // Extract beam spot parameters
- SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
-
- T2BeamSpot beamSpot(*beamSpotHandle);
-
- ATH_MSG_DEBUG( "Beamspot from BeamCondSvc: " << beamSpot);
-
const T2Vertex myVertex( *primaryVertex, vertexTracks, beamSpot, trackClusterer.seedZ0() );
// Monitor all vertices parameters
@@ -364,6 +365,13 @@ unsigned int T2VertexBeamSpotTool::reconstructVertices( ConstDataVector<TrackCol
auto deltaVtxZ = Monitored::Scalar<double> ( "ClusterDeltaVertexZ", trackClusterer.seedZ0() - myVertex.Z() );
auto mon = Monitored::Group(m_monTool, deltaVtxZ );
+ // monitor cluster-cluster delta Z0
+ for (double prevClusterZ0: clusterZ0) {
+ auto clusterClusterDeltaZ = Monitored::Scalar<double>("ClusterClusterDeltaZ0", trackClusterer.seedZ0() - prevClusterZ0);
+ auto mon = Monitored::Group(m_monTool, clusterClusterDeltaZ);
+ }
+ clusterZ0.push_back(trackClusterer.seedZ0());
+
// If the vertex is good, splits are requested, and we have enough tracks to split, split them!
if ( passVertex && m_nSplitVertices > 1 ) {
--
GitLab