diff --git a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt index e6f115dcf144c2023bdd2a18c6ef323157f052df..66b2d2037c1bb60e3051d12c334c2275ae47e2b5 100755 --- a/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt +++ b/Tracking/Acts/FaserActsKalmanFilter/CMakeLists.txt @@ -30,6 +30,7 @@ atlas_add_component(FaserActsKalmanFilter FaserActsKalmanFilter/ITrackFinderTool.h FaserActsKalmanFilter/TruthTrackFinderTool.h FaserActsKalmanFilter/Measurement.h + FaserActsKalmanFilter/SegmentFitTrackFinderTool.h FaserActsKalmanFilter/SimWriterTool.h FaserActsKalmanFilter/SPSeedBasedInitialParameterTool.h FaserActsKalmanFilter/SPSimpleInitialParameterTool.h @@ -39,6 +40,7 @@ atlas_add_component(FaserActsKalmanFilter FaserActsKalmanFilter/TruthSeededTrackFinderTool.h src/CombinatorialKalmbanFilterAlg.cxx src/FaserActsKalmanFilterAlg.cxx + src/SegmentFitTrackFinderTool.cxx src/SimWriterTool.cxx src/SPSeedBasedInitialParameterTool.cxx src/SPSimpleInitialParameterTool.cxx diff --git a/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SegmentFitTrackFinderTool.h b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SegmentFitTrackFinderTool.h new file mode 100644 index 0000000000000000000000000000000000000000..ff874e9374a7abd5561f74e1cb89a4eb5094894d --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/FaserActsKalmanFilter/SegmentFitTrackFinderTool.h @@ -0,0 +1,84 @@ +#ifndef FASERACTSKALMANFILTER_SEGMENTFITTRACKFINDERTOOL_H +#define FASERACTSKALMANFILTER_SEGMENTFITTRACKFINDERTOOL_H + +#include "AthenaBaseComps/AthAlgTool.h" +#include "Gaudi/Property.h" +#include "GaudiKernel/EventContext.h" +#include "GaudiKernel/IInterface.h" +#include "GaudiKernel/StatusCode.h" +#include "StoreGate/ReadHandleKey.h" +#include <string> +#include <vector> + +#include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h" +#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h" +#include "FaserActsKalmanFilter/ITrackFinderTool.h" +#include "TrkTrack/TrackCollection.h" + +class FaserSCT_ID; +namespace TrackerDD { +class SCT_DetectorManager; +} + + +class SegmentFitTrackFinderTool : public extends<AthAlgTool, ITrackFinderTool> { +public: + SegmentFitTrackFinderTool(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~SegmentFitTrackFinderTool() = default; + virtual StatusCode initialize() override; + virtual StatusCode finalize() override; + virtual StatusCode run() override; + + virtual const std::shared_ptr<const Acts::CurvilinearTrackParameters> initialTrackParameters() const override; + virtual const std::shared_ptr<const Acts::Surface> initialSurface() const override; + virtual const std::shared_ptr<std::vector<IndexSourceLink>> sourceLinks() const override; + virtual const std::shared_ptr<std::vector<Measurement>> measurements() const override; + +private: + std::shared_ptr<const Acts::CurvilinearTrackParameters> m_initialTrackParameters; + std::shared_ptr<const Acts::Surface> m_initialSurface; + std::shared_ptr<std::vector<IndexSourceLink>> m_sourceLinks {}; + std::shared_ptr<std::vector<Measurement>> m_measurements {}; + + const FaserSCT_ID* m_idHelper {nullptr}; + const TrackerDD::SCT_DetectorManager* m_detManager {nullptr}; + static double momentum(const std::map<int, Amg::Vector3D>& pos, double B=0.57) ; + + ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool { + this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"}; + SG::ReadHandleKey<TrackCollection> m_trackCollection { + this, "TrackCollection", "SegmentFit", "Input track collection name" }; + SG::ReadHandleKey<FaserSCT_SpacePointContainer> m_spacePointContainerKey { + this, "SpacePoints", "SCT_SpacePointContainer", "Space point container"}; + + // covariance of the initial parameters + Gaudi::Property<double> m_covLoc0 {this, "covLoc0", 1}; + Gaudi::Property<double> m_covLoc1 {this, "covLoc1", 1}; + Gaudi::Property<double> m_covPhi {this, "covPhi", 1}; + Gaudi::Property<double> m_covTheta {this, "covTheta", 1}; + Gaudi::Property<double> m_covQOverP {this, "covQOverP", 1}; + Gaudi::Property<double> m_covTime {this, "covTime", 1}; +}; + + +inline const std::shared_ptr<const Acts::CurvilinearTrackParameters> + SegmentFitTrackFinderTool::initialTrackParameters() const { + return m_initialTrackParameters; +} + +inline const std::shared_ptr<const Acts::Surface> + SegmentFitTrackFinderTool::initialSurface() const { +return m_initialSurface; +} + +inline const std::shared_ptr<std::vector<IndexSourceLink>> + SegmentFitTrackFinderTool::sourceLinks() const { + return m_sourceLinks; +} + +inline const std::shared_ptr<std::vector<Measurement>> + SegmentFitTrackFinderTool::measurements() const { + return m_measurements; +} + +#endif // FASERACTSKALMANFILTER_SEGMENTFITTRACKFINDERTOOL_H diff --git a/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py b/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py index aeb9be18ebf35e1b155c3c490475ab1a4d17e24c..0f10b935969a7112d69c4d54c496778b3f70513a 100644 --- a/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py +++ b/Tracking/Acts/FaserActsKalmanFilter/python/FaserActsKalmanFilterConfig.py @@ -55,19 +55,25 @@ def FaserActsKalmanFilterCfg(flags, **kwargs): acc = ComponentAccumulator() trajectory_writer_tool = CompFactory.TrajectoryWriterTool() trajectory_writer_tool.FilePath = "KalmanFilterTrajectories.root" - track_finder_tool = CompFactory.TruthSeededTrackFinderTool() + track_finder_tool = CompFactory.SegmentFitTrackFinderTool() + track_finder_tool.covLoc0 = 1e-4 + track_finder_tool.covLoc1 = 1e-4 + track_finder_tool.covPhi = 1e-4 + track_finder_tool.covTheta = 1e-4 + track_finder_tool.covQOverP = 1e-4 - track_finder_tool.useTrueInitialParameters = False - track_finder_tool.covMeas00 = 0.0004 - track_finder_tool.covMeas01 = 0.01 - track_finder_tool.covMeas10 = 0.01 - track_finder_tool.covMeas11 = 0.64 - - track_finder_tool.covLoc0 = 1e-1 - track_finder_tool.covLoc1 = 1e-1 - track_finder_tool.covPhi = 1e-1 - track_finder_tool.covTheta = 1e-1 - track_finder_tool.covQOverP = 1e-1 + # track_finder_tool = CompFactory.TruthSeededTrackFinderTool() + # track_finder_tool.useTrueInitialParameters = False + # track_finder_tool.covMeas00 = 0.0004 + # track_finder_tool.covMeas01 = 0.01 + # track_finder_tool.covMeas10 = 0.01 + # track_finder_tool.covMeas11 = 0.64 + # + # track_finder_tool.covLoc0 = 1e-1 + # track_finder_tool.covLoc1 = 1e-1 + # track_finder_tool.covPhi = 1e-1 + # track_finder_tool.covTheta = 1e-1 + # track_finder_tool.covQOverP = 1e-1 # track_finder_tool.sigmaLoc0 = 0.02 # track_finder_tool.sigmaLoc1 = 0.8 @@ -94,6 +100,6 @@ def FaserActsKalmanFilterCfg(flags, **kwargs): kalman_filter = CompFactory.FaserActsKalmanFilterAlg(**kwargs) kalman_filter.OutputTool = trajectory_writer_tool kalman_filter.TrackFinderTool = track_finder_tool - kalman_filter.ActsLogging = "VERBOSE" + kalman_filter.ActsLogging = "INFO" acc.addEventAlgo(kalman_filter) return acc diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/SegmentFitTrackFinderTool.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/SegmentFitTrackFinderTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..12ee0cae5609e33a477f2724f083a03909b32fad --- /dev/null +++ b/Tracking/Acts/FaserActsKalmanFilter/src/SegmentFitTrackFinderTool.cxx @@ -0,0 +1,149 @@ +#include "FaserActsKalmanFilter/SegmentFitTrackFinderTool.h" + +#include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/Measurement.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "FaserActsKalmanFilter/IndexSourceLink.h" +#include "Identifier/Identifier.h" +#include "TrackerIdentifier/FaserSCT_ID.h" +#include "TrackerReadoutGeometry/SCT_DetectorManager.h" +#include "TrkTrack/Track.h" +#include "TrkTrack/TrackStateOnSurface.h" +#include "TrackerRIO_OnTrack/FaserSCT_ClusterOnTrack.h" +#include "TrackerPrepRawData/FaserSCT_Cluster.h" +#include "TrackerSpacePoint/FaserSCT_SpacePoint.h" +#include "TrackerPrepRawData/FaserSCT_ClusterCollection.h" +#include <random> + + +SegmentFitTrackFinderTool::SegmentFitTrackFinderTool(const std::string& type, const std::string& name, const IInterface* parent) + : base_class(type, name, parent) {} + + +StatusCode SegmentFitTrackFinderTool::initialize() { + ATH_CHECK(detStore()->retrieve(m_idHelper, "FaserSCT_ID")); + ATH_CHECK(detStore()->retrieve(m_detManager, "SCT")); + ATH_CHECK(m_trackingGeometryTool.retrieve()); + ATH_CHECK(m_trackCollection.initialize()); + ATH_CHECK(m_spacePointContainerKey.initialize()); + return StatusCode::SUCCESS; +} + + +StatusCode SegmentFitTrackFinderTool::run() { + SG::ReadHandle<TrackCollection> trackCollection {m_trackCollection}; + ATH_CHECK(trackCollection.isValid()); + + SG::ReadHandle<FaserSCT_SpacePointContainer> spacePointContainer {m_spacePointContainerKey}; + ATH_CHECK(spacePointContainer.isValid()); + + using IdentifierMap = std::map<Identifier, Acts::GeometryIdentifier>; + std::shared_ptr<IdentifierMap> identifierMap = m_trackingGeometryTool->getIdentifierMap(); + + // create vector of cluster identifiers which have been used by any track fit + std::vector<Identifier> trackIds; + std::map<int, Amg::Vector3D> positions; + std::map<int, Amg::Vector3D> directions; + // TODO only take best track (most hits, smallest chi2) + // for now I assume that there is only one track in each station + for (const Trk::Track* track : *trackCollection) { + auto fitParameters = track->trackParameters()->front(); + const Amg::Vector3D fitPosition = fitParameters->position(); + const Amg::Vector3D fitMomentum = fitParameters->momentum(); + Identifier id; + for (const Trk::TrackStateOnSurface* state : *(track->trackStateOnSurfaces())) { + auto clusterOnTrack = dynamic_cast<const Tracker::FaserSCT_ClusterOnTrack*> ( + state->measurementOnTrack()); + if (clusterOnTrack) { + id = clusterOnTrack->identify(); + trackIds.push_back(id); + } + } + int station = m_idHelper->station(id); + positions[station] = fitPosition; + directions[station] = fitMomentum; + } + + // create source links and measurements + const int kSize = 2; + using ThisMeasurement = Acts::Measurement<IndexSourceLink, Acts::BoundIndices, kSize>; + std::array<Acts::BoundIndices, kSize> Indices = {Acts::eBoundLoc0, Acts::eBoundLoc1}; + std::vector<IndexSourceLink> sourceLinks; + std::vector<Measurement> measurements; + for (const FaserSCT_SpacePointCollection* spacePointCollection : *spacePointContainer) { + for (const Tracker::FaserSCT_SpacePoint *spacePoint: *spacePointCollection) { + Identifier id1 = spacePoint->cluster1()->identify(); + Identifier id2 = spacePoint->cluster2()->identify(); + // check if both clusters have been used to reconstruct track + if (std::find(trackIds.begin(), trackIds.end(), id1) != trackIds.end() && + std::find(trackIds.begin(), trackIds.end(), id2) != trackIds.end()) { + Identifier id = spacePoint->associatedSurface().associatedDetectorElementIdentifier(); + // create source link + Acts::GeometryIdentifier geoId = identifierMap->at(id); + IndexSourceLink sourceLink(geoId, measurements.size()); + // create measurement + const auto& par = spacePoint->localParameters(); + const auto& cov = spacePoint->localCovariance(); + ThisMeasurement meas(sourceLink, Indices, par, cov); + sourceLinks.push_back(sourceLink); + measurements.emplace_back(std::move(meas)); + } + } + } + m_sourceLinks = std::make_shared<std::vector<IndexSourceLink>>(sourceLinks); + m_measurements = std::make_shared<std::vector<Measurement>>(measurements); + + // TODO how should we handle events without a track in each station? + if (positions.size() != 3) { + return StatusCode::RECOVERABLE; + } + + // create initial surface + m_initialSurface = Acts::Surface::makeShared<Acts::PlaneSurface>( + Acts::Vector3 {0, 0, 0}, Acts::Vector3{0, 0, 1}); + + // create initial parameters + Acts::Vector3 pos = positions[1] - positions[1].z()/directions[1].z() * directions[1]; + Acts::Vector4 pos4 {pos.x(), pos.y(), pos.z(), 0}; + Acts::Vector3 dir = positions[2] - positions[1]; + Acts::Vector3 tmp {directions[1].x(), directions[1].y(), directions[1].z()}; + ATH_MSG_DEBUG(dir); + ATH_MSG_DEBUG(tmp); + double abs_momentum = momentum(positions); + // TODO get charge from momentum calculation + double charge = 1; + + Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero(); + cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = m_covLoc0; + cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = m_covLoc1; + cov(Acts::eBoundPhi, Acts::eBoundPhi) = m_covPhi; + cov(Acts::eBoundTheta, Acts::eBoundTheta) = m_covTheta; + cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = m_covQOverP; + cov(Acts::eBoundTime, Acts::eBoundTime) = m_covTime; + + auto initialParameters = Acts::CurvilinearTrackParameters( + pos4, dir, abs_momentum, charge, cov); + m_initialTrackParameters = std::make_shared<const Acts::CurvilinearTrackParameters>(initialParameters); + + return StatusCode::SUCCESS; +} + + +StatusCode SegmentFitTrackFinderTool::finalize() { + return StatusCode::SUCCESS; +} + + +double SegmentFitTrackFinderTool::momentum(const std::map<int, Amg::Vector3D>& pos, double B) { + Acts::Vector3 vec_l = pos.at(3) - pos.at(1); + double abs_l = std::sqrt(vec_l.y() * vec_l.y() + vec_l.z() * vec_l.z()); + double t = (pos.at(2).z() - pos.at(1).z()) / (pos.at(3).z() - pos.at(1).z()); + Acts::Vector3 vec_m = pos.at(1) + t * vec_l; + Acts::Vector3 vec_s = pos.at(2) - vec_m; + double abs_s = std::sqrt(vec_s.y() * vec_s.y() + vec_s.z() * vec_s.z()); + double p_yz = 0.3 * abs_l * abs_l * B / (8 * abs_s * 1000); + return p_yz; +} diff --git a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx index fe4f47bcb41697bc0f64992a360dd5d30ea16d61..f02dd329586e5016fc281bdaac1af1a7ef57fcdc 100755 --- a/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx +++ b/Tracking/Acts/FaserActsKalmanFilter/src/components/FaserActsKalmanFilter_entries.cxx @@ -12,6 +12,7 @@ #include "FaserActsKalmanFilter/SimWriterTool.h" #include "FaserActsKalmanFilter/TruthSeededTrackFinderTool.h" #include "FaserActsKalmanFilter/ProtoTrackWriterTool.h" +#include "FaserActsKalmanFilter/SegmentFitTrackFinderTool.h" DECLARE_COMPONENT(FaserActsKalmanFilterAlg) @@ -24,3 +25,4 @@ DECLARE_COMPONENT(TruthTrackFinderTool) DECLARE_COMPONENT(SimWriterTool) DECLARE_COMPONENT(TruthSeededTrackFinderTool) DECLARE_COMPONENT(ProtoTrackWriterTool) +DECLARE_COMPONENT(SegmentFitTrackFinderTool)