Acts/Seeding/GbtsTrackingFilter.hpp

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