Merge branch 'tickets/ATONLBS-30-master' into 'master'

Improve vertex finding / clustering in beamspot calibration.

Closes ATONLBS-30

See merge request atlas/athena!38957

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} )

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"};
 
    };
 

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

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):
 

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

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;
   };

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 ) {