Commit 34a19c3e authored by Alexander Leopold's avatar Alexander Leopold Committed by Tadej Novak
Browse files

Track extension tool

parent 7001d426
* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
* @file HGTD_RecToolInterfaces/IHGTD_TrackTimeExtensionTool.h
* @author Alexander Leopold <>
* @author
* @date August, 2021
* @brief Tracks reconstructed in ITk get extended to the surfaces of HGTD,
* where compatible hits get selected. In the extension process, the track
* is not altered. Currently the information from the selected hits in HGTD
* is retrieved and decorations added to the track. After this, the newly
* found track states are deleted.
* This might be changed in the future, where a dedicated object is written out,
* and a connection (link) between a track and its extension to HGTD stored.
* - FIXME: retrieveClusters should be removed from the interface and called
* once at the start of the event in the derived tool!!!
#include "GaudiKernel/IAlgTool.h"
#include "HGTD_PrepRawData/HGTD_ClusterContainer.h"
#include "TrkTrack/Track.h"
#include "TrkTrack/TrackStateOnSurface.h"
#include <memory>
class MsgStream;
namespace HGTD {
static const InterfaceID
IID_IHGTD_TrackExtensionTool("HGTD::IHGTD_TrackTimeExtensionTool", 1, 0);
class IHGTD_TrackTimeExtensionTool : virtual public IAlgTool {
static const InterfaceID& interfaceID();
* @brief Extends a track to (up to) 4 measurements in HGTD.
* @param [in] track Track built in the inner tracker to be extended to HGTD.
* @param [in] container Container of HGTD_Cluster objects.
* @return Returns 4 measurements in HGTD (can contain holes or outliers)
* associated to the track.
virtual std::array<std::unique_ptr<const Trk::TrackStateOnSurface>, 4>
extendTrackToHGTD(const Trk::Track& track,
const HGTD::HGTD_ClusterContainer* container) = 0;
/** Inline methods **/
inline const InterfaceID& IHGTD_TrackTimeExtensionTool::interfaceID() {
return IID_IHGTD_TrackExtensionTool;
} // namespace HGTD
# Package: HGTD_TrackTimeExtensionTools
# Declare the package name:
atlas_subdir( HGTD_TrackTimeExtensionTools )
# Component(s) in the package:
atlas_add_component( HGTD_TrackTimeExtensionTools
LINK_LIBRARIES AthenaBaseComps GaudiKernel InDetPrepRawData
TrkGeometry TrkEventPrimitives TrkSurfaces TrkTrack
HGTD_RecToolInterfaces TrkExInterfaces HGTD_PrepRawData
TrkDetDescrUtils HGTD_ReadoutGeometry HGTD_Identifier
HGTD_RIO_OnTrack TrkToolInterfaces)
# Install files from the package:
atlas_install_headers( HGTD_TrackTimeExtensionTools )
* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
* @file HGTD_TrackTimeExtensionTools/HGTD_IterativeExtensionTool.h
* @author Noemi Calace <>
* @author Alexander Leopold <>
* @date August, 2021
* @brief The HGTD_IterativeExtensionTool uses an iterative Kalman filter
* approach to find hits (HGTD_Cluster) in HGTD that are compatible with a track
* extrapolated to the surfaces of HGTD. During this process, the track itself
* is not altered, the found hits are returned as an array of
* HGTD_ClusterOnTrack objects.
* - exchange chi2/ndof with TMath::Prob?
* - add chi2 for the time -> get an estimate at layer 0 (TOF +- beamspot),
* decorate the TrackParameters with time and timeres before!
#include "AthenaBaseComps/AthAlgTool.h"
#include "HGTD_RecToolInterfaces/IHGTD_TrackTimeExtensionTool.h"
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
#include "HGTD_PrepRawData/HGTD_ClusterContainer.h"
#include "HGTD_RecToolInterfaces/IHGTD_TOFcorrectionTool.h"
#include "TrkExInterfaces/IExtrapolator.h"
#include "TrkToolInterfaces/IUpdator.h"
#include <memory>
class HGTD_ID;
class HGTD_DetectorManager;
namespace HGTD {
class HGTD_IterativeExtensionTool : virtual public IHGTD_TrackTimeExtensionTool,
public AthAlgTool {
HGTD_IterativeExtensionTool(const std::string&, const std::string&,
const IInterface*);
virtual ~HGTD_IterativeExtensionTool() = default;
virtual StatusCode initialize();
* @brief Finds the (up to) four measurements in HGTD that can be associated
* to the track. The Track itself is not altered during this process.
* Since HGTD contains four active layers per endcap, four hits is the
* maximum. An extension can be of hole-type
* @param [in] track Track reconstructed in ITk.
* @return Array of compatible HGTD hits in the form of HGTD_ClusterOnTrack.
virtual std::array<std::unique_ptr<const Trk::TrackStateOnSurface>, 4>
extendTrackToHGTD(const Trk::Track& track,
const HGTD::HGTD_ClusterContainer* container) override;
MsgStream& dump(MsgStream& out) const;
std::ostream& dump(std::ostream& out) const;
* @brief Retrieve the last hit on track stored in the Trk::Track.
* @param [in] track Track reconstructed in ITk.
* @return The last hit on track, i.e. the closest one to HGTD.
* FIXME: returns a const*, but the Track manages this TSOS. Could I return
* rather const&? But: I want to check if type is OK... should I return an
* "empty" TSOS if nothing passes the requirements?
const Trk::TrackStateOnSurface* getLastHitOnTrack(const Trk::Track& track);
* @brief Select all within the vincinity of the point of extrapolation.
* @param [in] param The already extrapolated track parameter.
* @param [in] layer Target layer of the extrapolation.
* @return Vector of compatible surfaces that represent modules.
* FIXME exchange the loop with the neighbour functionality (wherever that is)
std::vector<const Trk::Surface*>
getCompatibleSurfaces(const Trk::TrackParameters& param,
const Trk::Layer* layer);
* @brief After the compatible surfaces have been selected, the extrapolation
* is repeated to each of them.
* @param [in] param The last track parameter before extrapolation.
* @param [in] surfaces Target layer of the extrapolation.
* @return Vector of the updated track parameters.
* FIXME Is this really necessary to call the extrapolation on each target
* surface?
std::vector<std::unique_ptr<const Trk::TrackParameters>>
extrapolateToSurfaces(const Trk::TrackParameters& param,
const std::vector<const Trk::Surface*>& surfaces);
* @brief Finds the overall best fitting cluster by keeping the one that gave
* the lowest chi2 after testing each surface and returning it.
std::unique_ptr<const Trk::TrackStateOnSurface>
const std::vector<std::unique_ptr<const Trk::TrackParameters>>& params);
* @brief Find the cluster on a given surface that has the best chi2 passing
* the cut. Is called within updateStateWithBestFittingCluster.
std::unique_ptr<const Trk::TrackStateOnSurface>
findBestCompatibleCluster(const Trk::TrackParameters* param);
* @brief Calls the TOF correction tool to build an HGTD_ClusterOnTrack with a
* calibrated time and resolution and then calls the Kalman updator tools to
* update the track state.
std::unique_ptr<const Trk::TrackStateOnSurface>
updateState(const Trk::TrackParameters* param,
const HGTD::HGTD_Cluster* cluster);
// extrapolation tool
ToolHandle<Trk::IExtrapolator> m_extrapolator;
// kalman updator tool
ToolHandle<Trk::IUpdator> m_updator;
ToolHandle<HGTD::IHGTD_TOFcorrectionTool> m_tof_corr_tool;
// keep a pointer to the currently used track, but does not own it!
// FIXME: this is needed for the TOF correction. Maybe there is a smarter way
// without keeping this pointer and without pushing the Trk::Track object
// through everything ==>> move to a TOF corr. accepting just the perigee??
const Trk::Track* m_track;
* @brief Keep pointer to current cluster container.
const HGTD::HGTD_ClusterContainer* m_cluster_container;
const HGTD_DetectorManager* m_hgtd_det_mgr;
const HGTD_ID* m_hgtd_id_helper;
* @brief Track extensions are only kept if the chi2/ndof is lower than the
* defined cut.
float m_chi2_cut;
* @brief Particle hypothesis used for the extrapolation.
int m_particle_hypot;
} // namespace HGTD
* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
* @file HGTD_TrackTimeExtensionTools/src/HGTD_IterativeExtensionTool.cxx
* @author Noemi Calace <>
* @author Alexander Leopold <>
* @date August, 2021
* @brief
#include "HGTD_TrackTimeExtensionTools/HGTD_IterativeExtensionTool.h"
#include "HGTD_Identifier/HGTD_ID.h"
#include "HGTD_RIO_OnTrack/HGTD_ClusterOnTrack.h"
#include "HGTD_ReadoutGeometry/HGTD_DetectorManager.h"
#include "TrkDetDescrUtils/GeometrySignature.h"
#include "TrkGeometry/DiscLayer.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkTrack/TrackStateOnSurface.h"
#include <iostream>
#include <vector>
const std::string& t, const std::string& n, const IInterface* p)
: AthAlgTool(t, n, p),
m_particle_hypot(Trk::ParticleHypothesis::pion) {
declareProperty("ExtrapolatorTool", m_extrapolator);
declareProperty("UpdatorTool", m_updator);
declareProperty("TOFCorrTool", m_tof_corr_tool);
declareProperty("Chi2Cut", m_chi2_cut);
declareProperty("ParticleHypothesis", m_particle_hypot);
StatusCode HGTD::HGTD_IterativeExtensionTool::initialize() {
StatusCode sc = AthAlgTool::initialize();
// get HGTD Detector Description Manager and HGTD Helper
ATH_CHECK(detStore()->retrieve(m_hgtd_det_mgr, "HGTD"));
ATH_CHECK(detStore()->retrieve(m_hgtd_id_helper, "HGTD_ID"));
return sc;
std::array<std::unique_ptr<const Trk::TrackStateOnSurface>, 4>
const Trk::Track& track, const HGTD::HGTD_ClusterContainer* container) {
ATH_MSG_DEBUG("Start extending");
m_track = &track;
m_cluster_container = container;
std::array<std::unique_ptr<const Trk::TrackStateOnSurface>, 4> extension{
nullptr, nullptr, nullptr, nullptr};
// get last hit on track in ITk as a starting point for the extrapolation
const Trk::TrackStateOnSurface* last_hit = getLastHitOnTrack(track);
const Trk::TrackParameters* last_param = last_hit->trackParameters();
ATH_MSG_DEBUG("last param x: " << last_param->position().x()
<< ", y: " << last_param->position().y()
<< ", z: " << last_param->position().z());
ATH_MSG_DEBUG("last param px: " << last_param->momentum().x()
<< ", py: " << last_param->momentum().y()
<< ", pz: " << last_param->momentum().z());
// get the tracking geometry
const Trk::TrackingGeometry* trk_geom = m_extrapolator->trackingGeometry();
if (not trk_geom) {
ATH_MSG_DEBUG("trackingGeometry returns null");
return extension;
bool is_pos_endcap = last_param->eta() > 0;
// get the target volume
const Trk::TrackingVolume* hgtd_trk_volume = trk_geom->trackingVolume(
is_pos_endcap ? "HGTD::PositiveEndcap" : "HGTD::NegativeEndcap");
if (not hgtd_trk_volume) {
ATH_MSG_DEBUG("trackingVolume returns null");
return extension;
const auto* confined_layers =
hgtd_trk_volume->confinedLayers(); // BinnedArray<Trk::Layer>
// careful, this array is not ordered from inside out (only in pos endcap)
if (not confined_layers) {
ATH_MSG_DEBUG("confinedLayers returns null");
return extension;
// get the layers, traverse depending on endcap used
// since they are not in ascending z order!!
const auto layers_vec =
confined_layers->arrayObjects(); // is std::vector<Layer*>
size_t layer_vec_size = layers_vec.size();
short position = is_pos_endcap ? 0 : layer_vec_size - 1;
short step = is_pos_endcap ? 1 : -1;
short hgtd_layer_i = -1; // start at -1 since the ++ happens at the beginning
for (size_t i = 0; i < layer_vec_size - 1; i++, position = position + step) {
const Trk::Layer* layer =;
// only use it, if it is an active one
if (layer->layerType() != 1) {
const Trk::Surface& surf_obj = layer->surfaceRepresentation();
ATH_MSG_DEBUG("extrapolation surface z: " <<;
// if it is a good surface, extrapolate to this surface
// uses same particle hypothesis used for the track itself
// TODO: BoundaryCheck set to false as in 20.20 -> what does this do?
std::unique_ptr<const Trk::TrackParameters> extrap_result = nullptr;
*last_param, surf_obj, Trk::PropDirection::alongMomentum, false,;
ATH_MSG_DEBUG("extrap. params x: "
<< extrap_result->position().x()
<< ", y: " << extrap_result->position().y()
<< ", z: " << extrap_result->position().z());
ATH_MSG_DEBUG("extrap. params px: "
<< extrap_result->momentum().x()
<< ", py: " << extrap_result->momentum().y()
<< ", pz: " << extrap_result->momentum().z());
// find active surfaces in the vincinity
// const auto* surf_arr = layer->surfaceArray(); // these are the modules
// from this binned object, get the module closeby
auto compatible_surfaces = getCompatibleSurfaces(*extrap_result, layer);
auto extrapolated_params =
extrapolateToSurfaces(*last_param, compatible_surfaces);
std::unique_ptr<const Trk::TrackStateOnSurface> updated_state =
if (not updated_state) {
// if the state was updated with a measurement, the it becomes the new last
// parameter
last_param = updated_state->trackParameters();
// store the last track state to be returned = std::move(updated_state);
return extension;
const Trk::TrackStateOnSurface*
HGTD::HGTD_IterativeExtensionTool::getLastHitOnTrack(const Trk::Track& track) {
const DataVector<const Trk::TrackStateOnSurface>* tsos =
if (not tsos) {
ATH_MSG_ERROR("Failed to retrieve track state on surfaces");
return nullptr;
// loop over the associated hits in ITk in reverse order, since we want to
// select the one closest to HGTD to start the extrapolation
for (auto i = tsos->rbegin(); i != tsos->rend(); ++i) {
auto curr_last_tsos = *i;
if (not curr_last_tsos) {
if (curr_last_tsos->type(Trk::TrackStateOnSurface::Measurement) and
curr_last_tsos->trackParameters() and
curr_last_tsos->measurementOnTrack()) {
return curr_last_tsos;
return nullptr;
std::vector<const Trk::Surface*>
const Trk::TrackParameters& param, const Trk::Layer* layer) {
std::vector<const Trk::Surface*> surfaces;
// point of extrapolation as global position
Amg::Vector3D position = param.position();
// get the surface at the point of extrapolation
const auto* surface_arr = layer->surfaceArray();
const Trk::Surface* module_surface = surface_arr->object(position);
// pick up additional surfaces in a 4cm radius
short steps = 16;
float delta_angle = 6.2831853 / (float)steps;
float angle = 0.;
float radius = 20.0;
for (short i = 0; i < steps; i++, angle = angle + delta_angle) {
Amg::Vector3D delta(radius * std::cos(angle), radius * std::sin(angle), 0);
Amg::Vector3D result = position + delta;
const Trk::Surface* additional_surface = surface_arr->object(result);
// check if surface was added to avoid duplicates
if (std::find(surfaces.begin(), surfaces.end(), additional_surface) ==
surfaces.end()) {
return surfaces;
std::vector<std::unique_ptr<const Trk::TrackParameters>>
const Trk::TrackParameters& param,
const std::vector<const Trk::Surface*>& surfaces) {
std::vector<std::unique_ptr<const Trk::TrackParameters>> params;
for (const auto* surface : surfaces) {
std::unique_ptr<const Trk::TrackParameters> extrapolated_params = nullptr;
param, *surface, Trk::alongMomentum, false,
return params;
std::unique_ptr<const Trk::TrackStateOnSurface>
const std::vector<std::unique_ptr<const Trk::TrackParameters>>& params) {
ATH_MSG_DEBUG("[updateStateWithBestFittingCluster] start updating");
std::unique_ptr<const Trk::TrackStateOnSurface> updated_state = nullptr;
double lowest_chi2 = -1.;
// all compatible surfaces are tested for the best fitting cluster
for (auto&& param : params) {
std::unique_ptr<const Trk::TrackStateOnSurface> best_tsos =
if (not best_tsos) {
ATH_MSG_DEBUG("[updateStateWithBestFittingCluster] tsos is null");
ATH_MSG_DEBUG("[updateStateWithBestFittingCluster] tsos found");
double chi2 = best_tsos->fitQualityOnSurface()->chiSquared() /
"[updateStateWithBestFittingCluster] found state with chi2 of "
<< chi2);
// make sure only one TSOS is kept
if (!updated_state or chi2 < lowest_chi2) {
lowest_chi2 = chi2;
return std::move(updated_state);
std::unique_ptr<const Trk::TrackStateOnSurface>
const Trk::TrackParameters* param) {