Acts/Seeding/GbtsTrackingFilter.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // TODO: update to C++17 style
0012 #include "Acts/Seeding/GbtsDataStorage.hpp"  //includes geo which has trigindetsilayer, may move this to trigbase
0013 #include "Acts/Utilities/Logger.hpp"
0014
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <cstring>
0018 #include <iostream>
0019 #include <list>
0020 #include <vector>
0021
0022 namespace Acts::Experimental {
0023
0024 template <typename external_spacepoint_t>
0025 struct GbtsEdgeState {
0026  public:
0027   struct Compare {
0028     bool operator()(const struct GbtsEdgeState* s1,
0029                     const struct GbtsEdgeState* s2) {
0030       return s1->m_J > s2->m_J;
0031     }
0032   };
0033
0034   GbtsEdgeState() = default;
0035
0036   GbtsEdgeState(bool f) : m_initialized(f) {}
0037
0038   void initialize(GbtsEdge<external_spacepoint_t>* pS) {
0039     m_initialized = true;
0040
0041     m_J = 0.0;
0042     m_vs.clear();
0043
0044     // n2->n1
0045
0046     float dx = pS->m_n1->m_spGbts.SP->x() - pS->m_n2->m_spGbts.SP->x();
0047     float dy = pS->m_n1->m_spGbts.SP->y() - pS->m_n2->m_spGbts.SP->y();
0048     float L = std::sqrt(dx * dx + dy * dy);
0049
0050     m_s = dy / L;
0051     m_c = dx / L;
0052
0053     // transform for extrapolation and update
0054     //  x' =  x*m_c + y*m_s
0055     //  y' = -x*m_s + y*m_c
0056
0057     m_refY = pS->m_n2->m_spGbts.r();
0058     m_refX =
0059         pS->m_n2->m_spGbts.SP->x() * m_c + pS->m_n2->m_spGbts.SP->y() * m_s;
0060
0061     // X-state: y, dy/dx, d2y/dx2
0062
0063     m_X[0] =
0064         -pS->m_n2->m_spGbts.SP->x() * m_s + pS->m_n2->m_spGbts.SP->y() * m_c;
0065     m_X[1] = 0.0;
0066     m_X[2] = 0.0;
0067
0068     // Y-state: z, dz/dr
0069
0070     m_Y[0] = pS->m_n2->m_spGbts.SP->z();
0071     m_Y[1] = (pS->m_n1->m_spGbts.SP->z() - pS->m_n2->m_spGbts.SP->z()) /
0072              (pS->m_n1->m_spGbts.r() - pS->m_n2->m_spGbts.r());
0073
0074     memset(&m_Cx[0][0], 0, sizeof(m_Cx));
0075     memset(&m_Cy[0][0], 0, sizeof(m_Cy));
0076
0077     m_Cx[0][0] = 0.25;
0078     m_Cx[1][1] = 0.001;
0079     m_Cx[2][2] = 0.001;
0080
0081     m_Cy[0][0] = 1.5;
0082     m_Cy[1][1] = 0.001;
0083   }
0084
0085   void clone(const struct GbtsEdgeState& st) {
0086     memcpy(&m_X[0], &st.m_X[0], sizeof(m_X));
0087     memcpy(&m_Y[0], &st.m_Y[0], sizeof(m_Y));
0088     memcpy(&m_Cx[0][0], &st.m_Cx[0][0], sizeof(m_Cx));
0089     memcpy(&m_Cy[0][0], &st.m_Cy[0][0], sizeof(m_Cy));
0090     m_refX = st.m_refX;
0091     m_refY = st.m_refY;
0092     m_c = st.m_c;
0093     m_s = st.m_s;
0094     m_J = st.m_J;
0095     m_vs.clear();
0096     m_vs.reserve(st.m_vs.size());
0097     std::copy(st.m_vs.begin(), st.m_vs.end(), std::back_inserter(m_vs));
0098
0099     m_initialized = true;
0100   }
0101
0102   float m_J{};
0103
0104   std::vector<GbtsEdge<external_spacepoint_t>*> m_vs;
0105
0106   float m_X[3]{}, m_Y[2]{}, m_Cx[3][3]{}, m_Cy[2][2]{};
0107   float m_refX{}, m_refY{}, m_c{}, m_s{};
0108
0109   bool m_initialized{false};
0110 };
0111
0112 #define MAX_EDGE_STATE 2500
0113
0114 template <typename external_spacepoint_t>
0115 class GbtsTrackingFilter {
0116  public:
0117   GbtsTrackingFilter(const std::vector<TrigInDetSiLayer>& g,
0118                      std::vector<GbtsEdge<external_spacepoint_t>>& sb,
0119                      std::unique_ptr<const Acts::Logger> logger =
0120                          Acts::getDefaultLogger("Filter",
0121                                                 Acts::Logging::Level::INFO))
0122       : m_geo(g), m_segStore(sb), m_logger(std::move(logger)) {}
0123
0124   void followTrack(GbtsEdge<external_spacepoint_t>* pS,
0125                    GbtsEdgeState<external_spacepoint_t>& output) {
0126     if (pS->m_level == -1) {
0127       return;  // already collected
0128     }
0129     m_globalStateCounter = 0;
0130
0131     // create track state
0132
0133     GbtsEdgeState<external_spacepoint_t>* pInitState =
0134         &m_stateStore[m_globalStateCounter++];
0135
0136     pInitState->initialize(pS);
0137
0138     m_stateVec.clear();
0139
0140     // recursive branching and propagation
0141
0142     propagate(pS, *pInitState);
0143
0144     if (m_stateVec.empty()) {
0145       return;
0146     }
0147     std::ranges::sort(m_stateVec,
0148                       typename GbtsEdgeState<external_spacepoint_t>::Compare());
0149
0150     GbtsEdgeState<external_spacepoint_t>* best = (*m_stateVec.begin());
0151
0152     output.clone(*best);
0153
0154     m_globalStateCounter = 0;
0155   }
0156
0157  protected:
0158   void propagate(GbtsEdge<external_spacepoint_t>* pS,
0159                  GbtsEdgeState<external_spacepoint_t>& ts) {
0160     if (m_globalStateCounter >= MAX_EDGE_STATE) {
0161       return;
0162     }
0163     GbtsEdgeState<external_spacepoint_t>* p_new_ts =
0164         &m_stateStore[m_globalStateCounter++];
0165
0166     GbtsEdgeState<external_spacepoint_t>& new_ts = *p_new_ts;
0167     new_ts.clone(ts);
0168
0169     new_ts.m_vs.push_back(pS);
0170
0171     bool accepted = update(pS, new_ts);  // update using n1 of the segment
0172
0173     if (!accepted) {
0174       return;  // stop further propagation
0175     }
0176     int level = pS->m_level;
0177
0178     std::list<GbtsEdge<external_spacepoint_t>*> lCont;
0179
0180     for (int nIdx = 0; nIdx < pS->m_nNei;
0181          nIdx++) {  // loop over the neighbours of this segment
0182       unsigned int nextSegmentIdx = pS->m_vNei[nIdx];
0183
0184       GbtsEdge<external_spacepoint_t>* pN = &(m_segStore.at(nextSegmentIdx));
0185
0186       if (pN->m_level == -1) {
0187         continue;  // already collected
0188       }
0189       if (pN->m_level == level - 1) {
0190         lCont.push_back(pN);
0191       }
0192     }
0193     if (lCont.empty()) {  // the end of chain
0194
0195       // store in the vector
0196       if (m_globalStateCounter < MAX_EDGE_STATE) {
0197         if (m_stateVec.empty()) {  // add the first segment state
0198           GbtsEdgeState<external_spacepoint_t>* p =
0199               &m_stateStore[m_globalStateCounter++];
0200           p->clone(new_ts);
0201           m_stateVec.push_back(p);
0202         } else {  // compare with the best and add
0203           float best_so_far = (*m_stateVec.begin())->m_J;
0204           if (new_ts.m_J > best_so_far) {
0205             GbtsEdgeState<external_spacepoint_t>* p =
0206                 &m_stateStore[m_globalStateCounter++];
0207             p->clone(new_ts);
0208             m_stateVec.push_back(p);
0209           }
0210         }
0211       }
0212     } else {  // branching
0213       int nBranches = 0;
0214       for (typename std::list<GbtsEdge<external_spacepoint_t>*>::iterator sIt =
0215                lCont.begin();
0216            sIt != lCont.end(); ++sIt, nBranches++) {
0217         propagate((*sIt), new_ts);  // recursive call
0218       }
0219     }
0220   }
0221
0222   bool update(GbtsEdge<external_spacepoint_t>* pS,
0223               GbtsEdgeState<external_spacepoint_t>& ts) {
0224     const float sigma_t = 0.0003;
0225     const float sigma_w = 0.00009;
0226
0227     const float sigmaMS = 0.016;
0228
0229     const float sigma_x = 0.25;  // was 0.22
0230     const float sigma_y = 2.5;   // was 1.7
0231
0232     const float weight_x = 0.5;
0233     const float weight_y = 0.5;
0234
0235     const float maxDChi2_x = 60.0;  // was 35.0;
0236     const float maxDChi2_y = 60.0;  // was 31.0;
0237
0238     const float add_hit = 14.0;
0239
0240     if (ts.m_Cx[2][2] < 0.0 || ts.m_Cx[1][1] < 0.0 || ts.m_Cx[0][0] < 0.0) {
0241       ACTS_WARNING("Negative covariance detected in X components: "
0242                    << "cov[2][2]=" << ts.m_Cx[2][2] << " cov[1][1]="
0243                    << ts.m_Cx[1][1] << " cov[0][0]=" << ts.m_Cx[0][0]);
0244     }
0245
0246     if (ts.m_Cy[1][1] < 0.0 || ts.m_Cy[0][0] < 0.0) {
0247       ACTS_WARNING("Negative covariance detected in Y components: "
0248                    << "cov[1][1]=" << ts.m_Cy[1][1]
0249                    << " cov[0][0]=" << ts.m_Cy[0][0]);
0250     }
0251
0252     // add ms.
0253
0254     ts.m_Cx[2][2] += sigma_w * sigma_w;
0255     ts.m_Cx[1][1] += sigma_t * sigma_t;
0256
0257     int type1 = getLayerType(pS->m_n1->m_spGbts.combined_ID);
0258
0259     float t2 = type1 == 0 ? 1.0 + ts.m_Y[1] * ts.m_Y[1]
0260                           : 1.0 + 1.0 / (ts.m_Y[1] * ts.m_Y[1]);
0261     float s1 = sigmaMS * t2;
0262     float s2 = s1 * s1;
0263
0264     s2 *= std::sqrt(t2);
0265
0266     ts.m_Cy[1][1] += s2;
0267
0268     // extrapolation
0269
0270     float X[3], Y[2];
0271     float Cx[3][3], Cy[2][2];
0272
0273     float refX{}, refY{}, mx{}, my{};
0274
0275     float x{}, y{}, z{}, r{};
0276
0277     x = pS->m_n1->m_spGbts.SP->x();
0278     y = pS->m_n1->m_spGbts.SP->y();
0279     z = pS->m_n1->m_spGbts.SP->z();
0280     r = pS->m_n1->m_spGbts.r();
0281
0282     refX = x * ts.m_c + y * ts.m_s;
0283     mx = -x * ts.m_s + y * ts.m_c;  // measured X[0]
0284     refY = r;
0285     my = z;  // measured Y[0]
0286
0287     float A = refX - ts.m_refX;
0288     float B = 0.5 * A * A;
0289     float dr = refY - ts.m_refY;
0290
0291     X[0] = ts.m_X[0] + ts.m_X[1] * A + ts.m_X[2] * B;
0292     X[1] = ts.m_X[1] + ts.m_X[2] * A;
0293     X[2] = ts.m_X[2];
0294
0295     Cx[0][0] = ts.m_Cx[0][0] + 2 * ts.m_Cx[0][1] * A + 2 * ts.m_Cx[0][2] * B +
0296                A * A * ts.m_Cx[1][1] + 2 * A * B * ts.m_Cx[1][2] +
0297                B * B * ts.m_Cx[2][2];
0298     Cx[0][1] = Cx[1][0] = ts.m_Cx[0][1] + ts.m_Cx[1][1] * A +
0299                           ts.m_Cx[1][2] * B + ts.m_Cx[0][2] * A +
0300                           A * A * ts.m_Cx[1][2] + A * B * ts.m_Cx[2][2];
0301     Cx[0][2] = Cx[2][0] = ts.m_Cx[0][2] + ts.m_Cx[1][2] * A + ts.m_Cx[2][2] * B;
0302
0303     Cx[1][1] = ts.m_Cx[1][1] + 2 * A * ts.m_Cx[1][2] + A * A * ts.m_Cx[2][2];
0304     Cx[1][2] = Cx[2][1] = ts.m_Cx[1][2] + ts.m_Cx[2][2] * A;
0305
0306     Cx[2][2] = ts.m_Cx[2][2];
0307
0308     Y[0] = ts.m_Y[0] + ts.m_Y[1] * dr;
0309     Y[1] = ts.m_Y[1];
0310
0311     Cy[0][0] = ts.m_Cy[0][0] + 2 * ts.m_Cy[0][1] * dr + dr * dr * ts.m_Cy[1][1];
0312     Cy[0][1] = Cy[1][0] = ts.m_Cy[0][1] + dr * ts.m_Cy[1][1];
0313     Cy[1][1] = ts.m_Cy[1][1];
0314
0315     float resid_x = mx - X[0];
0316     float resid_y = my - Y[0];
0317
0318     float CHx[3] = {Cx[0][0], Cx[0][1], Cx[0][2]};
0319     float CHy[2] = {Cy[0][0], Cy[0][1]};
0320
0321     float sigma_rz = 0.0;
0322
0323     int type = getLayerType(pS->m_n1->m_spGbts.combined_ID);
0324
0325     if (type == 0) {  // barrel TODO: split into barrel Pixel and barrel SCT
0326       sigma_rz = sigma_y * sigma_y;
0327     } else {
0328       sigma_rz = sigma_y * ts.m_Y[1];
0329       sigma_rz = sigma_rz * sigma_rz;
0330     }
0331
0332     float Dx = 1.0 / (Cx[0][0] + sigma_x * sigma_x);
0333
0334     float Dy = 1.0 / (Cy[0][0] + sigma_rz);
0335
0336     float dchi2_x = resid_x * resid_x * Dx;
0337     float dchi2_y = resid_y * resid_y * Dy;
0338
0339     if (dchi2_x > maxDChi2_x || dchi2_y > maxDChi2_y) {
0340       return false;
0341     }
0342
0343     ts.m_J += add_hit - dchi2_x * weight_x - dchi2_y * weight_y;
0344
0345     // state update
0346     float Kx[3] = {Dx * Cx[0][0], Dx * Cx[0][1], Dx * Cx[0][2]};
0347     float Ky[2] = {Dy * Cy[0][0], Dy * Cy[0][1]};
0348
0349     for (int i = 0; i < 3; i++) {
0350       ts.m_X[i] = X[i] + Kx[i] * resid_x;
0351     }
0352     for (int i = 0; i < 2; i++) {
0353       ts.m_Y[i] = Y[i] + Ky[i] * resid_y;
0354     }
0355
0356     for (int i = 0; i < 3; i++) {
0357       for (int j = 0; j < 3; j++) {
0358         ts.m_Cx[i][j] = Cx[i][j] - Kx[i] * CHx[j];
0359       }
0360     }
0361
0362     for (int i = 0; i < 2; i++) {
0363       for (int j = 0; j < 2; j++) {
0364         ts.m_Cy[i][j] = Cy[i][j] - Ky[i] * CHy[j];
0365       }
0366     }
0367
0368     ts.m_refX = refX;
0369     ts.m_refY = refY;
0370
0371     return true;
0372   }
0373
0374   int getLayerType(int l) {
0375     auto iterator = find_if(m_geo.begin(), m_geo.end(), [l](auto n) {
0376       return n.m_subdet == l;
0377     });  // iterator to vector member with this id
0378     int index = std::distance(m_geo.begin(), iterator);
0379
0380     return m_geo.at(index).m_type;  // needs to be 0, 2, or -2
0381   }
0382
0383   const std::vector<TrigInDetSiLayer>& m_geo;
0384
0385   std::vector<GbtsEdge<external_spacepoint_t>>& m_segStore;
0386
0387   std::vector<GbtsEdgeState<external_spacepoint_t>*> m_stateVec;
0388
0389   GbtsEdgeState<external_spacepoint_t> m_stateStore[MAX_EDGE_STATE];
0390
0391   int m_globalStateCounter{0};
0392
0393   const Acts::Logger& logger() const { return *m_logger; }
0394   std::unique_ptr<const Acts::Logger> m_logger{nullptr};
0395 };
0396
0397 }  // namespace Acts::Experimental