diff --git a/Core/include/Acts/Seeding/GNN_DataStorage.hpp b/Core/include/Acts/Seeding/GNN_DataStorage.hpp new file mode 100644 index 00000000000..7759fdecb0d --- /dev/null +++ b/Core/include/Acts/Seeding/GNN_DataStorage.hpp @@ -0,0 +1,319 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Seeding/GNN_Geometry.hpp" + +#include +#include +#include + +namespace Acts { + +constexpr size_t MAX_SEG_PER_NODE = 1000; // was 30 +constexpr size_t N_SEG_CONNS = 6; // was 6 + +// new sp struct +template +struct FTF_SP { + const space_point_t *SP; // want inside to have pointer + int FTF_ID; + int combined_ID; + FTF_SP(const space_point_t *sp, int id, int combined_id) + : SP(sp), FTF_ID(id), combined_ID{combined_id} { + if (SP->sourceLinks().size() == 1) { // pixels have 1 SL + m_isPixel = true; + } else { + m_isPixel = false; + } + m_phi = std::atan(SP->x() / SP->y()); + }; + bool isPixel() const { return m_isPixel; } + bool isSCT() const { return !m_isPixel; } + float phi() const { return m_phi; } + bool m_isPixel; + float m_phi; +}; + +template +class TrigFTF_GNN_Node { + public: + struct CompareByPhi { + bool operator()(const TrigFTF_GNN_Node *n1, + const TrigFTF_GNN_Node *n2) { + return (n1->m_sp_FTF.phi() < n2->m_sp_FTF.phi()); + } + }; + + TrigFTF_GNN_Node(const FTF_SP &FTF_sp, float minT = -100.0, + float maxT = 100.0) + : m_sp_FTF(FTF_sp), m_minCutOnTau(minT), m_maxCutOnTau(maxT) {} + + ~TrigFTF_GNN_Node() {} + + inline void addIn(int i) { + if (m_in.size() < MAX_SEG_PER_NODE) { + m_in.push_back(i); + } + } + + inline void addOut(int i) { + if (m_out.size() < MAX_SEG_PER_NODE) { + m_out.push_back(i); + } + } + + inline bool isConnector() const { + if (m_in.empty() || m_out.empty()) + return false; + return true; + } + + inline bool isFull() const { + if (m_in.size() == MAX_SEG_PER_NODE && m_out.size() == MAX_SEG_PER_NODE) + return true; + else + return false; + } + + const FTF_SP &m_sp_FTF; + + std::vector m_in; // indices of the edges in the edge storage + std::vector m_out; + float m_minCutOnTau, m_maxCutOnTau; +}; + +template +class TrigFTF_GNN_EtaBin { + public: + TrigFTF_GNN_EtaBin() { m_vn.clear(); } + + ~TrigFTF_GNN_EtaBin() { + for (typename std::vector *>::iterator it = + m_vn.begin(); + it != m_vn.end(); ++it) { + delete (*it); + } + } + + void sortByPhi() { + std::sort(m_vn.begin(), m_vn.end(), + typename Acts::TrigFTF_GNN_Node::CompareByPhi()); + } + + bool empty() const { return m_vn.empty(); } + + void generatePhiIndexing(float dphi) { + for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) { + TrigFTF_GNN_Node *pN = m_vn.at(nIdx); + // float phi = pN->m_sp.phi(); + // float phi = (std::atan(pN->m_sp.x() / pN->m_sp.y())); + float phi = pN->m_sp_FTF.phi(); + if (phi <= M_PI - dphi) { + continue; + } + + m_vPhiNodes.push_back( + std::pair(phi - 2 * M_PI, nIdx)); + } + + for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) { + TrigFTF_GNN_Node *pN = m_vn.at(nIdx); + float phi = pN->m_sp_FTF.phi(); + m_vPhiNodes.push_back(std::pair(phi, nIdx)); + } + + for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) { + TrigFTF_GNN_Node *pN = m_vn.at(nIdx); + float phi = pN->m_sp_FTF.phi(); + if (phi >= -M_PI + dphi) { + break; + } + m_vPhiNodes.push_back( + std::pair(phi + 2 * M_PI, nIdx)); + } + } + + std::vector *> m_vn; + // TODO change to + // std::vector>> m_vn; + std::vector> m_vPhiNodes; +}; + +template +class TrigFTF_GNN_DataStorage { + public: + TrigFTF_GNN_DataStorage(const TrigFTF_GNN_Geometry &g) + : m_geo(g) { + m_etaBins.reserve(g.num_bins()); + for (int k = 0; k < g.num_bins(); k++) { + m_etaBins.emplace_back(TrigFTF_GNN_EtaBin()); + } + } + + ~TrigFTF_GNN_DataStorage() {} + + int addSpacePoint(const FTF_SP &sp, bool useClusterWidth) { + const TrigFTF_GNN_Layer *pL = + m_geo.getTrigFTF_GNN_LayerByKey(sp.combined_ID); + + if (pL == nullptr) { + return -1; + } + + int binIndex = pL->getEtaBin(sp.SP->z(), sp.SP->r()); + + if (binIndex == -1) { + return -2; + } + + bool isBarrel = (pL->m_layer.m_type == 0); + + if (isBarrel) { + float min_tau = -100.0; + float max_tau = 100.0; + // can't do this bit yet as dont have cluster width + if (useClusterWidth) { + // const Trk::SpacePoint* osp = sp.offlineSpacePoint(); + // const InDet::PixelCluster* pCL = dynamic_cast(osp->clusterList().first); + // float cluster_width = pCL->width().widthPhiRZ().y(); + float cluster_width = 1; // temporary while cluster width not available + min_tau = 6.7 * (cluster_width - 0.2); + max_tau = + 1.6 + 0.15 / (cluster_width + 0.2) + 6.1 * (cluster_width - 0.2); + } + + m_etaBins.at(binIndex).m_vn.push_back(new TrigFTF_GNN_Node( + sp, min_tau, max_tau)); // adding ftf member to nodes + } else { + if (useClusterWidth) { + // const Trk::SpacePoint* osp = sp.offlineSpacePoint(); + // const InDet::PixelCluster* pCL = dynamic_cast(osp->clusterList().first); + // float cluster_width = pCL->width().widthPhiRZ().y(); + float cluster_width = 1; // temporary while cluster width not available + if (cluster_width > 0.2) { + return -3; + } + } + m_etaBins.at(binIndex).m_vn.push_back( + new TrigFTF_GNN_Node(sp)); + } + + return 0; + } + + // for safety to prevent passing as copy + TrigFTF_GNN_DataStorage(const TrigFTF_GNN_DataStorage &) = delete; + TrigFTF_GNN_DataStorage &operator=(const TrigFTF_GNN_DataStorage &) = delete; + + unsigned int numberOfNodes() const { + unsigned int n = 0; + + for (auto &b : m_etaBins) { + n += b.m_vn.size(); + } + return n; + } + + void getConnectingNodes( + std::vector *> &vn) { + vn.clear(); + vn.reserve(numberOfNodes()); + for (const auto &b : m_etaBins) { + for (typename std::vector< + TrigFTF_GNN_Node *>::const_iterator nIt = + b.m_vn.begin(); + nIt != b.m_vn.end(); ++nIt) { + if ((*nIt)->m_in.empty()) { + continue; + } + if ((*nIt)->m_out.empty()) { + continue; + } + vn.push_back(*nIt); + } + } + } + + void sortByPhi() { + for (auto &b : m_etaBins) { + b.sortByPhi(); + } + } + + void generatePhiIndexing(float dphi) { + for (auto &b : m_etaBins) { + b.generatePhiIndexing(dphi); + } + } + + const TrigFTF_GNN_EtaBin &getEtaBin(int idx) const { + if (idx >= static_cast(m_etaBins.size())) { + idx = idx - 1; + } + return m_etaBins.at(idx); + } + + protected: + const TrigFTF_GNN_Geometry &m_geo; + + std::vector> m_etaBins; +}; + +template +class TrigFTF_GNN_Edge { + public: + struct CompareLevel { + public: + bool operator()(const TrigFTF_GNN_Edge *pS1, const TrigFTF_GNN_Edge *pS2) { + return pS1->m_level > pS2->m_level; + } + }; + + TrigFTF_GNN_Edge(TrigFTF_GNN_Node *n1, + TrigFTF_GNN_Node *n2, float p1, float p2, + float p3, float p4) + : m_n1(n1), m_n2(n2), m_level(1), m_next(1), m_nNei(0) { + m_p[0] = p1; + m_p[1] = p2; + m_p[2] = p3; + m_p[3] = p4; + } + + TrigFTF_GNN_Edge() + : m_n1(nullptr), m_n2(nullptr), m_level(-1), m_next(-1), m_nNei(0){}; + + // TrigFTF_GNN_Edge(const TrigFTF_GNN_Edge &e) + // : m_n1(e.m_n1), m_n2(e.m_n2){}; + + // inline void initialize(TrigFTF_GNN_Node *n1, + // TrigFTF_GNN_Node *n2) { + // m_n1 = n1; + // m_n2 = n2; + // m_level = 1; + // m_next = 1; + // m_nNei = 0; + // } + + TrigFTF_GNN_Node *m_n1{nullptr}; + TrigFTF_GNN_Node *m_n2{nullptr}; + + signed char m_level{-1}, m_next{-1}; + + unsigned char m_nNei{0}; + float m_p[4]{}; + + unsigned int m_vNei[N_SEG_CONNS]{}; // global indices of the connected edges +}; + +} // namespace Acts diff --git a/Core/include/Acts/Seeding/GNN_Geometry.hpp b/Core/include/Acts/Seeding/GNN_Geometry.hpp new file mode 100644 index 00000000000..c7f35145975 --- /dev/null +++ b/Core/include/Acts/Seeding/GNN_Geometry.hpp @@ -0,0 +1,386 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +#include "Acts/TrackFinding/FasTrackConnector.hpp" + +#include +#include +#include +#include +#include + +namespace Acts { +class TrigInDetSiLayer { + public: + int m_subdet; // combined ID + int m_type; // 0: barrel, +/-n : endcap + float m_refCoord; + float m_minBound, m_maxBound; + + TrigInDetSiLayer(int subdet, short int type, float center, float min, + float max) + : m_subdet(subdet), + m_type(type), + m_refCoord(center), + m_minBound(min), + m_maxBound(max) {} +}; + +template +class TrigFTF_GNN_Layer { + public: + TrigFTF_GNN_Layer(const TrigInDetSiLayer &ls, float ew, int bin0) + : m_layer(ls), m_etaBinWidth(ew) { + if (m_layer.m_type == 0) { // barrel + m_r1 = m_layer.m_refCoord; + m_r2 = m_layer.m_refCoord; + m_z1 = m_layer.m_minBound; + m_z2 = m_layer.m_maxBound; + } else { // endcap + m_r1 = m_layer.m_minBound; + m_r2 = m_layer.m_maxBound; + m_z1 = m_layer.m_refCoord; + m_z2 = m_layer.m_refCoord; + } + + float t1 = m_z1 / m_r1; + float eta1 = -std::log(sqrt(1 + t1 * t1) - t1); + + float t2 = m_z2 / m_r2; + float eta2 = -std::log(sqrt(1 + t2 * t2) - t2); + + m_minEta = eta1; + m_maxEta = eta2; + + if (m_maxEta < m_minEta) { + m_minEta = eta2; + m_maxEta = eta1; + } + + m_maxEta += 1e-6; // increasing them slightly to avoid range_check + // exceptions + m_minEta -= 1e-6; + + float deltaEta = m_maxEta - m_minEta; + + int binCounter = bin0; + + if (deltaEta < m_etaBinWidth) { + m_nBins = 1; + m_bins.push_back(binCounter++); + m_etaBin = deltaEta; + if (m_layer.m_type == 0) { // barrel + m_minRadius.push_back(m_layer.m_refCoord - 2.0); + m_maxRadius.push_back(m_layer.m_refCoord + 2.0); + m_minBinCoord.push_back(m_layer.m_minBound); + m_maxBinCoord.push_back(m_layer.m_maxBound); + } else { // endcap + m_minRadius.push_back(m_layer.m_minBound - 2.0); + m_maxRadius.push_back(m_layer.m_maxBound + 2.0); + m_minBinCoord.push_back(m_layer.m_minBound); + m_maxBinCoord.push_back(m_layer.m_maxBound); + } + } else { + float nB = static_cast(deltaEta / m_etaBinWidth); + m_nBins = nB; + if (deltaEta - m_etaBinWidth * nB > 0.5 * m_etaBinWidth) { + m_nBins++; + } + m_etaBin = deltaEta / m_nBins; + + if (m_nBins == 1) { + m_bins.push_back(binCounter++); + if (m_layer.m_type == 0) { // barrel + m_minRadius.push_back(m_layer.m_refCoord - 2.0); + m_maxRadius.push_back(m_layer.m_refCoord + 2.0); + m_minBinCoord.push_back(m_layer.m_minBound); + m_maxBinCoord.push_back(m_layer.m_maxBound); + } else { // endcap + m_minRadius.push_back(m_layer.m_minBound - 2.0); + m_maxRadius.push_back(m_layer.m_maxBound + 2.0); + m_minBinCoord.push_back(m_layer.m_minBound); + m_maxBinCoord.push_back(m_layer.m_maxBound); + } + } else { + float eta = m_minEta + 0.5 * m_etaBin; + + for (int i = 1; i <= m_nBins; i++) { + m_bins.push_back(binCounter++); + + float e1 = eta - 0.5 * m_etaBin; + float e2 = eta + 0.5 * m_etaBin; + + if (m_layer.m_type == 0) { // barrel + m_minRadius.push_back(m_layer.m_refCoord - 2.0); + m_maxRadius.push_back(m_layer.m_refCoord + 2.0); + float z1 = m_layer.m_refCoord * std::sinh(e1); + m_minBinCoord.push_back(z1); + float z2 = m_layer.m_refCoord * std::sinh(e2); + m_maxBinCoord.push_back(z2); + } else { // endcap + float r = m_layer.m_refCoord / std::sinh(e1); + m_minBinCoord.push_back(r); + m_minRadius.push_back(r - 2.0); + r = m_layer.m_refCoord / std::sinh(e2); + m_maxBinCoord.push_back(r); + m_maxRadius.push_back(r + 2.0); + } + + eta += m_etaBin; + } + } + } + } + + int getEtaBin(float zh, float rh) const { + if (m_bins.size() == 1) { + return m_bins.at(0); + } + float t1 = zh / rh; + float eta = -std::log(std::sqrt(1 + t1 * t1) - t1); + + int idx = static_cast((eta - m_minEta) / m_etaBin); + + if (idx < 0) { + idx = 0; + } + if (idx >= static_cast(m_bins.size())) { + idx = static_cast(m_bins.size()) - 1; + } + return m_bins.at(idx); // index in the global storage + } + + float getMinBinRadius(int idx) const { + if (idx >= static_cast(m_minRadius.size())) { + idx = idx - 1; + } + if (idx < 0) { + idx = 0; + } + return m_minRadius.at(idx); + } + + float getMaxBinRadius(int idx) const { + if (idx >= static_cast(m_maxRadius.size())) { + idx = idx - 1; + } + if (idx < 0) { + idx = 0; + } + return m_maxRadius.at(idx); + } + + int num_bins() const { return m_bins.size(); } + + bool verifyBin(const TrigFTF_GNN_Layer *pL, int b1, int b2, + float min_z0, float max_z0) const { + float z1min = m_minBinCoord.at(b1); + float z1max = m_maxBinCoord.at(b1); + float r1 = m_layer.m_refCoord; + + if (m_layer.m_type == 0 && pL->m_layer.m_type == 0) { // barrel -> barrel + + const float tol = 5.0; + + float min_b2 = pL->m_minBinCoord.at(b2); + float max_b2 = pL->m_maxBinCoord.at(b2); + + float r2 = pL->m_layer.m_refCoord; + + float A = r2 / (r2 - r1); + float B = r1 / (r2 - r1); + + float z0_min = z1min * A - max_b2 * B; + float z0_max = z1max * A - min_b2 * B; + + if (z0_max < min_z0 - tol || z0_min > max_z0 + tol) { + return false; + } + return true; + } + + if (m_layer.m_type == 0 && pL->m_layer.m_type != 0) { // barrel -> endcap + + const float tol = 10.0; + + float z2 = pL->m_layer.m_refCoord; + float r2max = pL->m_maxBinCoord.at(b2); + float r2min = pL->m_minBinCoord.at(b2); + + if (r2max <= r1) { + return false; + } + if (r2min <= r1) { + r2min = r1 + 1e-3; + } + + float z0_max = 0.0; + float z0_min = 0.0; + + if (z2 > 0) { + z0_max = (z1max * r2max - z2 * r1) / (r2max - r1); + z0_min = (z1min * r2min - z2 * r1) / (r2min - r1); + } else { + z0_max = (z1max * r2min - z2 * r1) / (r2min - r1); + z0_min = (z1min * r2max - z2 * r1) / (r2max - r1); + } + + if (z0_max < min_z0 - tol || z0_min > max_z0 + tol) { + return false; + } + return true; + } + + return true; + } + + const TrigInDetSiLayer &m_layer; + std::vector m_bins; // eta-bin indices + std::vector m_minRadius; + std::vector m_maxRadius; + std::vector m_minBinCoord; + std::vector m_maxBinCoord; + + float m_minEta, m_maxEta; + + protected: + float m_etaBinWidth, m_phiBinWidth; + + float m_r1, m_z1, m_r2, m_z2; + float m_nBins; + float m_etaBin; +}; + +template +class TrigFTF_GNN_Geometry { + public: + TrigFTF_GNN_Geometry(const std::vector &layers, + std::unique_ptr &conn) + + : m_nEtaBins(0), m_fastrack(move(conn)) { + const float min_z0 = -168.0; + const float max_z0 = 168.0; + + m_etaBinWidth = m_fastrack->m_etaBin; + for (const auto &layer : layers) { + const TrigFTF_GNN_Layer *pL = + addNewLayer(layer, m_nEtaBins); + m_nEtaBins += pL->num_bins(); + } + + // calculating bin tables in the connector... + + for (std::map>::const_iterator it = + m_fastrack->m_connMap.begin(); + it != m_fastrack->m_connMap.end(); ++it) { + const std::vector &vConn = (*it).second; + + for (std::vector::const_iterator cIt = + vConn.begin(); + cIt != vConn.end(); ++cIt) { + unsigned int src = (*cIt)->m_src; // n2 : the new connectors + unsigned int dst = (*cIt)->m_dst; // n1 + + const TrigFTF_GNN_Layer *pL1 = + getTrigFTF_GNN_LayerByKey(dst); + const TrigFTF_GNN_Layer *pL2 = + getTrigFTF_GNN_LayerByKey(src); + + if (pL1 == nullptr) { + std::cout << " skipping invalid dst layer " << dst << std::endl; + continue; + } + if (pL2 == nullptr) { + std::cout << " skipping invalid src layer " << src << std::endl; + continue; + } + int nSrcBins = pL2->m_bins.size(); + int nDstBins = pL1->m_bins.size(); + + (*cIt)->m_binTable.resize(nSrcBins * nDstBins, 0); + + for (int b1 = 0; b1 < nDstBins; b1++) { // loop over bins in Layer 1 + for (int b2 = 0; b2 < nSrcBins; b2++) { // loop over bins in Layer 2 + if (!pL1->verifyBin(pL2, b1, b2, min_z0, max_z0)) { + continue; + } + int address = b1 + b2 * nDstBins; + (*cIt)->m_binTable.at(address) = 1; + } + } + } + } + } + + TrigFTF_GNN_Geometry() : m_nEtaBins(0) {} + + // for safety to prevent passing as copy + TrigFTF_GNN_Geometry(const TrigFTF_GNN_Geometry &) = delete; + TrigFTF_GNN_Geometry &operator=(const TrigFTF_GNN_Geometry &) = delete; + + ~TrigFTF_GNN_Geometry() { + for (typename std::vector *>::iterator it = + m_layArray.begin(); + it != m_layArray.end(); ++it) { + delete (*it); + } + + m_layMap.clear(); + m_layArray.clear(); + } + + const TrigFTF_GNN_Layer *getTrigFTF_GNN_LayerByKey( + unsigned int key) const { + typename std::map *>::const_iterator it = + m_layMap.find(key); + if (it == m_layMap.end()) { + return nullptr; + } + + return (*it).second; + } + + const TrigFTF_GNN_Layer *getTrigFTF_GNN_LayerByIndex( + int idx) const { + return m_layArray.at(idx); + } + + int num_bins() const { return m_nEtaBins; } + + Acts::FasTrackConnector *fastrack() const { return m_fastrack.get(); } + + protected: + const TrigFTF_GNN_Layer *addNewLayer(const TrigInDetSiLayer &l, + int bin0) { + unsigned int layerKey = l.m_subdet; // this should be combined ID + float ew = m_etaBinWidth; + + TrigFTF_GNN_Layer *pHL = + new TrigFTF_GNN_Layer(l, ew, bin0); + + m_layMap.insert(std::pair *>( + layerKey, pHL)); + m_layArray.push_back(pHL); + return pHL; + } + + float m_etaBinWidth; + + std::map *> m_layMap; + std::vector *> m_layArray; + + int m_nEtaBins; + + std::unique_ptr m_fastrack; +}; + +} // namespace Acts diff --git a/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp b/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp new file mode 100644 index 00000000000..6df88c75e84 --- /dev/null +++ b/Core/include/Acts/Seeding/GNN_TrackingFilter.hpp @@ -0,0 +1,387 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +#include "Acts/Seeding/GNN_DataStorage.hpp" //includes geo which has trigindetsilayer, may move this to trigbase + +#include +#include +#include +#include +#include +#include + +template +struct TrigFTF_GNN_EdgeState { + public: + struct Compare { + bool operator()(const struct TrigFTF_GNN_EdgeState* s1, + const struct TrigFTF_GNN_EdgeState* s2) { + return s1->m_J > s2->m_J; + } + }; + + TrigFTF_GNN_EdgeState(){}; + + TrigFTF_GNN_EdgeState(bool f) : m_initialized(f){}; + + ~TrigFTF_GNN_EdgeState(){}; + + void initialize(Acts::TrigFTF_GNN_Edge* pS) { + m_initialized = true; + + m_J = 0.0; + m_vs.clear(); + + // n2->n1 + + float dx = pS->m_n1->m_sp_FTF.SP->x() - pS->m_n2->m_sp_FTF.SP->x(); + float dy = pS->m_n1->m_sp_FTF.SP->y() - pS->m_n2->m_sp_FTF.SP->y(); + float L = std::sqrt(dx * dx + dy * dy); + + m_s = dy / L; + m_c = dx / L; + + // transform for extrapolation and update + // x' = x*m_c + y*m_s + // y' = -x*m_s + y*m_c + + m_refY = pS->m_n2->m_sp_FTF.SP->r(); + m_refX = + pS->m_n2->m_sp_FTF.SP->x() * m_c + pS->m_n2->m_sp_FTF.SP->y() * m_s; + + // X-state: y, dy/dx, d2y/dx2 + + m_X[0] = + -pS->m_n2->m_sp_FTF.SP->x() * m_s + pS->m_n2->m_sp_FTF.SP->y() * m_c; + m_X[1] = 0.0; + m_X[2] = 0.0; + + // Y-state: z, dz/dr + + m_Y[0] = pS->m_n2->m_sp_FTF.SP->z(); + m_Y[1] = (pS->m_n1->m_sp_FTF.SP->z() - pS->m_n2->m_sp_FTF.SP->z()) / + (pS->m_n1->m_sp_FTF.SP->r() - pS->m_n2->m_sp_FTF.SP->r()); + + memset(&m_Cx[0][0], 0, sizeof(m_Cx)); + memset(&m_Cy[0][0], 0, sizeof(m_Cy)); + + m_Cx[0][0] = 0.25; + m_Cx[1][1] = 0.001; + m_Cx[2][2] = 0.001; + + m_Cy[0][0] = 1.5; + m_Cy[1][1] = 0.001; + } + + void clone(const struct TrigFTF_GNN_EdgeState& st) { + memcpy(&m_X[0], &st.m_X[0], sizeof(m_X)); + memcpy(&m_Y[0], &st.m_Y[0], sizeof(m_Y)); + memcpy(&m_Cx[0][0], &st.m_Cx[0][0], sizeof(m_Cx)); + memcpy(&m_Cy[0][0], &st.m_Cy[0][0], sizeof(m_Cy)); + m_refX = st.m_refX; + m_refY = st.m_refY; + m_c = st.m_c; + m_s = st.m_s; + m_J = st.m_J; + m_vs.clear(); + m_vs.reserve(st.m_vs.size()); + std::copy(st.m_vs.begin(), st.m_vs.end(), std::back_inserter(m_vs)); + + m_initialized = true; + } + + float m_J{}; + + std::vector*> m_vs; + + float m_X[3]{}, m_Y[2]{}, m_Cx[3][3]{}, m_Cy[2][2]{}; + float m_refX{}, m_refY{}, m_c{}, m_s{}; + + bool m_initialized{false}; +}; + +#define MAX_EDGE_STATE 2500 + +template +class TrigFTF_GNN_TrackingFilter { + public: + TrigFTF_GNN_TrackingFilter( + const std::vector& g, + std::vector>& sb) + : m_geo(g), m_segStore(sb) {} + ~TrigFTF_GNN_TrackingFilter(){}; + + void followTrack(Acts::TrigFTF_GNN_Edge* pS, + TrigFTF_GNN_EdgeState& output) { + if (pS->m_level == -1) { + return; // already collected + } + m_globalStateCounter = 0; + + // create track state + + TrigFTF_GNN_EdgeState* pInitState = + &m_stateStore[m_globalStateCounter++]; + + pInitState->initialize(pS); + + m_stateVec.clear(); + + // recursive branching and propagation + + propagate(pS, *pInitState); + + if (m_stateVec.empty()) { + return; + } + std::sort(m_stateVec.begin(), m_stateVec.end(), + typename TrigFTF_GNN_EdgeState::Compare()); + + TrigFTF_GNN_EdgeState* best = (*m_stateVec.begin()); + + output.clone(*best); + + m_globalStateCounter = 0; + } + + protected: + void propagate(Acts::TrigFTF_GNN_Edge* pS, + TrigFTF_GNN_EdgeState& ts) { + if (m_globalStateCounter >= MAX_EDGE_STATE) { + return; + } + TrigFTF_GNN_EdgeState* p_new_ts = + &m_stateStore[m_globalStateCounter++]; + + TrigFTF_GNN_EdgeState& new_ts = *p_new_ts; + new_ts.clone(ts); + + new_ts.m_vs.push_back(pS); + + bool accepted = update(pS, new_ts); // update using n1 of the segment + + if (!accepted) { + return; // stop further propagation + } + int level = pS->m_level; + + std::list*> lCont; + + for (int nIdx = 0; nIdx < pS->m_nNei; + nIdx++) { // loop over the neighbours of this segment + unsigned int nextSegmentIdx = pS->m_vNei[nIdx]; + + Acts::TrigFTF_GNN_Edge* pN = + &(m_segStore.at(nextSegmentIdx)); + + if (pN->m_level == -1) { + continue; // already collected + } + if (pN->m_level == level - 1) { + lCont.push_back(pN); + } + } + if (lCont.empty()) { // the end of chain + + // store in the vector + if (m_globalStateCounter < MAX_EDGE_STATE) { + if (m_stateVec.empty()) { // add the first segment state + TrigFTF_GNN_EdgeState* p = + &m_stateStore[m_globalStateCounter++]; + p->clone(new_ts); + m_stateVec.push_back(p); + } else { // compare with the best and add + float best_so_far = (*m_stateVec.begin())->m_J; + if (new_ts.m_J > best_so_far) { + TrigFTF_GNN_EdgeState* p = + &m_stateStore[m_globalStateCounter++]; + p->clone(new_ts); + m_stateVec.push_back(p); + } + } + } + } else { // branching + int nBranches = 0; + for (typename std::list< + Acts::TrigFTF_GNN_Edge*>::iterator sIt = + lCont.begin(); + sIt != lCont.end(); ++sIt, nBranches++) { + propagate((*sIt), new_ts); // recursive call + } + } + } + + bool update(Acts::TrigFTF_GNN_Edge* pS, + TrigFTF_GNN_EdgeState& ts) { + const float sigma_t = 0.0003; + const float sigma_w = 0.00009; + + const float sigmaMS = 0.016; + + const float sigma_x = 0.25; // was 0.22 + const float sigma_y = 2.5; // was 1.7 + + const float weight_x = 0.5; + const float weight_y = 0.5; + + const float maxDChi2_x = 60.0; // 35.0; + const float maxDChi2_y = 60.0; // 31.0; + + const float add_hit = 14.0; + + if (ts.m_Cx[2][2] < 0.0 || ts.m_Cx[1][1] < 0.0 || ts.m_Cx[0][0] < 0.0) { + std::cout << "Negative cov_x" << std::endl; + } + + if (ts.m_Cy[1][1] < 0.0 || ts.m_Cy[0][0] < 0.0) { + std::cout << "Negative cov_y" << std::endl; + } + + // add ms. + + ts.m_Cx[2][2] += sigma_w * sigma_w; + ts.m_Cx[1][1] += sigma_t * sigma_t; + + int type1 = getLayerType(pS->m_n1->m_sp_FTF.combined_ID); + + float t2 = type1 == 0 ? 1.0 + ts.m_Y[1] * ts.m_Y[1] + : 1.0 + 1.0 / (ts.m_Y[1] * ts.m_Y[1]); + float s1 = sigmaMS * t2; + float s2 = s1 * s1; + + s2 *= std::sqrt(t2); + + ts.m_Cy[1][1] += s2; + + // extrapolation + + float X[3], Y[2]; + float Cx[3][3], Cy[2][2]; + + float refX, refY, mx, my; + + float x, y, z, r; + + x = pS->m_n1->m_sp_FTF.SP->x(); + y = pS->m_n1->m_sp_FTF.SP->y(); + z = pS->m_n1->m_sp_FTF.SP->z(); + r = pS->m_n1->m_sp_FTF.SP->r(); + + refX = x * ts.m_c + y * ts.m_s; + mx = -x * ts.m_s + y * ts.m_c; // measured X[0] + refY = r; + my = z; // measured Y[0] + + float A = refX - ts.m_refX; + float B = 0.5 * A * A; + float dr = refY - ts.m_refY; + + X[0] = ts.m_X[0] + ts.m_X[1] * A + ts.m_X[2] * B; + X[1] = ts.m_X[1] + ts.m_X[2] * A; + X[2] = ts.m_X[2]; + + Cx[0][0] = ts.m_Cx[0][0] + 2 * ts.m_Cx[0][1] * A + 2 * ts.m_Cx[0][2] * B + + A * A * ts.m_Cx[1][1] + 2 * A * B * ts.m_Cx[1][2] + + B * B * ts.m_Cx[2][2]; + Cx[0][1] = Cx[1][0] = ts.m_Cx[0][1] + ts.m_Cx[1][1] * A + + ts.m_Cx[1][2] * B + ts.m_Cx[0][2] * A + + A * A * ts.m_Cx[1][2] + A * B * ts.m_Cx[2][2]; + Cx[0][2] = Cx[2][0] = ts.m_Cx[0][2] + ts.m_Cx[1][2] * A + ts.m_Cx[2][2] * B; + + Cx[1][1] = ts.m_Cx[1][1] + 2 * A * ts.m_Cx[1][2] + A * A * ts.m_Cx[2][2]; + Cx[1][2] = Cx[2][1] = ts.m_Cx[1][2] + ts.m_Cx[2][2] * A; + + Cx[2][2] = ts.m_Cx[2][2]; + + Y[0] = ts.m_Y[0] + ts.m_Y[1] * dr; + Y[1] = ts.m_Y[1]; + + Cy[0][0] = ts.m_Cy[0][0] + 2 * ts.m_Cy[0][1] * dr + dr * dr * ts.m_Cy[1][1]; + Cy[0][1] = Cy[1][0] = ts.m_Cy[0][1] + dr * ts.m_Cy[1][1]; + Cy[1][1] = ts.m_Cy[1][1]; + + // chi2 test + + float resid_x = mx - X[0]; + float resid_y = my - Y[0]; + + float CHx[3] = {Cx[0][0], Cx[0][1], Cx[0][2]}; + float CHy[2] = {Cy[0][0], Cy[0][1]}; + + float sigma_rz = 0.0; + + int type = getLayerType(pS->m_n1->m_sp_FTF.combined_ID); + + if (type == 0) { // barrel TO-DO: split into barrel Pixel and barrel SCT + sigma_rz = sigma_y * sigma_y; + } else { + sigma_rz = sigma_y * ts.m_Y[1]; + sigma_rz = sigma_rz * sigma_rz; + } + + float Dx = 1.0 / (Cx[0][0] + sigma_x * sigma_x); + + float Dy = 1.0 / (Cy[0][0] + sigma_rz); + + float dchi2_x = resid_x * resid_x * Dx; + float dchi2_y = resid_y * resid_y * Dy; + + if (dchi2_x > maxDChi2_x || dchi2_y > maxDChi2_y) { + return false; + } + + ts.m_J += add_hit - dchi2_x * weight_x - dchi2_y * weight_y; + + // state update + + float Kx[3] = {Dx * Cx[0][0], Dx * Cx[0][1], Dx * Cx[0][2]}; + float Ky[2] = {Dy * Cy[0][0], Dy * Cy[0][1]}; + + for (int i = 0; i < 3; i++) + ts.m_X[i] = X[i] + Kx[i] * resid_x; + for (int i = 0; i < 2; i++) + ts.m_Y[i] = Y[i] + Ky[i] * resid_y; + + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + ts.m_Cx[i][j] = Cx[i][j] - Kx[i] * CHx[j]; + } + } + + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + ts.m_Cy[i][j] = Cy[i][j] - Ky[i] * CHy[j]; + } + } + ts.m_refX = refX; + ts.m_refY = refY; + return true; + } + + int getLayerType(int l) { + auto iterator = find_if(m_geo.begin(), m_geo.end(), [l](auto n) { + return n.m_subdet == l; + }); // iterator to vector member with this id + int index = std::distance(m_geo.begin(), iterator); + + return m_geo.at(index).m_type; + } + + const std::vector& m_geo; + + std::vector>& m_segStore; + + std::vector*> m_stateVec; + + TrigFTF_GNN_EdgeState m_stateStore[MAX_EDGE_STATE]; + + int m_globalStateCounter{0}; +}; diff --git a/Core/include/Acts/Seeding/SeedFinderFTF.hpp b/Core/include/Acts/Seeding/SeedFinderFTF.hpp new file mode 100644 index 00000000000..5d407254022 --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderFTF.hpp @@ -0,0 +1,102 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style + +#include "Acts/EventData/SpacePointData.hpp" +#include "Acts/Seeding/InternalSeed.hpp" +#include "Acts/Seeding/InternalSpacePoint.hpp" +#include "Acts/Seeding/SeedFinderConfig.hpp" +#include "Acts/Seeding/SeedFinderFTFConfig.hpp" +#include "Acts/TrackFinding/RoiDescriptor.hpp" +#include "Acts/Utilities/KDTree.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Acts { + +template +struct GNN_TrigTracklet { + public: + GNN_TrigTracklet(std::vector *> &vSP, + std::vector> &tbuf) + : m_track(vSP), m_seeds(tbuf){}; + ~GNN_TrigTracklet(){}; + + std::vector *> m_track; + std::vector> m_seeds; +}; + +template +class SeedFinderFTF { + public: + static constexpr std::size_t NDims = 3; + + using seed_t = Seed; + // using internal_sp_t = InternalSpacePoint; + // using tree_t = KDTree; + + // constructors + SeedFinderFTF(const SeedFinderFTFConfig &config, + const TrigFTF_GNN_Geometry &GNNgeo); + + ~SeedFinderFTF(); //!!! is it dangerous not to use default? got def in ipp + SeedFinderFTF() = default; + SeedFinderFTF(const SeedFinderFTF &) = delete; + SeedFinderFTF &operator=( + const SeedFinderFTF &) = delete; + + void loadSpacePoints(const std::vector> &); + + void createSeeds(const Acts::RoiDescriptor &, + const Acts::TrigFTF_GNN_Geometry &); + + // create seeeds function + template + void createSeeds_old(const Acts::SeedFinderOptions &options, + const input_container_t &spacePoints, + output_container_t &out_cont, + callable_t &&extract_coordinates) const; + + template + std::vector createSeeds_old(const Acts::SeedFinderOptions &options, + const input_container_t &spacePoints, + callable_t &&extract_coordinates) const; + + private: + enum Dim { DimPhi = 0, DimR = 1, DimZ = 2 }; + + // config object + SeedFinderFTFConfig m_config; + + void runGNN_TrackFinder( + std::vector> &, + const Acts::RoiDescriptor &, + const Acts::TrigFTF_GNN_Geometry &); + + // needs to be member of class so can accessed by all member functions + TrigFTF_GNN_DataStorage *m_storage; + + // for create seeds: + std::vector> m_triplets; +}; + +} // namespace Acts + +#include "Acts/Seeding/SeedFinderFTF.ipp" diff --git a/Core/include/Acts/Seeding/SeedFinderFTF.ipp b/Core/include/Acts/Seeding/SeedFinderFTF.ipp new file mode 100644 index 00000000000..af7a54bb91b --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderFTF.ipp @@ -0,0 +1,731 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// SeedFinderFTF.ipp +// TODO: update to C++17 style + +#include "Acts/Definitions/Algebra.hpp" //for M_PI +#include "Acts/Geometry/Extent.hpp" +#include "Acts/Seeding/SeedFilter.hpp" +#include "Acts/Seeding/SeedFinder.hpp" +#include "Acts/Seeding/SeedFinderFTFConfig.hpp" +#include "Acts/Seeding/SeedFinderUtils.hpp" +#include "Acts/Utilities/BinningType.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +// core so in ACTS namespace + +namespace Acts { + +template +SeedFinderFTF::SeedFinderFTF( + const SeedFinderFTFConfig& config, + const TrigFTF_GNN_Geometry& GNNgeo) + : m_config(config) { + m_storage = new TrigFTF_GNN_DataStorage(GNNgeo); +} + +template +SeedFinderFTF::~SeedFinderFTF() { + delete m_storage; + + m_storage = nullptr; +} + +// define loadspace points function +template +void SeedFinderFTF::loadSpacePoints( + const std::vector>& FTF_SP_vect) { + for (const auto& FTF_sp : FTF_SP_vect) { + bool is_Pixel = FTF_sp.isPixel(); + if (!is_Pixel) { + continue; + } + m_storage->addSpacePoint(FTF_sp, (m_config.m_useClusterWidth > 0)); + } + + m_config.m_nMaxPhiSlice = 1; + m_config.m_phiSliceWidth = 2 * M_PI / m_config.m_nMaxPhiSlice; + + m_storage->sortByPhi(); + + m_storage->generatePhiIndexing(1.5 * m_config.m_phiSliceWidth); +} + +template +void SeedFinderFTF::runGNN_TrackFinder( + std::vector>& vTracks, + const Acts::RoiDescriptor& roi, + const Acts::TrigFTF_GNN_Geometry& gnngeo) { + // long term move these to ftf finder config, then m_config. to access them + const int MaxEdges = 2000000; + + const float cut_dphi_max = 0.012; + const float cut_dcurv_max = 0.001; + const float cut_tau_ratio_max = 0.007; + const float min_z0 = -2800; // roiDescriptor->zedMinus(); //blank for now, + // get from config eventually + const float max_z0 = 2800; // roiDescriptor->zedPlus(); + + const float maxOuterRadius = 550.0; + const float cut_zMinU = + min_z0 + + maxOuterRadius * roi.dzdrMinus(); // dzdr can only find =0 in athena + const float cut_zMaxU = max_z0 + maxOuterRadius * roi.dzdrPlus(); + + float m_minR_squ = 1; // set earlier + float m_maxCurv = 1; + + const float maxKappa_high_eta = 0.8 / m_minR_squ; + const float maxKappa_low_eta = 0.6 / m_minR_squ; + + // 1. loop over stages + + int currentStage = 0; + + const Acts::FasTrackConnector& fastrack = *(gnngeo.fastrack()); + + std::vector> edgeStorage; + + edgeStorage.reserve(MaxEdges); + + int nEdges = 0; + + for (std::map>::const_iterator + it = fastrack.m_layerGroups.begin(); + it != fastrack.m_layerGroups.end(); ++it, currentStage++) { + // loop over L1 layers for the current stage + + for (const auto& layerGroup : (*it).second) { + unsigned int dst = layerGroup.m_dst; // n1 : inner nodes + + const TrigFTF_GNN_Layer* pL1 = + gnngeo.getTrigFTF_GNN_LayerByKey(dst); + + if (pL1 == nullptr) { + continue; + } + + for (const auto& conn : + layerGroup.m_sources) { // loop over L2(L1) for the current stage + + unsigned int src = conn->m_src; // n2 : the new connectors + + const TrigFTF_GNN_Layer* pL2 = + gnngeo.getTrigFTF_GNN_LayerByKey(src); + + if (pL2 == nullptr) { + continue; + } + int nDstBins = pL1->m_bins.size(); + int nSrcBins = pL2->m_bins.size(); + + for (int b1 = 0; b1 < nDstBins; b1++) { // loop over bins in Layer 1 + + const TrigFTF_GNN_EtaBin& B1 = + m_storage->getEtaBin(pL1->m_bins.at(b1)); + + if (B1.empty()) { + continue; + } + + float rb1 = pL1->getMinBinRadius(b1); + + // 3. loops over source eta-bins + + for (int b2 = 0; b2 < nSrcBins; b2++) { // loop over bins in Layer 2 + + if (m_config.m_useEtaBinning && (nSrcBins + nDstBins > 2)) { + if (conn->m_binTable[b1 + b2 * nDstBins] != 1) { + continue; // using precomputed LUT + } + } + + const TrigFTF_GNN_EtaBin& B2 = + m_storage->getEtaBin(pL2->m_bins.at(b2)); + + if (B2.empty()) { + continue; + } + + float rb2 = pL2->getMaxBinRadius(b2); + + // calculated delta Phi for rb1 ---> rb2 extrapolation + + float deltaPhi = + 0.5f * + m_config + .m_phiSliceWidth; // the default sliding window along phi + + if (m_config.m_useEtaBinning) { + deltaPhi = 0.001f + m_maxCurv * std::fabs(rb2 - rb1); + } + + unsigned int first_it = 0; + for (typename std::vector< + TrigFTF_GNN_Node*>::const_iterator + n1It = B1.m_vn.begin(); + n1It != B1.m_vn.end(); ++n1It) { // loop over nodes in Layer 1 + + TrigFTF_GNN_Node* n1 = (*n1It); + + if (n1->m_in.size() >= MAX_SEG_PER_NODE) { + continue; + } + + float r1 = n1->m_sp_FTF.SP->r(); + float x1 = n1->m_sp_FTF.SP->x(); + float y1 = n1->m_sp_FTF.SP->y(); + float z1 = n1->m_sp_FTF.SP->z(); + float phi1 = std::atan(x1 / y1); + + float minPhi = phi1 - deltaPhi; + float maxPhi = phi1 + deltaPhi; + + for (unsigned int n2PhiIdx = first_it; + n2PhiIdx < B2.m_vPhiNodes.size(); + n2PhiIdx++) { // sliding window over nodes in Layer 2 + + float phi2 = B2.m_vPhiNodes.at(n2PhiIdx).first; + + if (phi2 < minPhi) { + first_it = n2PhiIdx; + continue; + } + if (phi2 > maxPhi) { + break; + } + + TrigFTF_GNN_Node* n2 = + B2.m_vn.at(B2.m_vPhiNodes.at(n2PhiIdx).second); + + if (n2->m_out.size() >= MAX_SEG_PER_NODE) { + continue; + } + if (n2->isFull()) { + continue; + } + + float r2 = n2->m_sp_FTF.SP->r(); + + float dr = r2 - r1; + + if (dr < m_config.m_minDeltaRadius) { + continue; + } + + float z2 = n2->m_sp_FTF.SP->z(); + + float dz = z2 - z1; + float tau = dz / dr; + float ftau = std::fabs(tau); + if (ftau > 36.0) { + continue; + } + + if (ftau < n1->m_minCutOnTau) { + continue; + } + if (ftau < n2->m_minCutOnTau) { + continue; + } + if (ftau > n1->m_maxCutOnTau) { + continue; + } + if (ftau > n2->m_maxCutOnTau) { + continue; + } + + if (m_config.m_doubletFilterRZ) { + float z0 = z1 - r1 * tau; + + if (z0 < min_z0 || z0 > max_z0) { + continue; + } + + float zouter = z0 + maxOuterRadius * tau; + + if (zouter < cut_zMinU || zouter > cut_zMaxU) { + continue; + } + } + + float dx = n2->m_sp_FTF.SP->x() - x1; + float dy = n2->m_sp_FTF.SP->y() - y1; + + float L2 = 1 / (dx * dx + dy * dy); + + float D = + (n2->m_sp_FTF.SP->y() * x1 - y1 * n2->m_sp_FTF.SP->x()) / + (r1 * r2); + + float kappa = D * D * L2; + + if (ftau < 4.0) { // eta = 2.1 + if (kappa > maxKappa_low_eta) { + continue; + } + + } else { + if (kappa > maxKappa_high_eta) { + continue; + } + } + + // match edge candidate against edges incoming to n2 + + float exp_eta = std::sqrt(1 + tau * tau) - tau; + + bool isGood = + n2->m_in.size() <= + 2; // we must have enough incoming edges to decide + + if (!isGood) { + float uat_1 = 1.0f / exp_eta; + + for (const auto& n2_in_idx : n2->m_in) { + float tau2 = edgeStorage.at(n2_in_idx).m_p[0]; + float tau_ratio = tau2 * uat_1 - 1.0f; + + if (std::fabs(tau_ratio) > cut_tau_ratio_max) { // bad + // match + continue; + } + isGood = true; // good match found + break; + } + } + if (!isGood) { + continue; // no moatch found, skip creating [n1 <- n2] edge + } + + float curv = D * std::sqrt(L2); // signed curvature + float dPhi2 = std::asin(curv * r2); + float dPhi1 = std::asin(curv * r1); + + if (nEdges < MaxEdges) { + edgeStorage.emplace_back(n1, n2, exp_eta, curv, phi1 + dPhi1, + phi2 + dPhi2); + + n1->addIn(nEdges); + n2->addOut(nEdges); + + nEdges++; + } + } // loop over n2 (outer) nodes + } // loop over n1 (inner) nodes + } // loop over source eta bins + } // loop over dst eta bins + } // loop over L2(L1) layers + } // loop over dst layers + } // loop over the stages of doublet making + + std::vector*> vNodes; + + m_storage->getConnectingNodes(vNodes); + + if (vNodes.empty()) { + return; + } + + int nNodes = vNodes.size(); + + for (int nodeIdx = 0; nodeIdx < nNodes; nodeIdx++) { + const TrigFTF_GNN_Node* pN = vNodes.at(nodeIdx); + + std::vector> in_sort, out_sort; + in_sort.resize(pN->m_in.size()); + out_sort.resize(pN->m_out.size()); + + for (int inIdx = 0; inIdx < static_cast(pN->m_in.size()); inIdx++) { + int inEdgeIdx = pN->m_in.at(inIdx); + Acts::TrigFTF_GNN_Edge* pS = + &(edgeStorage.at(inEdgeIdx)); + in_sort[inIdx].second = inEdgeIdx; + in_sort[inIdx].first = pS->m_p[0]; + } + for (int outIdx = 0; outIdx < static_cast(pN->m_out.size()); + outIdx++) { + int outEdgeIdx = pN->m_out.at(outIdx); + Acts::TrigFTF_GNN_Edge* pS = + &(edgeStorage.at(outEdgeIdx)); + out_sort[outIdx].second = outEdgeIdx; + out_sort[outIdx].first = pS->m_p[0]; + } + + std::sort(in_sort.begin(), in_sort.end()); + std::sort(out_sort.begin(), out_sort.end()); + + unsigned int last_out = 0; + + for (unsigned int in_idx = 0; in_idx < in_sort.size(); + in_idx++) { // loop over incoming edges + + int inEdgeIdx = in_sort[in_idx].second; + + Acts::TrigFTF_GNN_Edge* pS = + &(edgeStorage.at(inEdgeIdx)); + + pS->m_nNei = 0; + float tau1 = pS->m_p[0]; + float uat_1 = 1.0f / tau1; + float curv1 = pS->m_p[1]; + float Phi1 = pS->m_p[2]; + + for (unsigned int out_idx = last_out; out_idx < out_sort.size(); + out_idx++) { + int outEdgeIdx = out_sort[out_idx].second; + + Acts::TrigFTF_GNN_Edge* pNS = + &(edgeStorage.at(outEdgeIdx)); + + float tau2 = pNS->m_p[0]; + float tau_ratio = tau2 * uat_1 - 1.0f; + + if (tau_ratio < -cut_tau_ratio_max) { + last_out = out_idx; + continue; + } + if (tau_ratio > cut_tau_ratio_max) { + break; + } + + float dPhi = pNS->m_p[3] - Phi1; + + if (dPhi < -M_PI) { + dPhi += 2 * M_PI; + } else if (dPhi > M_PI) { + dPhi -= 2 * M_PI; + } + + if (dPhi < -cut_dphi_max || dPhi > cut_dphi_max) { + continue; + } + + float curv2 = pNS->m_p[1]; + float dcurv = curv2 - curv1; + + if (dcurv < -cut_dcurv_max || dcurv > cut_dcurv_max) { + continue; + } + + pS->m_vNei[pS->m_nNei++] = outEdgeIdx; + if (pS->m_nNei >= N_SEG_CONNS) { + break; + } + } + } + } + + const int maxIter = 15; + + int maxLevel = 0; + + int iter = 0; + + std::vector*> v_old; + + for (int edgeIndex = 0; edgeIndex < nEdges; edgeIndex++) { + Acts::TrigFTF_GNN_Edge* pS = + &(edgeStorage.at(edgeIndex)); + if (pS->m_nNei == 0) { + continue; + } + v_old.push_back(pS); // TO-DO: increment level for segments as they already + // have at least one neighbour + } + + for (; iter < maxIter; iter++) { + // generate proposals + std::vector*> v_new; + v_new.clear(); + + for (auto pS : v_old) { + int next_level = pS->m_level; + + for (int nIdx = 0; nIdx < pS->m_nNei; nIdx++) { + unsigned int nextEdgeIdx = pS->m_vNei[nIdx]; + + Acts::TrigFTF_GNN_Edge* pN = + &(edgeStorage.at(nextEdgeIdx)); + + if (pS->m_level == pN->m_level) { + next_level = pS->m_level + 1; + v_new.push_back(pS); + break; + } + } + + pS->m_next = next_level; // proposal + } + // update + + int nChanges = 0; + + for (auto pS : v_new) { + if (pS->m_next != pS->m_level) { + nChanges++; + pS->m_level = pS->m_next; + if (maxLevel < pS->m_level) + maxLevel = pS->m_level; + } + } + + if (nChanges == 0) { + break; + } + + v_old = std::move(v_new); + v_new.clear(); + } + + int minLevel = 3; // a triplet + 2 confirmation + + std::vector*> vSeeds; + + vSeeds.reserve(MaxEdges / 2); + + for (int edgeIndex = 0; edgeIndex < nEdges; edgeIndex++) { + Acts::TrigFTF_GNN_Edge* pS = + &(edgeStorage.at(edgeIndex)); + + if (pS->m_level < minLevel) { + continue; + } + + vSeeds.push_back(pS); + } + + m_triplets.clear(); + + std::sort( + vSeeds.begin(), vSeeds.end(), + typename Acts::TrigFTF_GNN_Edge::CompareLevel()); + + if (vSeeds.empty()) { + return; + } + + // backtracking + + TrigFTF_GNN_TrackingFilter tFilter( + m_config.m_layerGeometry, edgeStorage); + + for (auto pS : vSeeds) { + if (pS->m_level == -1) { + continue; + } + + TrigFTF_GNN_EdgeState rs(false); + + tFilter.followTrack(pS, rs); + + if (!rs.m_initialized) { + continue; + } + + if (static_cast(rs.m_vs.size()) < minLevel) { + continue; + } + + std::vector*> vSP; + + for (typename std::vector*>:: + reverse_iterator sIt = rs.m_vs.rbegin(); + sIt != rs.m_vs.rend(); ++sIt) { + (*sIt)->m_level = -1; // mark as collected + + if (sIt == rs.m_vs.rbegin()) { + vSP.push_back(&(*sIt)->m_n1->m_sp_FTF); + } + vSP.push_back(&(*sIt)->m_n2->m_sp_FTF); + } + + if (vSP.size() < 3) { + continue; + } + + // making triplets + + unsigned int nTriplets = 0; + + std::vector> output; + + for (unsigned int idx_m = 1; idx_m < vSP.size() - 1; idx_m++) { + const FTF_SP& spM = *vSP.at(idx_m); + const double pS_r = spM.SP->r(); + const double pS_x = spM.SP->x(); + const double pS_y = spM.SP->y(); + const double cosA = pS_x / pS_r; + const double sinA = pS_y / pS_r; + + for (unsigned int idx_o = idx_m + 1; idx_o < vSP.size(); idx_o++) { + const FTF_SP& spO = *vSP.at(idx_o); + + double dx = spO.SP->x() - pS_x; + double dy = spO.SP->y() - pS_y; + double R2inv = 1.0 / (dx * dx + dy * dy); + double xn = dx * cosA + dy * sinA; + double yn = -dx * sinA + dy * cosA; + + const double uo = xn * R2inv; + const double vo = yn * R2inv; + + for (unsigned int idx_i = 0; idx_i < idx_m; idx_i++) { + const FTF_SP& spI = *vSP.at(idx_i); + + dx = spI.SP->x() - pS_x; + dy = spI.SP->y() - pS_y; + R2inv = 1.0 / (dx * dx + dy * dy); + + xn = dx * cosA + dy * sinA; + yn = -dx * sinA + dy * cosA; + + const double ui = xn * R2inv; + const double vi = yn * R2inv; + + // 1. pT estimate + + const double du = uo - ui; + if (du == 0.0) { + continue; + } + const double A = (vo - vi) / du; + const double B = vi - A * ui; + const double R_squ = (1 + A * A) / (B * B); + + if (R_squ < m_minR_squ) { + continue; + } + + // 2. d0 cut + + const double fabs_d0 = std::abs(pS_r * (B * pS_r - A)); + + if (fabs_d0 > m_config.m_tripletD0Max) { + continue; + } + + // 3. phi0 cut + + // if (!roi.isFullscan()) { + // const double uc = 2 * B * pS_r - A; + // // const double phi0 = std::atan2(sinA - uc * cosA, cosA + uc * + // sinA); + // // if ( !RoiUtil::containsPhi( *roiDescriptor, phi0 ) ) + // { + // // continue; + // // } + // } + + // 4. add new triplet + + const double Q = fabs_d0 * fabs_d0; + + output.emplace_back(spI, spM, spO, Q); + + nTriplets++; + + if (nTriplets >= m_config.m_maxTripletBufferLength) { + break; + } + } + if (nTriplets >= m_config.m_maxTripletBufferLength) { + break; + } + } + if (nTriplets >= m_config.m_maxTripletBufferLength) { + break; + } + } + + if (output.empty()) { + continue; + } + + vTracks.emplace_back(vSP, output); + } +} + +template +void SeedFinderFTF::createSeeds( + const Acts::RoiDescriptor& roi, + const Acts::TrigFTF_GNN_Geometry& gnngeo) { + std::vector> + vTracks; // make empty vector + + vTracks.reserve(5000); + + runGNN_TrackFinder(vTracks, roi, gnngeo); // returns filled vector + + if (vTracks.empty()) { + return; + } + + m_triplets.clear(); // member of class , saying not declared, maybe public? + + for (auto& track : vTracks) { + for (auto& seed : track.m_seeds) { // access mmeber of GNN_TrigTracklet + + float newQ = seed.Q(); // function of TrigInDetTriplet + if (m_config.m_LRTmode) { + // In LRT mode penalize pixels in Triplets + if (seed.s1().isPixel()) { + newQ += 1000; // functions of TrigSiSpacePointBase + } + if (seed.s2().isPixel()) { + newQ += 1000; + } + if (seed.s3().isPixel()) { + newQ += 1000; + } + } else { + // In normal (non LRT) mode penalise SSS by 1000, PSS (if enabled) and + // PPS by 10000 + if (seed.s3().isSCT()) { + newQ += seed.s1().isSCT() ? 1000.0 : 10000.0; + } + } + seed.Q(newQ); + m_triplets.emplace_back(seed); + } + } + vTracks.clear(); +} + +// // still to be developed +template +template +void SeedFinderFTF::createSeeds_old( + const Acts::SeedFinderOptions& /*options*/, + const input_container_t& /*spacePoints*/, output_container_t& /*out_cont*/, + callable_t&& /*extract_coordinates*/) const {} + +template +template +std::vector> +SeedFinderFTF::createSeeds_old( + const Acts::SeedFinderOptions& options, + const input_container_t& spacePoints, + callable_t&& extract_coordinates) const { + std::vector r; + createSeeds_old(options, spacePoints, r, + std::forward(extract_coordinates)); + return r; +} + +} // namespace Acts diff --git a/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp b/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp new file mode 100644 index 00000000000..8738a01af1f --- /dev/null +++ b/Core/include/Acts/Seeding/SeedFinderFTFConfig.hpp @@ -0,0 +1,90 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/Seeding/SeedConfirmationRangeConfig.hpp" +#include "Acts/Seeding/TrigBase.hpp" //definition of Trigsispacepoint base and trigtriplets + +#include + +// core algorithm so in acts namespace +namespace Acts { + +template +class SeedFilter; + +template +struct SeedFinderFTFConfig { + // // how many sigmas of scattering angle should be considered? + float sigmaScattering = 5; + + // Seed cut + float minPt = 400. * Acts::UnitConstants::MeV; + + ///////////some declared not filled in by reco: ////// + std::shared_ptr> seedFilter; + + // //detector ROI + // // derived values, set on SeedFinder construction + float highland = 0; + float maxScatteringAngle2 = 0; + // bool isInInternalUnits = false; + /// for load space points + unsigned int maxSeedsPerSpM = 5; + + float m_phiSliceWidth; + float m_nMaxPhiSlice; + bool m_useClusterWidth = false; + std::string fastrack_input_file; + std::vector m_layerGeometry; + + // for run function + // m_settings: + bool m_LRTmode = true; // eventually want to set from full chain + bool m_useEtaBinning = true; + bool m_doubletFilterRZ = true; + float m_minDeltaRadius = 5.0; // eventually set in config or to equivalent + // acts 2.0 but increasing to test loops + // float m_maxDeltaRadius = 270.0 ; + float m_tripletD0Max = 4.0; // m_settings + unsigned int m_maxTripletBufferLength = 3; + + // ROI: + bool containsPhi() { + return false; + // need to implement this function + } + + //// + // 2 member functions + SeedFinderFTFConfig calculateDerivedQuantities() const { + // thorw statement if the isInternalUnits member is false, ie if dont call + // this function + SeedFinderFTFConfig config = *this; + // use a formula to calculate scattering + + return config; + } + + SeedFinderFTFConfig toInternalUnits() const { + // throw statement if the isInternalUnits member is false, ie if dont call + // this function + SeedFinderFTFConfig config = *this; + // divides inputs by 1mm, all ones input + // changes member inInInternalUnits to true + return config; + } + +}; // end of config struct + +} // namespace Acts diff --git a/Core/include/Acts/Seeding/TrigBase.hpp b/Core/include/Acts/Seeding/TrigBase.hpp new file mode 100644 index 00000000000..d407afe4736 --- /dev/null +++ b/Core/include/Acts/Seeding/TrigBase.hpp @@ -0,0 +1,45 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +#include "Acts/Seeding/GNN_TrackingFilter.hpp" + +#include + +#define MAX_SILICON_LAYER_NUM 19 +#define OffsetEndcapPixels 7 +#define OffsetBarrelSCT 3 +#define OffsetEndcapSCT 10 + +template +class TrigInDetTriplet { + public: + TrigInDetTriplet() = delete; // to prevent creation w/o initialization + + TrigInDetTriplet(Acts::FTF_SP s1, + Acts::FTF_SP s2, + Acts::FTF_SP s3, float Q) + : m_s1(std::move(s1)), m_s2(std::move(s2)), m_s3(std::move(s3)), m_Q(Q){}; + + TrigInDetTriplet(TrigInDetTriplet* t) + : m_s1(t->m_s1), m_s2(t->m_s2), m_s3(t->m_s3), m_Q(t->m_Q){}; + + const Acts::FTF_SP& s1() const { return m_s1; } + const Acts::FTF_SP& s2() const { return m_s2; } + const Acts::FTF_SP& s3() const { return m_s3; } + float Q() const { return m_Q; } + void Q(double newQ) { m_Q = newQ; } + + protected: + Acts::FTF_SP m_s1; + Acts::FTF_SP m_s2; + Acts::FTF_SP m_s3; + float m_Q; // Quality +}; diff --git a/Core/include/Acts/TrackFinding/FasTrackConnector.hpp b/Core/include/Acts/TrackFinding/FasTrackConnector.hpp new file mode 100644 index 00000000000..96b9e415651 --- /dev/null +++ b/Core/include/Acts/TrackFinding/FasTrackConnector.hpp @@ -0,0 +1,54 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +// Consider to moving to detail subdirectory +#include +#include +#include +#include + +namespace Acts { + +struct FasTrackConnection { + public: + FasTrackConnection(unsigned int, unsigned int); + ~FasTrackConnection(){}; + + unsigned int m_src, m_dst; + std::vector m_binTable; +}; + +class FasTrackConnector { + public: + struct LayerGroup { + LayerGroup(unsigned int l1Key, + const std::vector &v) + : m_dst(l1Key), m_sources(v){}; + + unsigned int m_dst; // the target layer of the group + std::vector + m_sources; // the source layers of the group + }; + + FasTrackConnector(std::ifstream &); + + ~FasTrackConnector(); + + float m_etaBin; + + std::map> m_layerGroups; + std::map> m_connMap; + // TODO: change to std::map > + // m_connMap; or std::map> > m_connMap; +}; + +} // namespace Acts diff --git a/Core/include/Acts/TrackFinding/RoiDescriptor.hpp b/Core/include/Acts/TrackFinding/RoiDescriptor.hpp new file mode 100644 index 00000000000..da5b5a26c5c --- /dev/null +++ b/Core/include/Acts/TrackFinding/RoiDescriptor.hpp @@ -0,0 +1,209 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +// TODO: update to C++17 style +#include +#include +#include +#include +#include + +#include + +namespace Acts { + +class RoiDescriptor { + public: + // iterator + typedef std::vector::const_iterator roi_iterator; + /// convenient + static constexpr bool FULLSCAN = true; + static constexpr bool ROI = false; + + /** + * @brief constructor + * @param eta eta of RoI + * @param etaMinus eta at rear of RoI + * @param etaPlus eta at front of RoI + * @param phi phi of RoI + * @param phiMinus minimum phi of RoI + * @param phiPlus maximum phi of RoI + * @param zed zed of RoI + * @param zedMinus zed at rear of RoI + * @param zedPlus zed at front of RoI + */ + RoiDescriptor(double eta, double etaMinus, double etaPlus, double phi, + double phiMinus, double phiPlus, double zed = 0, + double zedMinus = -s_zedWidthDefault, + double zedPlus = s_zedWidthDefault); + // zedminus - s_zedWidthDefault = 225 //from ROIDescriptor + + /* + * need an explicit class copy constructor + */ + RoiDescriptor(const RoiDescriptor& roi); + RoiDescriptor& operator=(const RoiDescriptor& r); + + // Destructor + ~RoiDescriptor(); + + // Methods to retrieve data members + + double phi() const { return m_phi; } + double eta() const { return m_eta; } + double zed() const { return m_zed; } + + /// these quantities probably don't need to be used any more + /// - they are implemented here only because we had them in + /// the original legacy interface + + double zedPlus() const { + return m_zedPlus; + } //!< z at the most forward end of the RoI + double zedMinus() const { + return m_zedMinus; + } //!< z at the most backward end of the RoI + + double etaPlus() const { return m_etaPlus; } //!< gets eta at zedPlus + double etaMinus() const { return m_etaMinus; } //!< gets eta at zMinus + + double phiPlus() const { return m_phiPlus; } //!< gets phiPlus + double phiMinus() const { return m_phiMinus; } //!< gets phiMinus + + /// versioning + int version() const { return m_version; } + void version(int v) { m_version = v; } + + /// output + // virtual operator std::string() const ; + + /// is this a full scan RoI? + bool isFullscan() const { return m_fullscan; } + + /// SuperRoI compatibility methods + + /// am I a SuperRoi? + bool composite() const { return m_composite; } + void setComposite(bool b = true) { m_composite = b; } + + /// always manage constituents ??? + bool manageConstituents() const { return m_manageConstituents; } + void manageConstituents(bool b) { m_manageConstituents = b; } + + /// number of constituents + unsigned size() const { return m_roiDescriptors.size(); } + + /// find an RoiDescriptor constituent + const RoiDescriptor* at(int i) const { return m_roiDescriptors.at(i); } + + /// clear the vector + void clear() { m_roiDescriptors.clear(); } // setComposite(false); } + + /// reserve elements in vector + void reserve(size_t s) { m_roiDescriptors.reserve(s); } + + /// add a RoiDescriptor + void push_back(const RoiDescriptor* roi) { + m_roiDescriptors.push_back(roi); + setComposite(true); + } + + /// iterators + roi_iterator begin() const { return m_roiDescriptors.begin(); } + roi_iterator end() const { return m_roiDescriptors.end(); } + + /// return the gradients + double dzdrMinus() const { + return m_dzdrMinus; + } //!< dz/dr at the rear of the RoI + double dzdrPlus() const { + return m_dzdrPlus; + } //!< dz/dr at the front of the RoI + + double drdzMinus() const { + return m_drdzMinus; + } //!< dr/dz at the rear of the RoI + double drdzPlus() const { + return m_drdzPlus; + } //!< dr/dz at the front of the RoI + + /// methods to calculate z position at the RoI boundary + /// at a given radius + double zedMin(double r) const; + double zedMax(double r) const; + + double zedOuterPlus() const { + return m_zedOuterPlus; + } //!< z at the most forward end of the RoI + double zedOuterMinus() const { + return m_zedOuterMinus; + } //!< z at the most backward end of the RoI + + double rhoMin(double z) const; + double rhoMax(double z) const; + + static double zedWidthDefault() { return s_zedWidthDefault; } + + /// set default z-width (but only before any RoiDescriptor has been created) + static void zedWidthDefault(double d); + + // fromn trig + unsigned int roiId() const { return m_roiId; } + unsigned int l1Id() const { return m_l1Id; } + unsigned int roiWord() const { return m_roiWord; } + + private: + /// default parameters - there may be better ways, but this will do + static std::atomic s_zedWidthDefault; + /// to ensure default width is only set once at job startup + static std::atomic s_firstInstanceCreated; + + float m_phi; //!< phi of RoI center + float m_eta; //!< eta of RoI center + float m_zed; //!< zed of RoI center + + float m_phiMinus; //!< most negative RoI in azimuthal + float m_phiPlus; //!< most positive RoI in azimuthal + float m_etaMinus; //!< eta of RoI at zedMinus + float m_etaPlus; //!< eta of RoI at zedPlus + float + m_zedMinus; //!< z position at most negative position along the beamline + float m_zedPlus; //!< z position at most positive position along the beamline + + float m_dzdrMinus; //!< dz/dr at the rear of the RoI + float m_dzdrPlus; //!< dz/dr at the front of the RoI + + float m_drdzMinus; //!< dr/dz at the rear of the RoI + float m_drdzPlus; //!< dr/dz at the front of the RoI + + float m_zedOuterMinus; //!< z at rear of RoI at the outer radius ( = 1100 mm) + float m_zedOuterPlus; //!< z at front of RoI at the outer radius ( = 1100 mm) + + bool m_fullscan; //!< flag this as a full detector RoI + bool m_composite; //!< flag this as a composite RoI + bool m_manageConstituents; //!< flag to determine whether constituents should + //!< be managed + + int m_version; //!< transient version identifier + + std::vector m_roiDescriptors; //!< roi constituents + + // from trig + unsigned int m_l1Id; //!< lvl1 event number + unsigned int m_roiId; //!< RoI number + unsigned int m_roiWord; //!< lvl1 RoI word from which this RoI was initially + //!< constructed + + // std::string str( const RoiDescriptor& d ); // +#include +#include +#include +#include + +namespace Acts { + +FasTrackConnection::FasTrackConnection(unsigned int s, unsigned int d) + : m_src(s), m_dst(d) {} + +FasTrackConnector::FasTrackConnector(std::ifstream &inFile) { + m_layerGroups.clear(); + + int nLinks; + + inFile >> nLinks >> m_etaBin; + + for (int l = 0; l < nLinks; l++) { + unsigned int stage, lIdx, src, dst, nEntries; + int height, width; + + inFile >> lIdx >> stage >> src >> dst >> height >> width >> nEntries; + + FasTrackConnection *pC = new FasTrackConnection(src, dst); + + int dummy; + + for (int i = 0; i < height; i++) { + for (int j = 0; j < width; j++) { + inFile >> dummy; // pC->m_binTable[j+i*width]; + } + } + + int vol_id = src / 1000; + + if (vol_id == 13 || vol_id == 12 || vol_id == 14) { + delete pC; + continue; + } + + vol_id = dst / 1000; + + if (vol_id == 13 || vol_id == 12 || vol_id == 14) { + delete pC; + continue; + } + + std::map>::iterator it = + m_connMap.find(stage); + + if (it == m_connMap.end()) { + std::vector v(1, pC); + m_connMap.insert(std::make_pair(stage, v)); + } else { + (*it).second.push_back(pC); + } + } + + // re-arrange the connection stages + + std::list lConns; + + std::map> newConnMap; + + for (const auto &conn : m_connMap) { + std::copy(conn.second.begin(), conn.second.end(), + std::back_inserter(lConns)); + } + + int stageCounter = 0; + + while (!lConns.empty()) { + std::unordered_map> + mCounter; // layerKey, nDst, nSrc + + for (const auto &conn : lConns) { + auto entryIt = mCounter.find(conn->m_dst); + if (entryIt != mCounter.end()) { + (*entryIt).second.first++; + } else { + int nDst = 1; + int nSrc = 0; + mCounter.insert( + std::make_pair(conn->m_dst, std::make_pair(nDst, nSrc))); + } + + entryIt = mCounter.find(conn->m_src); + if (entryIt != mCounter.end()) { + (*entryIt).second.second++; + } else { + int nDst = 0; + int nSrc = 1; + mCounter.insert( + std::make_pair(conn->m_src, std::make_pair(nDst, nSrc))); + } + } + + // find layers with nSrc = 0 + + std::set zeroLayers; + + for (const auto &layerCounts : mCounter) { + if (layerCounts.second.second != 0) { + continue; + } + + zeroLayers.insert(layerCounts.first); + } + + // remove connections which use zeroLayer as destination + + std::vector theStage; + + std::list::iterator cIt = lConns.begin(); + + while (cIt != lConns.end()) { + if (zeroLayers.find((*cIt)->m_dst) != + zeroLayers.end()) { // check if contains + theStage.push_back(*cIt); + cIt = lConns.erase(cIt); + continue; + } + cIt++; + } + newConnMap.insert(std::make_pair(stageCounter, theStage)); + stageCounter++; + } + + // create layer groups + + int currentStage = 0; + + // the doublet making is done using "outside-in" approach hence the reverse + // iterations + + for (std::map>::reverse_iterator + it = newConnMap.rbegin(); + it != newConnMap.rend(); ++it, currentStage++) { + const std::vector &vConn = (*it).second; + + // loop over links, extract all connections for the stage, group sources by + // L1 (dst) index + + std::map> l1ConnMap; + + for (const auto *conn : vConn) { + unsigned int dst = conn->m_dst; + + std::map>::iterator + l1MapIt = l1ConnMap.find(dst); + if (l1MapIt != l1ConnMap.end()) { + (*l1MapIt).second.push_back(conn); + } else { + std::vector v = {conn}; + l1ConnMap.insert(std::make_pair(dst, v)); + } + } + + std::vector lgv; + + lgv.reserve(l1ConnMap.size()); + + for (const auto &l1Group : l1ConnMap) { + lgv.push_back(LayerGroup(l1Group.first, l1Group.second)); + } + + m_layerGroups.insert(std::make_pair(currentStage, lgv)); + } + + newConnMap.clear(); +} + +FasTrackConnector::~FasTrackConnector() { + m_layerGroups.clear(); + for (std::map>::iterator it = + m_connMap.begin(); + it != m_connMap.end(); ++it) { + for (std::vector::iterator cIt = (*it).second.begin(); + cIt != (*it).second.end(); ++cIt) { + delete (*cIt); + } + } +} +} // namespace Acts diff --git a/Core/src/TrackFinding/RoiDescriptor.cpp b/Core/src/TrackFinding/RoiDescriptor.cpp new file mode 100644 index 00000000000..3ba80329fa7 --- /dev/null +++ b/Core/src/TrackFinding/RoiDescriptor.cpp @@ -0,0 +1,35 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// TODO: update to C++17 style +#include "Acts/TrackFinding/RoiDescriptor.hpp" + +#include +#include + +namespace Acts { + +Acts::RoiDescriptor::RoiDescriptor(double eta, double etaMinus, double etaPlus, + double phi, double phiMinus, double phiPlus, + double zed, double zedMinus, double zedPlus) + : m_phi(phi), + m_eta(eta), + m_zed(zed), + m_phiMinus(phiMinus), + m_phiPlus(phiPlus), + m_etaMinus(etaMinus), + m_etaPlus(etaPlus), + m_zedMinus(zedMinus), + m_zedPlus(zedPlus), + m_l1Id(0), + m_roiId(0), + m_roiWord(0) {} + +Acts::RoiDescriptor::~RoiDescriptor() {} + +} // namespace Acts diff --git a/Examples/Algorithms/TrackFinding/CMakeLists.txt b/Examples/Algorithms/TrackFinding/CMakeLists.txt index f55bcda945b..ab7b025c8a2 100644 --- a/Examples/Algorithms/TrackFinding/CMakeLists.txt +++ b/Examples/Algorithms/TrackFinding/CMakeLists.txt @@ -7,6 +7,7 @@ add_library( src/TrackFindingAlgorithmFunction.cpp src/HoughTransformSeeder.cpp src/TrackParamsEstimationAlgorithm.cpp + src/SeedingFTFAlgorithm.cpp ) target_include_directories( diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp new file mode 100644 index 00000000000..a1022d61566 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp @@ -0,0 +1,92 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// basing off of SeedingOrtho + +#pragma once + +#include "Acts/Seeding/InternalSeed.hpp" +#include "Acts/Seeding/SeedFilterConfig.hpp" +#include "Acts/Seeding/SeedFinderFTF.hpp" +#include "Acts/Seeding/SeedFinderFTFConfig.hpp" +// in core +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Seeding/SeedFinderConfig.hpp" +#include "Acts/Seeding/SpacePointGrid.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/EventData/SimSpacePoint.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +#include +#include + +namespace ActsExamples { + +class SeedingFTFAlgorithm final : public IAlgorithm { + public: + struct Config { + std::vector inputSpacePoints; + + std::string outputSeeds; + + Acts::SeedFilterConfig seedFilterConfig; + Acts::SeedFinderFTFConfig seedFinderConfig; + Acts::SeedFinderOptions seedFinderOptions; + + std::string layerMappingFile; + + std::vector geometrySelection; + + std::string inputSourceLinks; + + std::shared_ptr trackingGeometry; + + std::map, std::pair> ACTS_FTF_Map; + + bool fill_module_csv = false; + }; + + // constructor: + /// @param cfg is the algorithm configuration + /// @param lvl is the logging level + SeedingFTFAlgorithm(Config cfg, Acts::Logging::Level lvl); + + // code to make the algorithm run + /// @param txt is the algorithm context with event information + /// @return a process code indication success or failure + ProcessCode execute(const AlgorithmContext &ctx) const override; + + // access to config + const Config &config() const { return m_cfg; } + + // own class functions + // make the map + std::map, std::pair> Make_ACTS_FTF_Map() const; + // make the vector of space points with FTF Info + std::vector> Make_FTF_spacePoints( + const AlgorithmContext &ctx, + std::map, std::pair> map) const; + // layer numbering + std::vector LayerNumbering() const; + + private: + Config m_cfg; + + std::unique_ptr> mGNNgeo; + + std::vector>> + m_inputSpacePoints{}; + + WriteDataHandle m_outputSeeds{this, "OutputSeeds"}; + + ReadDataHandle m_inputSourceLinks{ + this, "InputSourceLinks"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp new file mode 100644 index 00000000000..76e8b3fb0a5 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/src/SeedingFTFAlgorithm.cpp @@ -0,0 +1,375 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp" + +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Seeding/Seed.hpp" +#include "Acts/Seeding/SeedFilter.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/ProtoTrack.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" + +#include +#include +#include +#include +#include +#include + +using namespace std; + +template class Acts::TrigFTF_GNN_Layer; +template class Acts::TrigFTF_GNN_Geometry; +template class Acts::TrigFTF_GNN_Node; +template class Acts::TrigFTF_GNN_EtaBin; +template struct Acts::FTF_SP; +template class Acts::TrigFTF_GNN_DataStorage; +template class Acts::TrigFTF_GNN_Edge; + +// constructor: +ActsExamples::SeedingFTFAlgorithm::SeedingFTFAlgorithm( + ActsExamples::SeedingFTFAlgorithm::Config cfg, Acts::Logging::Level lvl) + : ActsExamples::IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) { + // fill config struct + m_cfg.layerMappingFile = m_cfg.layerMappingFile; + + m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits(); + + m_cfg.seedFinderConfig = + m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities(); + + m_cfg.seedFinderOptions = + m_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities( + m_cfg.seedFinderConfig); + + for (const auto &spName : m_cfg.inputSpacePoints) { + if (spName.empty()) { + throw std::invalid_argument("Invalid space point input collection"); + } + + auto &handle = m_inputSpacePoints.emplace_back( + std::make_unique>( + this, + "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size()))); + handle->initialize(spName); + } + + m_outputSeeds.initialize(m_cfg.outputSeeds); + + m_inputSourceLinks.initialize(m_cfg.inputSourceLinks); + + m_cfg.seedFinderConfig.seedFilter = + std::make_unique>( + Acts::SeedFilter(m_cfg.seedFilterConfig)); + + // map + m_cfg.ACTS_FTF_Map = Make_ACTS_FTF_Map(); + // input trig vector + m_cfg.seedFinderConfig.m_layerGeometry = LayerNumbering(); + + std::ifstream input_ifstream( + m_cfg.seedFinderConfig.fastrack_input_file.c_str(), std::ifstream::in); + + // fastrack + std::unique_ptr input_fastrack = + std::make_unique(input_ifstream); + + mGNNgeo = std::make_unique>( + m_cfg.seedFinderConfig.m_layerGeometry, input_fastrack); + +} // this is not FTF config type because it is a member of the algs config, + // which is of type FTF cofig + +// execute: +ActsExamples::ProcessCode ActsExamples::SeedingFTFAlgorithm::execute( + const AlgorithmContext &ctx) const { + std::vector> FTF_spacePoints = + Make_FTF_spacePoints(ctx, m_cfg.ACTS_FTF_Map); + + for (auto sp : FTF_spacePoints) { + ACTS_DEBUG("FTF space points: " + << " FTF_id: " << sp.FTF_ID << " z: " << sp.SP->z() + << " r: " << sp.SP->r() << " ACTS volume: " + << sp.SP->sourceLinks() + .front() + .get() + .geometryId() + .volume() + << "\n"); + } + + // this is now calling on a core algorithm + Acts::SeedFinderFTF finder(m_cfg.seedFinderConfig, *mGNNgeo); + + // need this function as create_coords is needed for seeds + std::function( + const SimSpacePoint *sp)> + create_coordinates = [](const SimSpacePoint *sp) { + Acts::Vector3 position(sp->x(), sp->y(), sp->z()); + Acts::Vector2 variance(sp->varianceR(), sp->varianceZ()); + return std::make_pair(position, variance); + }; + // output of function needed for seed + + finder.loadSpacePoints(FTF_spacePoints); + + Acts::RoiDescriptor internalRoi(0, -4.5, 4.5, 0, -M_PI, M_PI, 0, -150.0, + 150.0); + + finder.createSeeds(internalRoi, *mGNNgeo); + + // still to develop + SimSeedContainer seeds = finder.createSeeds_old( + m_cfg.seedFinderOptions, FTF_spacePoints, create_coordinates); + + m_outputSeeds(ctx, std::move(seeds)); + + return ActsExamples::ProcessCode::SUCCESS; +} + +std::map, std::pair> +ActsExamples::SeedingFTFAlgorithm::Make_ACTS_FTF_Map() const { + map, std::pair> ACTS_FTF; + std::ifstream data( + m_cfg.layerMappingFile); // 0 in this file refers to no FTF ID + std::string line; + std::vector> parsedCsv; + while (std::getline(data, line)) { + std::stringstream lineStream(line); + std::string cell; + std::vector parsedRow; + while (std::getline(lineStream, cell, ',')) { + parsedRow.push_back(cell); + } + + parsedCsv.push_back(parsedRow); + } + // file in format ACTS_vol,ACTS_lay,ACTS_mod,FTF_id + for (auto i : parsedCsv) { + int ACTS_vol = stoi(i[0]); + int ACTS_lay = stoi(i[1]); + int ACTS_mod = stoi(i[2]); + int FTF = stoi(i[5]); + int eta_mod = stoi(i[6]); + int ACTS_joint = ACTS_vol * 100 + ACTS_lay; + ACTS_FTF.insert({{ACTS_joint, ACTS_mod}, {FTF, eta_mod}}); + } + + return ACTS_FTF; +} + +std::vector> +ActsExamples::SeedingFTFAlgorithm::Make_FTF_spacePoints( + const AlgorithmContext &ctx, + std::map, std::pair> map) const { + // create space point vectors + std::vector spacePoints; + std::vector> FTF_spacePoints; + FTF_spacePoints.reserve( + m_inputSpacePoints.size()); // not sure if this is enough + + // for loop filling space + for (const auto &isp : m_inputSpacePoints) { + for (const auto &spacePoint : (*isp)(ctx)) { + // fill original space point vector + spacePoints.push_back(&spacePoint); + + // FTF space point vector + // loop over space points, call on map + const auto &source_link = spacePoint.sourceLinks(); + const auto &index_source_link = + source_link.front().get(); + + // warning if source link empty + if (source_link.empty()) { + // warning in officaial acts format + ACTS_WARNING("warning source link vector is empty"); + continue; + } + int ACTS_vol_id = index_source_link.geometryId().volume(); + int ACTS_lay_id = index_source_link.geometryId().layer(); + int ACTS_mod_id = index_source_link.geometryId().sensitive(); + + // dont want strips or HGTD + if (ACTS_vol_id == 2 or ACTS_vol_id == 22 or ACTS_vol_id == 23 or + ACTS_vol_id == 24) { + continue; + } + + // Search for vol, lay and module=0, if doesn't esist (end) then search + // for full thing vol*100+lay as first number in pair then 0 or mod id + auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id; + auto key = std::make_pair( + ACTS_joint_id, + 0); // here the key needs to be pair of(vol*100+lay, 0) + auto Find = map.find(key); + + if (Find == + map.end()) { // if end then make new key of (vol*100+lay, modid) + key = std::make_pair(ACTS_joint_id, ACTS_mod_id); // mod ID + Find = map.find(key); + } + + // warning if key not in map + if (Find == map.end()) { + ACTS_WARNING("Key not found in FTF map for volume id: " + << ACTS_vol_id << " and layer id: " << ACTS_lay_id); + continue; + } + + // now should be pixel with FTF ID: + int FTF_id = + Find->second + .first; // new map the item is a pair so want first from it + + if (FTF_id == 0) { + ACTS_WARNING("No assigned FTF ID for key for volume id: " + << ACTS_vol_id << " and layer id: " << ACTS_lay_id); + } + + // access IDs from map + int eta_mod = Find->second.second; + int combined_id = FTF_id * 1000 + eta_mod; + + // fill FTF vector with current sapce point and ID + FTF_spacePoints.emplace_back(&spacePoint, FTF_id, combined_id); + } + } + ACTS_VERBOSE("Space points successfully assigned FTF ID"); + + return FTF_spacePoints; +} + +std::vector +ActsExamples::SeedingFTFAlgorithm::LayerNumbering() const { + std::vector input_vector; + std::vector count_vector; + + m_cfg.trackingGeometry->visitSurfaces([this, &input_vector, &count_vector]( + const Acts::Surface *surface) { + Acts::GeometryIdentifier geoid = surface->geometryId(); + auto ACTS_vol_id = geoid.volume(); + auto ACTS_lay_id = geoid.layer(); + auto mod_id = geoid.sensitive(); + auto bounds_vect = surface->bounds().values(); + auto center = surface->center(Acts::GeometryContext()); + + // make bounds global + Acts::Vector3 globalFakeMom(1, 1, 1); + Acts::Vector2 min_bound_local = + Acts::Vector2(bounds_vect[0], bounds_vect[1]); + Acts::Vector2 max_bound_local = + Acts::Vector2(bounds_vect[2], bounds_vect[3]); + Acts::Vector3 min_bound_global = surface->localToGlobal( + Acts::GeometryContext(), min_bound_local, globalFakeMom); + Acts::Vector3 max_bound_global = surface->localToGlobal( + Acts::GeometryContext(), max_bound_local, globalFakeMom); + + if (min_bound_global(0) > + max_bound_global(0)) { // checking that not wrong way round + min_bound_global.swap(max_bound_global); + } + + float rc = 0.0; + float minBound = 100000.0; + float maxBound = -100000.0; + + // convert to FTF ID + auto ACTS_joint_id = ACTS_vol_id * 100 + ACTS_lay_id; + auto key = + std::make_pair(ACTS_joint_id, + 0); // here the key needs to be pair of(vol*100+lay, 0) + auto Find = m_cfg.ACTS_FTF_Map.find(key); + int FTF_id = Find->second.first; // new map, item is pair want first + if (Find == + m_cfg.ACTS_FTF_Map + .end()) { // if end then make new key of (vol*100+lay, modid) + key = std::make_pair(ACTS_joint_id, mod_id); // mod ID + Find = m_cfg.ACTS_FTF_Map.find(key); + FTF_id = Find->second.first; + } + + short barrel_ec = 0; // a variable that says if barrrel, 0 = barrel + int eta_mod = Find->second.second; + + // assign barrel_ec depending on FTF_layer + if (79 < FTF_id && FTF_id < 85) { // 80s, barrel + barrel_ec = 0; + } else if (89 < FTF_id && FTF_id < 99) { // 90s positive + barrel_ec = 2; + } else { // 70s negative + barrel_ec = -2; + } + + if (barrel_ec == 0) { + rc = sqrt(center(0) * center(0) + + center(1) * center(1)); // barrel center in r + // bounds of z + if (min_bound_global(2) < minBound) + minBound = min_bound_global(2); + if (max_bound_global(2) > maxBound) + maxBound = max_bound_global(2); + } else { + rc = center(2); // not barrel center in Z + // bounds of r + float min = sqrt(min_bound_global(0) * min_bound_global(0) + + min_bound_global(1) * min_bound_global(1)); + float max = sqrt(max_bound_global(0) * max_bound_global(0) + + max_bound_global(1) * max_bound_global(1)); + if (min < minBound) + minBound = min; + if (max > maxBound) + maxBound = max; + } + + int combined_id = FTF_id * 1000 + eta_mod; + auto current_index = + find_if(input_vector.begin(), input_vector.end(), + [combined_id](auto n) { return n.m_subdet == combined_id; }); + if (current_index != input_vector.end()) { // not end so does exist + size_t index = std::distance(input_vector.begin(), current_index); + input_vector[index].m_refCoord += rc; + input_vector[index].m_minBound += minBound; + input_vector[index].m_maxBound += maxBound; + count_vector[index] += 1; // increase count at the index + + } else { // end so doesn't exists + // make new if one with FTF ID doesn't exist: + Acts::TrigInDetSiLayer new_FTF_ID(combined_id, barrel_ec, rc, minBound, + maxBound); + input_vector.push_back(new_FTF_ID); + count_vector.push_back( + 1); // so the element exists and not divinding by 0 + } + + if (m_cfg.fill_module_csv) { + fstream fout; + fout.open("ACTS_modules.csv", + ios::out | ios::app); // add to file each time + // print to csv for each module, no repeats so dont need to make + // map for averaging + fout << ACTS_vol_id << ", " // vol + << ACTS_lay_id << ", " // lay + << mod_id << ", " // module + << FTF_id << "," // FTF id + << eta_mod << "," // eta_mod + << center(2) << ", " // z + << sqrt(center(0) * center(0) + center(1) * center(1)) // r + << "\n"; + } + }); + + for (long unsigned int i = 0; i < input_vector.size(); i++) { + input_vector[i].m_refCoord = input_vector[i].m_refCoord / count_vector[i]; + } + + return input_vector; +} diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 55a5c6e657c..42dfa4f4ab6 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -10,7 +10,8 @@ u = acts.UnitConstants SeedingAlgorithm = Enum( - "SeedingAlgorithm", "Default TruthSmeared TruthEstimated Orthogonal HoughTransform" + "SeedingAlgorithm", + "Default TruthSmeared TruthEstimated Orthogonal HoughTransform FTF", ) TruthSeedRanges = namedtuple( @@ -167,6 +168,8 @@ def addSeeding( trackingGeometry: acts.TrackingGeometry, field: acts.MagneticFieldProvider, geoSelectionConfigFile: Optional[Union[Path, str]] = None, + layerMappingConfigFile: Optional[Union[Path, str]] = None, + fastrack_inputConfigFile: Optional[Union[Path, str]] = None, seedingAlgorithm: SeedingAlgorithm = SeedingAlgorithm.Default, truthSeedRanges: Optional[TruthSeedRanges] = TruthSeedRanges(), particleSmearingSigmas: ParticleSmearingSigmas = ParticleSmearingSigmas(), @@ -314,6 +317,21 @@ def addSeeding( houghTransformConfig.outputSeeds = "seeds" houghTransformConfig.trackingGeometry = trackingGeometry seeds = addHoughTransformSeeding(s, houghTransformConfig, logLevel) + elif seedingAlgorithm == SeedingAlgorithm.FTF: + logger.info("Using FTF seeding") + # output of algs changed, only one output now + seeds = addFTFSeeding( + s, + spacePoints, + seedFinderConfigArg, + seedFinderOptionsArg, + seedFilterConfigArg, + trackingGeometry, + logLevel, + layerMappingConfigFile, + geoSelectionConfigFile, + fastrack_inputConfigFile, + ) else: logger.fatal("unknown seedingAlgorithm %s", seedingAlgorithm) @@ -778,6 +796,84 @@ def addHoughTransformSeeding( return ht.config.outputSeeds +def addFTFSeeding( + sequence: acts.examples.Sequencer, + spacePoints: str, + seedFinderConfigArg: SeedFinderConfigArg, + seedFinderOptionsArg: SeedFinderOptionsArg, + seedFilterConfigArg: SeedFilterConfigArg, + trackingGeometry: acts.TrackingGeometry, + logLevel: acts.logging.Level = None, + layerMappingConfigFile: Union[Path, str] = None, + geoSelectionConfigFile: Union[Path, str] = None, + fastrack_inputConfigFile: Union[Path, str] = None, +): + """FTF seeding""" + + logLevel = acts.examples.defaultLogging(sequence, logLevel)() + layerMappingFile = str(layerMappingConfigFile) # turn path into string + fastrack_inputFile = str(fastrack_inputConfigFile) + seedFinderConfig = acts.SeedFinderFTFConfig( + **acts.examples.defaultKWArgs( + sigmaScattering=seedFinderConfigArg.sigmaScattering, + maxSeedsPerSpM=seedFinderConfigArg.maxSeedsPerSpM, + minPt=seedFinderConfigArg.minPt, + fastrack_input_file=fastrack_inputFile, + m_useClusterWidth=False, + ), + ) + seedFinderOptions = acts.SeedFinderOptions( + **acts.examples.defaultKWArgs( + beamPos=acts.Vector2(0.0, 0.0) + if seedFinderOptionsArg.beamPos == (None, None) + else acts.Vector2( + seedFinderOptionsArg.beamPos[0], seedFinderOptionsArg.beamPos[1] + ), + bFieldInZ=seedFinderOptionsArg.bFieldInZ, + ) + ) + seedFilterConfig = acts.SeedFilterConfig( + **acts.examples.defaultKWArgs( + maxSeedsPerSpM=seedFinderConfig.maxSeedsPerSpM, + deltaRMin=( + seedFinderConfigArg.deltaR[0] + if seedFilterConfigArg.deltaRMin is None + else seedFilterConfigArg.deltaRMin + ), + impactWeightFactor=seedFilterConfigArg.impactWeightFactor, + zOriginWeightFactor=seedFilterConfigArg.zOriginWeightFactor, + compatSeedWeight=seedFilterConfigArg.compatSeedWeight, + compatSeedLimit=seedFilterConfigArg.compatSeedLimit, + numSeedIncrement=seedFilterConfigArg.numSeedIncrement, + seedWeightIncrement=seedFilterConfigArg.seedWeightIncrement, + seedConfirmation=seedFilterConfigArg.seedConfirmation, + # curvatureSortingInFilter=seedFilterConfigArg.curvatureSortingInFilter, + maxSeedsPerSpMConf=seedFilterConfigArg.maxSeedsPerSpMConf, + maxQualitySeedsPerSpMConf=seedFilterConfigArg.maxQualitySeedsPerSpMConf, + useDeltaRorTopRadius=seedFilterConfigArg.useDeltaRorTopRadius, + ) + ) + + seedingAlg = acts.examples.SeedingFTFAlgorithm( + level=logLevel, + inputSpacePoints=[spacePoints], + outputSeeds="seeds", + seedFilterConfig=seedFilterConfig, + seedFinderConfig=seedFinderConfig, + seedFinderOptions=seedFinderOptions, + layerMappingFile=layerMappingFile, + geometrySelection=acts.examples.readJsonGeometryList( + str(geoSelectionConfigFile) + ), + inputSourceLinks="sourcelinks", + trackingGeometry=trackingGeometry, + fill_module_csv=False, + ) + + sequence.addAlgorithm(seedingAlg) + return seedingAlg.config.outputSeeds + + def addSeedPerformanceWriters( sequence: acts.examples.Sequencer, outputDirRoot: Union[Path, str], diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index 68f941734b5..66eab7d4157 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -12,6 +12,7 @@ #include "Acts/Seeding/SeedConfirmationRangeConfig.hpp" #include "Acts/Seeding/SeedFilterConfig.hpp" #include "Acts/Seeding/SeedFinderConfig.hpp" +#include "Acts/Seeding/SeedFinderFTFConfig.hpp" #include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp" #include "Acts/Seeding/SpacePointGrid.hpp" #include "Acts/TrackFinding/MeasurementSelector.hpp" @@ -20,6 +21,7 @@ #include "ActsExamples/EventData/Track.hpp" #include "ActsExamples/TrackFinding/HoughTransformSeeder.hpp" #include "ActsExamples/TrackFinding/SeedingAlgorithm.hpp" +#include "ActsExamples/TrackFinding/SeedingFTFAlgorithm.hpp" #include "ActsExamples/TrackFinding/SeedingOrthogonalAlgorithm.hpp" #include "ActsExamples/TrackFinding/SpacePointMaker.hpp" #include "ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp" @@ -192,6 +194,24 @@ void addTrackFinding(Context& ctx) { patchKwargsConstructor(c); } + { + using Config = Acts::SeedFinderFTFConfig; + auto c = py::class_(m, "SeedFinderFTFConfig").def(py::init<>()); + ACTS_PYTHON_STRUCT_BEGIN(c, Config); + ACTS_PYTHON_MEMBER(minPt); + ACTS_PYTHON_MEMBER(sigmaScattering); + ACTS_PYTHON_MEMBER(highland); + ACTS_PYTHON_MEMBER(maxScatteringAngle2); + ACTS_PYTHON_MEMBER(fastrack_input_file); + ACTS_PYTHON_MEMBER(m_phiSliceWidth); + ACTS_PYTHON_MEMBER(m_nMaxPhiSlice); + ACTS_PYTHON_MEMBER(m_useClusterWidth); + ACTS_PYTHON_MEMBER(m_layerGeometry); + ACTS_PYTHON_MEMBER(maxSeedsPerSpM); + ACTS_PYTHON_STRUCT_END(); + patchKwargsConstructor(c); + } + { using seedConf = Acts::SeedConfirmationRangeConfig; auto c = py::class_(m, "SeedConfirmationRangeConfig") @@ -250,6 +270,12 @@ void addTrackFinding(Context& ctx) { outputSeeds, seedFilterConfig, seedFinderConfig, seedFinderOptions); + ACTS_PYTHON_DECLARE_ALGORITHM( + ActsExamples::SeedingFTFAlgorithm, mex, "SeedingFTFAlgorithm", + inputSpacePoints, outputSeeds, seedFilterConfig, seedFinderConfig, + seedFinderOptions, layerMappingFile, geometrySelection, inputSourceLinks, + trackingGeometry, ACTS_FTF_Map, fill_module_csv); + ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::HoughTransformSeeder, mex, "HoughTransformSeeder", inputSpacePoints, outputProtoTracks, inputSourceLinks, trackingGeometry, diff --git a/Examples/Scripts/Python/full_chain_itk_FTF.py b/Examples/Scripts/Python/full_chain_itk_FTF.py new file mode 100755 index 00000000000..83f98a666a6 --- /dev/null +++ b/Examples/Scripts/Python/full_chain_itk_FTF.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +import pathlib, acts, acts.examples, acts.examples.itk +from acts.examples.simulation import ( + addParticleGun, + MomentumConfig, + EtaConfig, + ParticleConfig, + addPythia8, + addFatras, + ParticleSelectorConfig, + addDigitization, +) +from acts.examples.reconstruction import ( + addSeeding, + SeedingAlgorithm, + TruthSeedRanges, + addCKFTracks, + TrackSelectorConfig, + addAmbiguityResolution, + AmbiguityResolutionConfig, + addVertexFitting, + VertexFinder, +) + +ttbar_pu200 = False +u = acts.UnitConstants +geo_dir = pathlib.Path("acts-itk") +outputDir = pathlib.Path.cwd() / "itk_output" +# acts.examples.dump_args_calls(locals()) # show acts.examples python binding calls + +detector, trackingGeometry, decorators = acts.examples.itk.buildITkGeometry(geo_dir) +field = acts.examples.MagneticFieldMapXyz(str(geo_dir / "bfield/ATLAS-BField-xyz.root")) +rnd = acts.examples.RandomNumbers(seed=42) + +s = acts.examples.Sequencer(events=100, numThreads=1, outputDir=str(outputDir)) + +if not ttbar_pu200: + addParticleGun( + s, + MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), + EtaConfig(-4.0, 4.0, uniform=True), + ParticleConfig(2, acts.PdgParticle.eMuon, randomizeCharge=True), + rnd=rnd, + ) +else: + addPythia8( + s, + hardProcess=["Top:qqbar2ttbar=on"], + npileup=200, + vtxGen=acts.examples.GaussianVertexGenerator( + stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 5.0 * u.ns), + mean=acts.Vector4(0, 0, 0, 0), + ), + rnd=rnd, + outputDirRoot=outputDir, + ) + +addFatras( + s, + trackingGeometry, + field, + rnd=rnd, + preSelectParticles=ParticleSelectorConfig( + rho=(0.0 * u.mm, 28.0 * u.mm), + absZ=(0.0 * u.mm, 1.0 * u.m), + eta=(-4.0, 4.0), + pt=(150 * u.MeV, None), + removeNeutral=True, + ) + if ttbar_pu200 + else ParticleSelectorConfig(), + outputDirRoot=outputDir, +) + +addDigitization( + s, + trackingGeometry, + field, + digiConfigFile=geo_dir + / "itk-hgtd/itk-smearing-config.json", # change this file to make it do digitization + outputDirRoot=outputDir, + rnd=rnd, +) + +addSeeding( + s, + trackingGeometry, + field, + TruthSeedRanges(pt=(1.0 * u.GeV, None), eta=(-4.0, 4.0), nHits=(9, None)) + if ttbar_pu200 + else TruthSeedRanges(), + seedingAlgorithm=SeedingAlgorithm.FTF, + *acts.examples.itk.itkSeedingAlgConfig( + acts.examples.itk.InputSpacePointsType.PixelSpacePoints + ), + geoSelectionConfigFile=geo_dir / "itk-hgtd/geoSelection-ITk.json", + layerMappingConfigFile=geo_dir / "itk-hgtd/ACTS_FTF_mapinput.csv", + fastrack_inputConfigFile=geo_dir / "itk-hgtd/binTables_ITK_RUN4.txt", + outputDirRoot=outputDir, +) + +s.run()