Commit e61e4e77 authored by Tadej Novak's avatar Tadej Novak
Browse files

Merge branch 'hgtd_extension_tool_21p9' into '21.9'

Track extension tool

See merge request !46020
parents d9ef2c1a 34a19c3e
/**
* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
*
* @file HGTD_RecToolInterfaces/IHGTD_TrackTimeExtensionTool.h
* @author Alexander Leopold <alexander.leopold@cern.ch>
* @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!!!
*/
#ifndef IHGTD_TRACKTIMEEXTENSIONTOOL_H
#define IHGTD_TRACKTIMEEXTENSIONTOOL_H
#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 {
public:
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
#endif // IHGTD_TRACKTIMEEXTENSIONTOOL_H
################################################################################
# Package: HGTD_TrackTimeExtensionTools
################################################################################
# Declare the package name:
atlas_subdir( HGTD_TrackTimeExtensionTools )
# Component(s) in the package:
atlas_add_component( HGTD_TrackTimeExtensionTools
src/*.cxx
src/components/*.cxx
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 <noemi.calace@cern.ch>
* @author Alexander Leopold <alexander.leopold@cern.ch>
*
* @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.
*
* TODO:
* - 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!
*/
#ifndef HGTD_ITERATIVEEXTENTSIONTOOL_H
#define HGTD_ITERATIVEEXTENTSIONTOOL_H
#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 {
public:
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;
private:
/**
* @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>
updateStateWithBestFittingCluster(
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
#endif // HGTD_ITERATIVEEXTENTSIONTOOL_H
/**
* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
*
* @file HGTD_TrackTimeExtensionTools/src/HGTD_IterativeExtensionTool.cxx
* @author Noemi Calace <noemi.calace@cern.ch>
* @author Alexander Leopold <alexander.leopold@cern.ch>
* @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>
HGTD::HGTD_IterativeExtensionTool::HGTD_IterativeExtensionTool(
const std::string& t, const std::string& n, const IInterface* p)
: AthAlgTool(t, n, p),
m_extrapolator("Trk::Extrapolator/AtlasExtrapolator"),
m_updator("Trk::KalmanUpdator/KalmanUpdator"),
m_tof_corr_tool(
"StraightLineTOFcorrectionTool/StraightLineTOFcorrectionTool"),
m_track(nullptr),
m_cluster_container(nullptr),
m_hgtd_det_mgr(nullptr),
m_hgtd_id_helper(nullptr),
m_chi2_cut(5.0),
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();
ATH_CHECK(m_extrapolator.retrieve());
ATH_CHECK(m_updator.retrieve());
// 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>
HGTD::HGTD_IterativeExtensionTool::extendTrackToHGTD(
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 = layers_vec.at(position);
// only use it, if it is an active one
if (layer->layerType() != 1) {
continue;
}
hgtd_layer_i++;
const Trk::Surface& surf_obj = layer->surfaceRepresentation();
ATH_MSG_DEBUG("extrapolation surface z: " << surf_obj.center().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;
extrap_result.reset(m_extrapolator->extrapolate(
*last_param, surf_obj, Trk::PropDirection::alongMomentum, false,
track.info().particleHypothesis()));
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 =
updateStateWithBestFittingCluster(extrapolated_params);
if (not updated_state) {
continue;
}
// 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
extension.at(hgtd_layer_i) = std::move(updated_state);
}
return extension;
}
const Trk::TrackStateOnSurface*
HGTD::HGTD_IterativeExtensionTool::getLastHitOnTrack(const Trk::Track& track) {
const DataVector<const Trk::TrackStateOnSurface>* tsos =
track.trackStateOnSurfaces();
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) {
continue;
}
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*>
HGTD::HGTD_IterativeExtensionTool::getCompatibleSurfaces(
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);
surfaces.push_back(module_surface);
// pick up additional surfaces in a 4cm radius
// TODO REPLACE THIS BY NEIGHBOUR FUNCTIONALITY IN FUTURE!!
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()) {
surfaces.push_back(additional_surface);
}
}
return surfaces;
}
std::vector<std::unique_ptr<const Trk::TrackParameters>>
HGTD::HGTD_IterativeExtensionTool::extrapolateToSurfaces(
const Trk::TrackParameters& param,
const std::vector<const Trk::Surface*>& surfaces) {
std::vector<std::unique_ptr<const Trk::TrackParameters>> params;
params.reserve(surfaces.size());
for (const auto* surface : surfaces) {
std::unique_ptr<const Trk::TrackParameters> extrapolated_params = nullptr;
extrapolated_params.reset(m_extrapolator->extrapolate(
param, *surface, Trk::alongMomentum, false,
static_cast<Trk::ParticleHypothesis>(m_particle_hypot)));
params.push_back(std::move(extrapolated_params));
}
return params;
}
std::unique_ptr<const Trk::TrackStateOnSurface>
HGTD::HGTD_IterativeExtensionTool::updateStateWithBestFittingCluster(
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 =
findBestCompatibleCluster(param.get());
if (not best_tsos) {
ATH_MSG_DEBUG("[updateStateWithBestFittingCluster] tsos is null");
continue;
}
ATH_MSG_DEBUG("[updateStateWithBestFittingCluster] tsos found");
double chi2 = best_tsos->fitQualityOnSurface()->chiSquared() /
best_tsos->fitQualityOnSurface()->doubleNumberDoF();
ATH_MSG_DEBUG(
"[updateStateWithBestFittingCluster] found state with chi2 of "
<< chi2);
// make sure only one TSOS is kept
if (!updated_state or chi2 < lowest_chi2) {
updated_state.swap(best_tsos);
lowest_chi2 = chi2;
}
}
return std::move(updated_state);
}
std::unique_ptr<const Trk::TrackStateOnSurface>
HGTD::HGTD_IterativeExtensionTool::findBestCompatibleCluster(
const Trk::TrackParameters* param) {