Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-02-14 04:12:07

0001 #include "MkFinder.h"
0002 
0003 #include "RecoTracker/MkFitCore/interface/cms_common_macros.h"
0004 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0005 #include "CandCloner.h"
0006 #include "FindingFoos.h"
0007 #include "KalmanUtilsMPlex.h"
0008 #include "MatriplexPackers.h"
0009 #include "MiniPropagators.h"
0010 
0011 //#define DEBUG
0012 #include "Debug.h"
0013 
0014 #if defined(MKFIT_STANDALONE)
0015 #include "RecoTracker/MkFitCore/standalone/Event.h"
0016 #endif
0017 
0018 #ifdef RNT_DUMP_MkF_SelHitIdcs
0019 // declares struct RntIfc_selectHitIndices rnt_shi in unnamed namespace;
0020 #include "RecoTracker/MkFitCore/standalone/RntDumper/MkFinder_selectHitIndices.icc"
0021 #endif
0022 
0023 #include "vdt/atan2.h"
0024 
0025 #include <algorithm>
0026 #include <queue>
0027 
0028 namespace mkfit {
0029 
0030   void MkFinder::setup(const PropagationConfig &pc,
0031                        const IterationConfig &ic,
0032                        const IterationParams &ip,
0033                        const IterationLayerConfig &ilc,
0034                        const SteeringParams &sp,
0035                        const std::vector<bool> *ihm,
0036                        const Event *ev,
0037                        int region,
0038                        bool infwd) {
0039     m_prop_config = &pc;
0040     m_iteration_config = &ic;
0041     m_iteration_params = &ip;
0042     m_iteration_layer_config = &ilc;
0043     m_steering_params = &sp;
0044     m_iteration_hit_mask = ihm;
0045     m_event = ev;
0046     m_current_region = region;
0047     m_in_fwd = infwd;
0048   }
0049 
0050   void MkFinder::setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp, const Event *ev) {
0051     m_prop_config = &pc;
0052     m_steering_params = &sp;
0053     m_event = ev;
0054   }
0055 
0056   void MkFinder::release() {
0057     m_prop_config = nullptr;
0058     m_iteration_config = nullptr;
0059     m_iteration_params = nullptr;
0060     m_iteration_layer_config = nullptr;
0061     m_steering_params = nullptr;
0062     m_iteration_hit_mask = nullptr;
0063     m_event = nullptr;
0064     m_current_region = -1;
0065     m_in_fwd = true;
0066   }
0067 
0068   void MkFinder::begin_layer(const LayerOfHits &layer_of_hits) {
0069 #ifdef RNT_DUMP_MkF_SelHitIdcs
0070     const LayerOfHits &L = layer_of_hits;
0071     const LayerInfo &LI = *L.layer_info();
0072     rnt_shi.ResetH();
0073     rnt_shi.ResetF();
0074     *rnt_shi.h = {m_event->evtID(),
0075                   m_iteration_config->m_iteration_index,
0076                   m_iteration_config->m_track_algorithm,
0077                   m_current_region,
0078                   L.layer_id(),
0079                   L.is_barrel() ? LI.rin() : LI.zmin(),
0080                   LI.is_barrel() ? LI.rout() : LI.zmax(),
0081                   L.is_barrel(),
0082                   L.is_pixel(),
0083                   L.is_stereo()};
0084     *rnt_shi.f = *rnt_shi.h;
0085 #endif
0086   }
0087 
0088   void MkFinder::end_layer() {
0089 #ifdef RNT_DUMP_MkF_SelHitIdcs
0090     rnt_shi.FillH();
0091     rnt_shi.FillF();
0092 #endif
0093   }
0094 
0095   //==============================================================================
0096   // Input / Output TracksAndHitIdx
0097   //==============================================================================
0098 
0099   void MkFinder::inputTracksAndHitIdx(const std::vector<Track> &tracks, int beg, int end, bool inputProp) {
0100     // Assign track parameters to initial state and copy hit values in.
0101 
0102     // This might not be true for the last chunk!
0103     // assert(end - beg == NN);
0104 
0105     const int iI = inputProp ? iP : iC;
0106 
0107     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0108       copy_in(tracks[i], imp, iI);
0109     }
0110   }
0111 
0112   void MkFinder::inputTracksAndHitIdx(
0113       const std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool inputProp, int mp_offset) {
0114     // Assign track parameters to initial state and copy hit values in.
0115 
0116     // This might not be true for the last chunk!
0117     // assert(end - beg == NN);
0118 
0119     const int iI = inputProp ? iP : iC;
0120 
0121     for (int i = beg, imp = mp_offset; i < end; ++i, ++imp) {
0122       copy_in(tracks[idxs[i]], imp, iI);
0123     }
0124   }
0125 
0126   void MkFinder::inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0127                                       const std::vector<std::pair<int, int>> &idxs,
0128                                       int beg,
0129                                       int end,
0130                                       bool inputProp) {
0131     // Assign track parameters to initial state and copy hit values in.
0132 
0133     // This might not be true for the last chunk!
0134     // assert(end - beg == NN);
0135 
0136     const int iI = inputProp ? iP : iC;
0137 
0138     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0139       const TrackCand &trk = tracks[idxs[i].first][idxs[i].second];
0140 
0141       copy_in(trk, imp, iI);
0142 
0143       m_SeedIdx(imp, 0, 0) = idxs[i].first;
0144       m_CandIdx(imp, 0, 0) = idxs[i].second;
0145       m_SeedOriginIdx[imp] = tracks[idxs[i].first].seed_origin_index();
0146     }
0147   }
0148 
0149   void MkFinder::inputTracksAndHits(const std::vector<CombCandidate> &tracks,
0150                                     const LayerOfHits &layer_of_hits,
0151                                     const std::vector<UpdateIndices> &idxs,
0152                                     int beg,
0153                                     int end,
0154                                     bool inputProp) {
0155     // Assign track parameters to initial state and copy hit values in.
0156 
0157     // This might not be true for the last chunk!
0158     // assert(end - beg == NN);
0159 
0160     const int iI = inputProp ? iP : iC;
0161 
0162     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0163       const TrackCand &trk = tracks[idxs[i].seed_idx][idxs[i].cand_idx];
0164 
0165       copy_in(trk, imp, iI);
0166 
0167       m_SeedIdx(imp, 0, 0) = idxs[i].seed_idx;
0168       m_CandIdx(imp, 0, 0) = idxs[i].cand_idx;
0169       m_SeedOriginIdx[imp] = tracks[idxs[i].seed_idx].seed_origin_index();
0170 
0171       // Reuse selectHitIndices() arrays -- used also in packModuleNormDir()
0172       m_XHitArr(imp, 0, 0) = idxs[i].hit_idx;
0173       m_XHitSize(imp, 0, 0) = 1;
0174 
0175       const Hit &hit = layer_of_hits.refHit(idxs[i].hit_idx);
0176       m_msErr.copyIn(imp, hit.errArray());
0177       m_msPar.copyIn(imp, hit.posArray());
0178     }
0179   }
0180 
0181   void MkFinder::inputOverlapHits(const LayerOfHits &layer_of_hits,
0182                                   const std::vector<UpdateIndices> &idxs,
0183                                   int beg,
0184                                   int end) {
0185     // Copy overlap hit values in.
0186 
0187     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0188       const Hit &hit = layer_of_hits.refHit(idxs[i].ovlp_idx);
0189       m_msErr.copyIn(imp, hit.errArray());
0190       m_msPar.copyIn(imp, hit.posArray());
0191     }
0192   }
0193 
0194   void MkFinder::inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0195                                       const std::vector<std::pair<int, IdxChi2List>> &idxs,
0196                                       int beg,
0197                                       int end,
0198                                       bool inputProp) {
0199     // Assign track parameters to initial state and copy hit values in.
0200 
0201     // This might not be true for the last chunk!
0202     // assert(end - beg == NN);
0203 
0204     const int iI = inputProp ? iP : iC;
0205 
0206     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0207       const TrackCand &trk = tracks[idxs[i].first][idxs[i].second.trkIdx];
0208 
0209       copy_in(trk, imp, iI);
0210 
0211       m_SeedIdx(imp, 0, 0) = idxs[i].first;
0212       m_CandIdx(imp, 0, 0) = idxs[i].second.trkIdx;
0213       m_SeedOriginIdx[imp] = tracks[idxs[i].first].seed_origin_index();
0214     }
0215   }
0216 
0217   void MkFinder::outputTracksAndHitIdx(std::vector<Track> &tracks, int beg, int end, bool outputProp) const {
0218     // Copies requested track parameters into Track objects.
0219     // The tracks vector should be resized to allow direct copying.
0220 
0221     const int iO = outputProp ? iP : iC;
0222 
0223     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0224       copy_out(tracks[i], imp, iO);
0225     }
0226   }
0227 
0228   void MkFinder::outputTracksAndHitIdx(
0229       std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool outputProp) const {
0230     // Copies requested track parameters into Track objects.
0231     // The tracks vector should be resized to allow direct copying.
0232 
0233     const int iO = outputProp ? iP : iC;
0234 
0235     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0236       copy_out(tracks[idxs[i]], imp, iO);
0237     }
0238   }
0239 
0240   void MkFinder::packModuleNormDir(
0241       const LayerOfHits &layer_of_hits, int hit_cnt, MPlexHV &norm, MPlexHV &dir, int N_proc) const {
0242     for (int itrack = 0; itrack < NN; ++itrack) {
0243       if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
0244         const auto &hit = layer_of_hits.refHit(m_XHitArr.constAt(itrack, hit_cnt, 0));
0245         unsigned int mid = hit.detIDinLayer();
0246         const ModuleInfo &mi = layer_of_hits.layer_info()->module_info(mid);
0247         norm.At(itrack, 0, 0) = mi.zdir[0];
0248         norm.At(itrack, 1, 0) = mi.zdir[1];
0249         norm.At(itrack, 2, 0) = mi.zdir[2];
0250         dir.At(itrack, 0, 0) = mi.xdir[0];
0251         dir.At(itrack, 1, 0) = mi.xdir[1];
0252         dir.At(itrack, 2, 0) = mi.xdir[2];
0253       }
0254     }
0255   }
0256 
0257   //==============================================================================
0258   // getHitSelDynamicWindows
0259   //==============================================================================
0260   // From HitSelectionWindows.h: track-related config on hit selection windows
0261 
0262   void MkFinder::getHitSelDynamicWindows(
0263       const float invpt, const float theta, float &min_dq, float &max_dq, float &min_dphi, float &max_dphi) {
0264     float max_invpt = std::min(invpt, 10.0f);  // => pT>0.1 GeV
0265 
0266     enum SelWinParameters_e { dp_sf = 0, dp_0, dp_1, dp_2, dq_sf, dq_0, dq_1, dq_2 };
0267     auto &v = m_iteration_layer_config->get_window_params(m_in_fwd, true);
0268 
0269     if (!v.empty()) {
0270       // dq hit selection window
0271       const float this_dq = v[dq_sf] * (v[dq_0] * max_invpt + v[dq_1] * theta + v[dq_2]);
0272       // In case value is below 0 (bad window derivation or other reasons), leave original limits
0273       if (this_dq > 0.f) {
0274         min_dq = this_dq;
0275         max_dq = 2.0f * min_dq;
0276       }
0277 
0278       // dphi hit selection window
0279       const float this_dphi = v[dp_sf] * (v[dp_0] * max_invpt + v[dp_1] * theta + v[dp_2]);
0280       // In case value is too low (bad window derivation or other reasons), leave original limits
0281       if (this_dphi > min_dphi) {
0282         min_dphi = this_dphi;
0283         max_dphi = 2.0f * min_dphi;
0284       }
0285     }
0286   }
0287 
0288   //==============================================================================
0289   // getHitSelDynamicChi2Cut
0290   //==============================================================================
0291   // From HitSelectionWindows.h: track-related config on hit selection windows
0292 
0293   inline float MkFinder::getHitSelDynamicChi2Cut(const int itrk, const int ipar) {
0294     const float minChi2Cut = m_iteration_params->chi2Cut_min;
0295     const float invpt = m_Par[ipar].At(itrk, 3, 0);
0296     const float theta = std::abs(m_Par[ipar].At(itrk, 5, 0) - Const::PIOver2);
0297 
0298     float max_invpt = std::min(invpt, 10.0f);  // => pT>0.1 GeV
0299 
0300     enum SelWinParameters_e { c2_sf = 8, c2_0, c2_1, c2_2 };
0301     auto &v = m_iteration_layer_config->get_window_params(m_in_fwd, true);
0302 
0303     if (!v.empty()) {
0304       float this_c2 = v[c2_sf] * (v[c2_0] * max_invpt + v[c2_1] * theta + v[c2_2]);
0305       // In case value is too low (bad window derivation or other reasons), leave original limits
0306       if (this_c2 > minChi2Cut)
0307         return this_c2;
0308     }
0309     return minChi2Cut;
0310   }
0311 
0312   //==============================================================================
0313   // SelectHitIndices
0314   //==============================================================================
0315 
0316   void MkFinder::selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc, bool fill_binsearch_only) {
0317     // bool debug = true;
0318     using bidx_t = LayerOfHits::bin_index_t;
0319     using bcnt_t = LayerOfHits::bin_content_t;
0320     const LayerOfHits &L = layer_of_hits;
0321     const IterationLayerConfig &ILC = *m_iteration_layer_config;
0322 
0323     const int iI = iP;
0324     const float nSigmaPhi = 3;
0325     const float nSigmaZ = 3;
0326     const float nSigmaR = 3;
0327 
0328     dprintf("LayerOfHits::SelectHitIndices %s layer=%d N_proc=%d\n",
0329             L.is_barrel() ? "barrel" : "endcap",
0330             L.layer_id(),
0331             N_proc);
0332 
0333     float dqv[NN], dphiv[NN], qv[NN], phiv[NN];
0334     bidx_t qb1v[NN], qb2v[NN], qbv[NN], pb1v[NN], pb2v[NN];
0335 
0336     const auto assignbins = [&](int itrack,
0337                                 float q,
0338                                 float dq,
0339                                 float phi,
0340                                 float dphi,
0341                                 float min_dq,
0342                                 float max_dq,
0343                                 float min_dphi,
0344                                 float max_dphi) {
0345       dphi = std::clamp(std::abs(dphi), min_dphi, max_dphi);
0346       dq = std::clamp(dq, min_dq, max_dq);
0347       //
0348       qv[itrack] = q;
0349       phiv[itrack] = phi;
0350       dphiv[itrack] = dphi;
0351       dqv[itrack] = dq;
0352       //
0353       qbv[itrack] = L.qBinChecked(q);
0354       qb1v[itrack] = L.qBinChecked(q - dq);
0355       qb2v[itrack] = L.qBinChecked(q + dq) + 1;
0356       pb1v[itrack] = L.phiBinChecked(phi - dphi);
0357       pb2v[itrack] = L.phiMaskApply(L.phiBin(phi + dphi) + 1);
0358     };
0359 
0360     const auto calcdphi2 = [&](int itrack, float dphidx, float dphidy) {
0361       return dphidx * dphidx * m_Err[iI].constAt(itrack, 0, 0) + dphidy * dphidy * m_Err[iI].constAt(itrack, 1, 1) +
0362              2 * dphidx * dphidy * m_Err[iI].constAt(itrack, 0, 1);
0363     };
0364 
0365     const auto calcdphi = [&](float dphi2, float min_dphi) {
0366       return std::max(nSigmaPhi * std::sqrt(std::abs(dphi2)), min_dphi);
0367     };
0368 
0369     if (L.is_barrel()) {
0370       // Pull out the part of the loop that vectorizes with icc and gcc
0371       // In llvm16 clang issues a warning that it can't vectorize
0372       // the loop.  Unfortunately, there doesn't seem to be a
0373       // pragma to suppress the warning, so we ifdef it out.  This
0374       // should be rechecked if llvm vectorization improves.
0375 #if !defined(__clang__)
0376 #pragma omp simd
0377 #endif
0378       for (int itrack = 0; itrack < NN; ++itrack) {
0379         if (itrack >= N_proc)
0380           continue;
0381         m_XHitSize[itrack] = 0;
0382 
0383         float min_dq = ILC.min_dq();
0384         float max_dq = ILC.max_dq();
0385         float min_dphi = ILC.min_dphi();
0386         float max_dphi = ILC.max_dphi();
0387 
0388         const float invpt = m_Par[iI].At(itrack, 3, 0);
0389         const float theta = std::fabs(m_Par[iI].At(itrack, 5, 0) - Const::PIOver2);
0390         getHitSelDynamicWindows(invpt, theta, min_dq, max_dq, min_dphi, max_dphi);
0391 
0392         const float x = m_Par[iI].constAt(itrack, 0, 0);
0393         const float y = m_Par[iI].constAt(itrack, 1, 0);
0394         const float r2 = x * x + y * y;
0395         const float dphidx = -y / r2, dphidy = x / r2;
0396         const float dphi2 = calcdphi2(itrack, dphidx, dphidy);
0397 #ifdef HARD_CHECK
0398         assert(dphi2 >= 0);
0399 #endif
0400 
0401         const float phi = getPhi(x, y);
0402         float dphi = calcdphi(dphi2, min_dphi);
0403 
0404         const float z = m_Par[iI].constAt(itrack, 2, 0);
0405         const float dz = std::abs(nSigmaZ * std::sqrt(m_Err[iI].constAt(itrack, 2, 2)));
0406         const float edgeCorr = std::abs(0.5f * (L.layer_info()->rout() - L.layer_info()->rin()) /
0407                                         std::tan(m_Par[iI].constAt(itrack, 5, 0)));
0408         // XXX-NUM-ERR above, m_Err(2,2) gets negative!
0409 
0410         m_XWsrResult[itrack] = L.is_within_z_sensitive_region(z, std::sqrt(dz * dz + edgeCorr * edgeCorr));
0411         assignbins(itrack, z, dz, phi, dphi, min_dq, max_dq, min_dphi, max_dphi);
0412 
0413         // Relax propagation-fail detection to be in line with pre-43145.
0414         if (m_FailFlag[itrack] && std::sqrt(r2) >= L.layer_info()->rin()) {
0415           m_FailFlag[itrack] = 0;
0416         }
0417       }
0418     } else  // endcap
0419     {
0420       //layer half-thikness for dphi spread calculation; only for very restrictive iters
0421       const float layerD = std::abs(L.layer_info()->zmax() - L.layer_info()->zmin()) * 0.5f *
0422                            (m_iteration_params->maxConsecHoles == 0 || m_iteration_params->maxHolesPerCand == 0);
0423       // Pull out the part of the loop that vectorizes with icc and gcc
0424 #if !defined(__clang__)
0425 #pragma omp simd
0426 #endif
0427       for (int itrack = 0; itrack < NN; ++itrack) {
0428         m_XHitSize[itrack] = 0;
0429 
0430         float min_dq = ILC.min_dq();
0431         float max_dq = ILC.max_dq();
0432         float min_dphi = ILC.min_dphi();
0433         float max_dphi = ILC.max_dphi();
0434 
0435         const float invpt = m_Par[iI].At(itrack, 3, 0);
0436         const float theta = std::fabs(m_Par[iI].At(itrack, 5, 0) - Const::PIOver2);
0437         getHitSelDynamicWindows(invpt, theta, min_dq, max_dq, min_dphi, max_dphi);
0438 
0439         const float x = m_Par[iI].constAt(itrack, 0, 0);
0440         const float y = m_Par[iI].constAt(itrack, 1, 0);
0441         const float r2 = x * x + y * y;
0442         const float r2Inv = 1.f / r2;
0443         const float dphidx = -y * r2Inv, dphidy = x * r2Inv;
0444         const float phi = getPhi(x, y);
0445         const float dphi2 =
0446             calcdphi2(itrack, dphidx, dphidy)
0447             //range from finite layer thickness
0448             + std::pow(layerD * std::tan(m_Par[iI].At(itrack, 5, 0)) * std::sin(m_Par[iI].At(itrack, 4, 0) - phi), 2) *
0449                   r2Inv;
0450 #ifdef HARD_CHECK
0451         assert(dphi2 >= 0);
0452 #endif
0453 
0454         float dphi = calcdphi(dphi2, min_dphi);
0455 
0456         const float r = std::sqrt(r2);
0457         const float dr = nSigmaR * std::sqrt(std::abs(x * x * m_Err[iI].constAt(itrack, 0, 0) +
0458                                                       y * y * m_Err[iI].constAt(itrack, 1, 1) +
0459                                                       2 * x * y * m_Err[iI].constAt(itrack, 0, 1)) /
0460                                              r2);
0461         const float edgeCorr = std::abs(0.5f * (L.layer_info()->zmax() - L.layer_info()->zmin()) *
0462                                         std::tan(m_Par[iI].constAt(itrack, 5, 0)));
0463 
0464         m_XWsrResult[itrack] = L.is_within_r_sensitive_region(r, std::sqrt(dr * dr + edgeCorr * edgeCorr));
0465         assignbins(itrack, r, dr, phi, dphi, min_dq, max_dq, min_dphi, max_dphi);
0466       }
0467     }
0468 
0469 #ifdef RNT_DUMP_MkF_SelHitIdcs
0470     if (fill_binsearch_only) {
0471       // XXX loop over good indices (prepared in V2) and put in V1 BinSearch results
0472       for (auto i : rnt_shi.f_h_idcs) {
0473         CandInfo &ci = (*rnt_shi.ci)[rnt_shi.f_h_remap[i]];
0474         ci.bso = BinSearch({phiv[i],
0475                             dphiv[i],
0476                             qv[i],
0477                             dqv[i],
0478                             pb1v[i],
0479                             pb2v[i],
0480                             qb1v[i],
0481                             qb2v[i],
0482                             m_XWsrResult[i].m_wsr,
0483                             m_XWsrResult[i].m_in_gap,
0484                             false});
0485       }
0486       return;
0487     }
0488 #endif
0489 
0490     // Vectorizing this makes it run slower!
0491     //#pragma omp simd
0492     for (int itrack = 0; itrack < NN; ++itrack) {
0493       if (itrack >= N_proc) {
0494         continue;
0495       }
0496       // PROP-FAIL-ENABLE The following to be enabled when propagation failure
0497       // detection is properly implemented in propagate-to-R/Z.
0498       if (m_FailFlag[itrack]) {
0499         m_XWsrResult[itrack].m_wsr = WSR_Failed;
0500         continue;
0501       }
0502 
0503       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
0504         continue;
0505       }
0506 
0507       const bidx_t qb = qbv[itrack];
0508       const bidx_t qb1 = qb1v[itrack];
0509       const bidx_t qb2 = qb2v[itrack];
0510       const bidx_t pb1 = pb1v[itrack];
0511       const bidx_t pb2 = pb2v[itrack];
0512 
0513       const float q = qv[itrack];
0514       const float phi = phiv[itrack];
0515       const float dphi = dphiv[itrack];
0516       const float dq = dqv[itrack];
0517 
0518       // clang-format off
0519       dprintf("  %2d/%2d: %6.3f %6.3f %6.6f %7.5f %3u %3u %4u %4u\n",
0520               L.layer_id(), itrack, q, phi, dq, dphi,
0521               qb1, qb2, pb1, pb2);
0522       // clang-format on
0523 
0524 #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE)
0525       const auto ngr = [](float f) { return isFinite(f) ? f : -999.0f; };
0526 
0527       const int seed_lbl = m_event->currentSeed(m_SeedOriginIdx[itrack]).label();
0528       Event::SimLabelFromHits slfh = m_event->simLabelForCurrentSeed(m_SeedOriginIdx[itrack]);
0529       const int seed_mcid = (slfh.is_set() && slfh.good_frac() > 0.7f) ? slfh.label : -999999;
0530 #endif
0531 
0532       for (bidx_t qi = qb1; qi != qb2; ++qi) {
0533         for (bidx_t pi = pb1; pi != pb2; pi = L.phiMaskApply(pi + 1)) {
0534           // Limit to central Q-bin
0535           if (qi == qb && L.isBinDead(pi, qi) == true) {
0536             dprint("dead module for track in layer=" << L.layer_id() << " qb=" << qi << " pi=" << pi << " q=" << q
0537                                                      << " phi=" << phi);
0538             m_XWsrResult[itrack].m_in_gap = true;
0539           }
0540 
0541           // MT: The following line is the biggest hog (4% total run time).
0542           // This comes from cache misses, I presume.
0543           // It might make sense to make first loop to extract bin indices
0544           // and issue prefetches at the same time.
0545           // Then enter vectorized loop to actually collect the hits in proper order.
0546 
0547           //SK: ~20x1024 bin sizes give mostly 1 hit per bin. Commented out for 128 bins or less
0548           // #pragma nounroll
0549           auto pbi = L.phiQBinContent(pi, qi);
0550           for (bcnt_t hi = pbi.begin(); hi < pbi.end(); ++hi) {
0551             // MT: Access into m_hit_zs and m_hit_phis is 1% run-time each.
0552 
0553             const unsigned int hi_orig = L.getOriginalHitIndex(hi);
0554 
0555             if (m_iteration_hit_mask && (*m_iteration_hit_mask)[hi_orig]) {
0556               dprintf(
0557                   "Yay, denying masked hit on layer %u, hi %u, orig idx %u\n", L.layer_info()->layer_id(), hi, hi_orig);
0558               continue;
0559             }
0560 
0561             if (Config::usePhiQArrays) {
0562               if (m_XHitSize[itrack] >= MPlexHitIdxMax)
0563                 break;
0564 
0565               const float ddq = std::abs(q - L.hit_q(hi));
0566               const float ddphi = cdist(std::abs(phi - L.hit_phi(hi)));
0567 
0568               // clang-format off
0569               dprintf("     SHI %3u %4u %5u  %6.3f %6.3f %6.4f %7.5f   %s\n",
0570                       qi, pi, hi, L.hit_q(hi), L.hit_phi(hi),
0571                       ddq, ddphi, (ddq < dq && ddphi < dphi) ? "PASS" : "FAIL");
0572               // clang-format on
0573 
0574 #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE)
0575               // clang-format off
0576               MPlexQF thisOutChi2;
0577               {
0578                 const MCHitInfo &mchinfo = m_event->simHitsInfo_[L.refHit(hi).mcHitID()];
0579                 int mchid = mchinfo.mcTrackID();
0580                 int st_isfindable = 0;
0581                 int st_label = -999999;
0582                 int st_prodtype = 0;
0583                 int st_nhits = -1;
0584                 int st_charge = 0;
0585                 float st_r = -999.;
0586                 float st_z = -999.;
0587                 float st_pt = -999.;
0588                 float st_eta = -999.;
0589                 float st_phi = -999.;
0590                 if (mchid >= 0) {
0591                   Track simtrack = m_event->simTracks_[mchid];
0592                   st_isfindable = (int)simtrack.isFindable();
0593                   st_label = simtrack.label();
0594                   st_prodtype = (int)simtrack.prodType();
0595                   st_pt = simtrack.pT();
0596                   st_eta = simtrack.momEta();
0597                   st_phi = simtrack.momPhi();
0598                   st_nhits = simtrack.nTotalHits();
0599                   st_charge = simtrack.charge();
0600                   st_r = simtrack.posR();
0601                   st_z = simtrack.z();
0602                 }
0603 
0604                 const Hit &thishit = L.refHit(hi_orig);
0605                 m_msErr.copyIn(itrack, thishit.errArray());
0606                 m_msPar.copyIn(itrack, thishit.posArray());
0607 
0608                 MPlexQI propFail;
0609                 MPlexLV tmpPropPar;
0610                 const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(L.is_barrel());
0611                 (*fnd_foos.m_compute_chi2_foo)(m_Err[iI],
0612                                                m_Par[iI],
0613                                                m_Chg,
0614                                                m_msErr,
0615                                                m_msPar,
0616                                                thisOutChi2,
0617                                                tmpPropPar,
0618                                                propFail,
0619                                                N_proc,
0620                                                m_prop_config->finding_intra_layer_pflags,
0621                                                m_prop_config->finding_requires_propagation_to_hit_pos);
0622 
0623                 float hx = thishit.x();
0624                 float hy = thishit.y();
0625                 float hz = thishit.z();
0626                 float hr = std::hypot(hx, hy);
0627                 float hphi = std::atan2(hy, hx);
0628                 float hex = ngr( std::sqrt(thishit.exx()) );
0629                 float hey = ngr( std::sqrt(thishit.eyy()) );
0630                 float hez = ngr( std::sqrt(thishit.ezz()) );
0631                 float her = ngr( std::sqrt(
0632                     (hx * hx * thishit.exx() + hy * hy * thishit.eyy() + 2.0f * hx * hy * m_msErr.At(itrack, 0, 1)) /
0633                     (hr * hr)) );
0634                 float hephi = ngr( std::sqrt(thishit.ephi()) );
0635                 float hchi2 = ngr( thisOutChi2[itrack] );
0636                 float tx = m_Par[iI].At(itrack, 0, 0);
0637                 float ty = m_Par[iI].At(itrack, 1, 0);
0638                 float tz = m_Par[iI].At(itrack, 2, 0);
0639                 float tr = std::hypot(tx, ty);
0640                 float tphi = std::atan2(ty, tx);
0641                 // float tchi2 = ngr( m_Chi2(itrack, 0, 0) ); // unused
0642                 float tex = ngr( std::sqrt(m_Err[iI].At(itrack, 0, 0)) );
0643                 float tey = ngr( std::sqrt(m_Err[iI].At(itrack, 1, 1)) );
0644                 float tez = ngr( std::sqrt(m_Err[iI].At(itrack, 2, 2)) );
0645                 float ter = ngr( std::sqrt(
0646                     (tx * tx * tex * tex + ty * ty * tey * tey + 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) /
0647                     (tr * tr)) );
0648                 float tephi = ngr( std::sqrt(
0649                     (ty * ty * tex * tex + tx * tx * tey * tey - 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) /
0650                     (tr * tr * tr * tr)) );
0651                 float ht_dxy = std::hypot(hx - tx, hy - ty);
0652                 float ht_dz = hz - tz;
0653                 float ht_dphi = cdist(std::abs(hphi - tphi));
0654 
0655                 static bool first = true;
0656                 if (first) {
0657                   printf(
0658                       "HITWINDOWSEL "
0659                       "evt_id/I:track_algo/I:"
0660                       "lyr_id/I:lyr_isbrl/I:hit_idx/I:"
0661                       "trk_cnt/I:trk_idx/I:trk_label/I:"
0662                       "trk_pt/F:trk_eta/F:trk_mphi/F:trk_chi2/F:"
0663                       "nhits/I:"
0664                       "seed_idx/I:seed_label/I:seed_mcid/I:"
0665                       "hit_mcid/I:"
0666                       "st_isfindable/I:st_prodtype/I:st_label/I:"
0667                       "st_pt/F:st_eta/F:st_phi/F:"
0668                       "st_nhits/I:st_charge/I:st_r/F:st_z/F:"
0669                       "trk_q/F:hit_q/F:dq_trkhit/F:dq_cut/F:trk_phi/F:hit_phi/F:dphi_trkhit/F:dphi_cut/F:"
0670                       "t_x/F:t_y/F:t_r/F:t_phi/F:t_z/F:"
0671                       "t_ex/F:t_ey/F:t_er/F:t_ephi/F:t_ez/F:"
0672                       "h_x/F:h_y/F:h_r/F:h_phi/F:h_z/F:"
0673                       "h_ex/F:h_ey/F:h_er/F:h_ephi/F:h_ez/F:"
0674                       "ht_dxy/F:ht_dz/F:ht_dphi/F:"
0675                       "h_chi2/F"
0676                       "\n");
0677                   first = false;
0678                 }
0679 
0680                 if (!(std::isnan(phi)) && !(std::isnan(getEta(m_Par[iI].At(itrack, 5, 0))))) {
0681                   //|| std::isnan(ter) || std::isnan(her) || std::isnan(m_Chi2(itrack, 0, 0)) || std::isnan(hchi2)))
0682                   printf("HITWINDOWSEL "
0683                          "%d %d"
0684                          "%d %d %d "
0685                          "%d %d %d "
0686                          "%6.3f %6.3f %6.3f %6.3f "
0687                          "%d "
0688                          "%d %d %d "
0689                          "%d "
0690                          "%d %d %d "
0691                          "%6.3f %6.3f %6.3f "
0692                          "%d %d %6.3f %6.3f "
0693                          "%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f "
0694                          "%6.3f %6.3f %6.3f %6.3f %6.3f "
0695                          "%6.6f %6.6f %6.6f %6.6f %6.6f "
0696                          "%6.3f %6.3f %6.3f %6.3f %6.3f "
0697                          "%6.6f %6.6f %6.6f %6.6f %6.6f "
0698                          "%6.3f %6.3f %6.3f "
0699                          "%6.3f"
0700                          "\n",
0701                          m_event->evtID(), m_iteration_config->m_track_algorithm,
0702                          L.layer_id(), L.is_barrel(), hi_orig,
0703                          itrack, m_CandIdx(itrack, 0, 0), m_Label(itrack, 0, 0),
0704                          1.0f / m_Par[iI].At(itrack, 3, 0), getEta(m_Par[iI].At(itrack, 5, 0)), m_Par[iI].At(itrack, 4, 0), m_Chi2(itrack, 0, 0),
0705                          m_NFoundHits(itrack, 0, 0),
0706                          m_SeedOriginIdx(itrack, 0, 0), seed_lbl, seed_mcid,
0707                          mchid,
0708                          st_isfindable, st_prodtype, st_label,
0709                          st_pt, st_eta, st_phi,
0710                          st_nhits, st_charge, st_r, st_z,
0711                          q, L.hit_q(hi), ddq, dq, phi, L.hit_phi(hi), ddphi, dphi,
0712                          tx, ty, tr, tphi, tz,
0713                          tex, tey, ter, tephi, tez,
0714                          hx, hy, hr, hphi, hz,
0715                          hex, hey, her, hephi, hez,
0716                          ht_dxy, ht_dz, ht_dphi,
0717                          hchi2);
0718                 }
0719               }
0720               // clang-format on
0721 #endif
0722 
0723               if (ddq >= dq)
0724                 continue;
0725               if (ddphi >= dphi)
0726                 continue;
0727 
0728               // MT: Removing extra check gives full efficiency ...
0729               //     and means our error estimations are wrong!
0730               // Avi says we should have *minimal* search windows per layer.
0731               // Also ... if bins are sufficiently small, we do not need the extra
0732               // checks, see above.
0733               m_XHitArr.At(itrack, m_XHitSize[itrack]++, 0) = hi_orig;
0734             } else {
0735               // MT: The following check alone makes more sense with spiral traversal,
0736               // we'd be taking in closest hits first.
0737 
0738               // Hmmh -- there used to be some more checks here.
0739               // Or, at least, the phi binning was much smaller and no further checks were done.
0740               assert(false && "this code has not been used in a while -- see comments in code");
0741 
0742               if (m_XHitSize[itrack] < MPlexHitIdxMax) {
0743                 m_XHitArr.At(itrack, m_XHitSize[itrack]++, 0) = hi_orig;
0744               }
0745             }
0746           }  //hi
0747         }    //pi
0748       }      //qi
0749     }        //itrack
0750   }
0751 
0752   //==============================================================================
0753   // SelectHitIndicesV2
0754   //==============================================================================
0755 
0756   void MkFinder::selectHitIndicesV2(const LayerOfHits &layer_of_hits, const int N_proc) {
0757     // bool debug = true;
0758     using bidx_t = LayerOfHits::bin_index_t;
0759     using bcnt_t = LayerOfHits::bin_content_t;
0760     const LayerOfHits &L = layer_of_hits;
0761     const LayerInfo &LI = *L.layer_info();
0762 
0763     const int iI = iP;
0764 
0765     dprintf("LayerOfHits::SelectHitIndicesV2 %s layer=%d N_proc=%d\n",
0766             L.is_barrel() ? "barrel" : "endcap",
0767             L.layer_id(),
0768             N_proc);
0769 
0770 #ifdef RNT_DUMP_MkF_SelHitIdcs
0771     rnt_shi.InnerIdcsReset(N_proc);
0772     for (int i = 0; i < N_proc; ++i) {
0773       auto slfh = m_event->simLabelForCurrentSeed(m_SeedOriginIdx[i]);
0774       if (m_FailFlag[i]) {
0775         rnt_shi.RegisterFailedProp(i, m_Par[1 - iI], m_Par[iI], m_event, m_SeedOriginIdx[i]);
0776       } else if (slfh.is_set()) {
0777         rnt_shi.RegisterGoodProp(i, m_Par[iI], m_event, m_SeedOriginIdx[i]);
0778         // get BinSearch result from V1.
0779         selectHitIndices(layer_of_hits, N_proc, true);
0780       }  // else ... could do something about the bad seeds ... probably better to collect elsewhere.
0781     }
0782 #endif
0783 
0784     constexpr int NEW_MAX_HIT = 6;  // 4 - 6 give about the same # of tracks in quality-val
0785     constexpr float DDPHI_PRESEL_FAC = 2.0f;
0786     constexpr float DDQ_PRESEL_FAC = 1.2f;
0787     constexpr float PHI_BIN_EXTRA_FAC = 2.75f;
0788     constexpr float Q_BIN_EXTRA_FAC = 1.6f;
0789 
0790     namespace mp = mini_propagators;
0791     struct Bins {
0792       MPlexQUH q0, q1, q2, p1, p2;
0793       mp::InitialStatePlex isp;
0794       mp::StatePlex sp1, sp2;
0795       int n_proc;
0796 
0797       MPlexQF dphi_track, dq_track;  // 3 sigma track errors at initial state
0798 
0799       // debug & ntuple dump -- to be local in functions
0800       MPlexQF phi_c, dphi;
0801       MPlexQF q_c, qmin, qmax;
0802 
0803       Bins(const MPlexLV &par, const MPlexQI &chg, int np = NN) : isp(par, chg), n_proc(np) {}
0804 
0805       void prop_to_limits(const LayerInfo &li) {
0806         // Positions 1 and 2 should really be by "propagation order", 1 is the closest/
0807         // This should also work for backward propagation so not exactly trivial.
0808         // Also, do not really need propagation to center.
0809         if (li.is_barrel()) {
0810           isp.propagate_to_r(mp::PA_Exact, li.rin(), sp1, true, n_proc);
0811           isp.propagate_to_r(mp::PA_Exact, li.rout(), sp2, true, n_proc);
0812         } else {
0813           isp.propagate_to_z(mp::PA_Exact, li.zmin(), sp1, true, n_proc);
0814           isp.propagate_to_z(mp::PA_Exact, li.zmax(), sp2, true, n_proc);
0815         }
0816       }
0817 
0818       void find_bin_ranges(const LayerInfo &li, const LayerOfHits &loh, const MPlexLS &err) {
0819         // Below made members for debugging
0820         // MPlexQF phi_c, dphi_min, dphi_max;
0821         // phi_c = mp::fast_atan2(isp.y, isp.x);  // calculated below as difference
0822 
0823         // Matriplex::min_max(sp1.dphi, sp2.dphi, dphi_min, dphi_max);
0824         // the above is wrong: dalpha is not dphi --> renamed variable in State
0825         MPlexQF xp1, xp2, pmin, pmax;
0826         xp1 = mp::fast_atan2(sp1.y, sp1.x);
0827         xp2 = mp::fast_atan2(sp2.y, sp2.x);
0828         Matriplex::min_max(xp1, xp2, pmin, pmax);
0829         // Matriplex::min_max(mp::fast_atan2(sp1.y, sp1.x), smp::fast_atan2(sp2.y, sp2.x), pmin, pmax);
0830         MPlexQF dp = pmax - pmin;
0831         phi_c = 0.5f * (pmax + pmin);
0832         for (int ii = 0; ii < NN; ++ii) {
0833           if (ii < n_proc) {
0834             if (dp[ii] > Const::PI) {
0835               std::swap(pmax[ii], pmin[ii]);
0836               dp[ii] = Const::TwoPI - dp[ii];
0837               phi_c[ii] = Const::PI - phi_c[ii];
0838             }
0839             dphi[ii] = 0.5f * dp[ii];
0840             // printf("phic: %f  p1: %f  p2: %f   pmin: %f  pmax: %f   dphi: %f\n",
0841             //       phi_c[ii], xp1[ii], xp2[ii], pmin[ii], pmax[ii], dphi[ii]);
0842           }
0843         }
0844 
0845         const auto calc_err_xy = [&](const MPlexQF &x, const MPlexQF &y) {
0846           return x * x * err.ReduceFixedIJ(0, 0) + y * y * err.ReduceFixedIJ(1, 1) +
0847                  2.0f * x * y * err.ReduceFixedIJ(0, 1);
0848         };
0849 
0850         // Calculate dphi_track, dq_track differs for barrel/endcap
0851         MPlexQF r2_c = isp.x * isp.x + isp.y * isp.y;
0852         MPlexQF r2inv_c = 1.0f / r2_c;
0853         MPlexQF dphidx_c = -isp.y * r2inv_c;
0854         MPlexQF dphidy_c = isp.x * r2inv_c;
0855         dphi_track = 3.0f * calc_err_xy(dphidx_c, dphidy_c).abs().sqrt();
0856 
0857         // MPlexQF qmin, qmax;
0858         if (li.is_barrel()) {
0859           Matriplex::min_max(sp1.z, sp2.z, qmin, qmax);
0860           q_c = isp.z;
0861           dq_track = 3.0f * err.ReduceFixedIJ(2, 2).abs().sqrt();
0862         } else {
0863           Matriplex::min_max(Matriplex::hypot(sp1.x, sp1.y), Matriplex::hypot(sp2.x, sp2.y), qmin, qmax);
0864           q_c = Matriplex::sqrt(r2_c);
0865           dq_track = 3.0f * (r2inv_c * calc_err_xy(isp.x, isp.y).abs()).sqrt();
0866         }
0867 
0868         for (int i = 0; i < p1.kTotSize; ++i) {
0869           // Clamp crazy sizes. This actually only happens when prop-fail flag is set.
0870           // const float dphi_clamp = 0.1;
0871           // if (dphi_min[i] > 0.0f || dphi_min[i] < -dphi_clamp) dphi_min[i] = -dphi_clamp;
0872           // if (dphi_max[i] < 0.0f || dphi_max[i] > dphi_clampf) dphi_max[i] = dphi_clamp;
0873           p1[i] = loh.phiBinChecked(pmin[i] - dphi_track[i] - PHI_BIN_EXTRA_FAC * 0.0123f);
0874           p2[i] = loh.phiBinChecked(pmax[i] + dphi_track[i] + PHI_BIN_EXTRA_FAC * 0.0123f);
0875 
0876           q0[i] = loh.qBinChecked(q_c[i]);
0877           q1[i] = loh.qBinChecked(qmin[i] - dq_track[i] - Q_BIN_EXTRA_FAC * 0.5f * li.q_bin());
0878           q2[i] = loh.qBinChecked(qmax[i] + dq_track[i] + Q_BIN_EXTRA_FAC * 0.5f * li.q_bin()) + 1;
0879         }
0880       }
0881     };
0882 
0883     Bins B(m_Par[iI], m_Chg, N_proc);
0884     B.prop_to_limits(LI);
0885     B.find_bin_ranges(LI, L, m_Err[iI]);
0886 
0887     for (int i = 0; i < NN; ++i) {
0888       if (i < N_proc) {
0889         m_XHitSize[i] = 0;
0890         // Notify failure. Ideally should be detected before selectHitIndices().
0891         if (m_FailFlag[i]) {
0892           m_XWsrResult[i].m_wsr = WSR_Failed;
0893         } else {
0894           if (LI.is_barrel()) {
0895             m_XWsrResult[i] = L.is_within_z_sensitive_region(B.q_c[i], 0.5f * (B.q2[i] - B.q1[i]));
0896           } else {
0897             m_XWsrResult[i] = L.is_within_r_sensitive_region(B.q_c[i], 0.5f * (B.q2[i] - B.q1[i]));
0898           }
0899         }
0900       }
0901     }
0902 
0903     // for (int i = 0; i < N_proc; ++i) {
0904     //   printf("BinCheck %c %+8.6f %+8.6f | %3d %3d - %3d %3d ||  | %2d %2d - %2d %2d\n", LI.is_barrel() ? 'B' : 'E',
0905     //          B.phi[i], B.dphi[i], B.p1[i], B.p2[i], pb1v[i], pb2v[i],
0906     //          B.q[i], B.dq[i], B.q1[i], B.q2[i], qb1v[i], qb2v[i]);
0907     // }
0908 
0909 #ifdef RNT_DUMP_MkF_SelHitIdcs
0910     for (auto i : rnt_shi.f_h_idcs) {
0911       CandInfo &ci = (*rnt_shi.ci)[rnt_shi.f_h_remap[i]];
0912       ci.bsn = BinSearch({B.phi_c[i],
0913                           B.dphi[i],
0914                           B.q_c[i],
0915                           0.5f * (B.q2[i] - B.q1[i]),
0916                           B.p1[i],
0917                           B.p2[i],
0918                           B.q1[i],
0919                           B.q2[i],
0920                           m_XWsrResult[i].m_wsr,
0921                           m_XWsrResult[i].m_in_gap,
0922                           false});
0923       ci.ps_min = statep2propstate(B.sp1, i);
0924       ci.ps_max = statep2propstate(B.sp2, i);
0925     }
0926 #endif
0927 
0928     struct PQE {
0929       float score;
0930       unsigned int hit_index;
0931     };
0932     auto pqe_cmp = [](const PQE &a, const PQE &b) { return a.score < b.score; };
0933     std::priority_queue<PQE, std::vector<PQE>, decltype(pqe_cmp)> pqueue(pqe_cmp);
0934     int pqueue_size = 0;
0935 
0936     // Vectorizing this makes it run slower!
0937     //#pragma omp simd
0938     for (int itrack = 0; itrack < NN; ++itrack) {
0939       if (itrack >= N_proc) {
0940         continue;
0941       }
0942 
0943       if (m_FailFlag[itrack]) {
0944         m_XWsrResult[itrack].m_wsr = WSR_Failed;
0945         continue;
0946       }
0947 
0948       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
0949         continue;
0950       }
0951 
0952       // New binning -- known to be too restrictive, so scaled up. Probably esp. in stereo layers.
0953       // Also, could take track covariance dphi / dq extras + known tilt stuff.
0954       const bidx_t qb = B.q0[itrack];
0955       const bidx_t qb1 = B.q1[itrack];
0956       const bidx_t qb2 = B.q2[itrack];
0957       const bidx_t pb1 = B.p1[itrack];
0958       const bidx_t pb2 = B.p2[itrack];
0959 
0960       // clang-format off
0961       dprintf("  %2d/%2d: %6.3f %6.3f %6.6f %7.5f %3u %3u %4u %4u\n",
0962               L.layer_id(), itrack, qv[itrack], phi[itrack], dqv[itrack], dphiv[itrack],
0963               qb1, qb2, pb1, pb2);
0964       // clang-format on
0965 
0966       mp::InitialState mp_is(m_Par[iI], m_Chg, itrack);
0967       mp::State mp_s;
0968 
0969       for (bidx_t qi = qb1; qi != qb2; ++qi) {
0970         for (bidx_t pi = pb1; pi != pb2; pi = L.phiMaskApply(pi + 1)) {
0971           // Limit to central Q-bin
0972           if (qi == qb && L.isBinDead(pi, qi) == true) {
0973             dprint("dead module for track in layer=" << L.layer_id() << " qb=" << qi << " pi=" << pi << " q=" << q
0974                                                      << " phi=" << phi);
0975             m_XWsrResult[itrack].m_in_gap = true;
0976           }
0977 
0978           // It might make sense to make first loop to extract bin indices
0979           // and issue prefetches at the same time.
0980           // Then enter vectorized loop to actually collect the hits in proper order.
0981 
0982           //SK: ~20x1024 bin sizes give mostly 1 hit per bin. Commented out for 128 bins or less
0983           // #pragma nounroll
0984           auto pbi = L.phiQBinContent(pi, qi);
0985           for (bcnt_t hi = pbi.begin(); hi < pbi.end(); ++hi) {
0986             // MT: Access into m_hit_zs and m_hit_phis is 1% run-time each.
0987 
0988             const unsigned int hi_orig = L.getOriginalHitIndex(hi);
0989 
0990             if (m_iteration_hit_mask && (*m_iteration_hit_mask)[hi_orig]) {
0991               dprintf(
0992                   "Yay, denying masked hit on layer %u, hi %u, orig idx %u\n", L.layer_info()->layer_id(), hi, hi_orig);
0993               continue;
0994             }
0995 
0996             if (m_XHitSize[itrack] >= MPlexHitIdxMax)
0997               break;
0998 
0999             float new_q, new_phi, new_ddphi, new_ddq;
1000             bool prop_fail;
1001 
1002             if (L.is_barrel()) {
1003               prop_fail = mp_is.propagate_to_r(mp::PA_Exact, L.hit_qbar(hi), mp_s, true);
1004               new_q = mp_s.z;
1005             } else {
1006               prop_fail = mp_is.propagate_to_z(mp::PA_Exact, L.hit_qbar(hi), mp_s, true);
1007               new_q = std::hypot(mp_s.x, mp_s.y);
1008             }
1009 
1010             new_phi = vdt::fast_atan2f(mp_s.y, mp_s.x);
1011             new_ddphi = cdist(std::abs(new_phi - L.hit_phi(hi)));
1012             new_ddq = std::abs(new_q - L.hit_q(hi));
1013 
1014             bool dqdphi_presel = new_ddq < B.dq_track[itrack] + DDQ_PRESEL_FAC * L.hit_q_half_length(hi) &&
1015                                  new_ddphi < B.dphi_track[itrack] + DDPHI_PRESEL_FAC * 0.0123f;
1016 
1017             // clang-format off
1018             dprintf("     SHI %3u %4u %5u  %6.3f %6.3f %6.4f %7.5f  PROP-%s  %s\n",
1019                     qi, pi, hi, L.hit_q(hi), L.hit_phi(hi),
1020                     ddq, ddphi, prop_fail ? "FAIL" : "OK", dqdphi_presel ? "PASS" : "REJECT");
1021             // clang-format on
1022 
1023             if (prop_fail || !dqdphi_presel)
1024               continue;
1025             if (pqueue_size < NEW_MAX_HIT) {
1026               pqueue.push({new_ddphi, hi_orig});
1027               ++pqueue_size;
1028             } else if (new_ddphi < pqueue.top().score) {
1029               pqueue.pop();
1030               pqueue.push({new_ddphi, hi_orig});
1031             }
1032           }  //hi
1033         }    //pi
1034       }      //qi
1035 
1036       dprintf(" PQUEUE (%d)", pqueue_size);
1037       // Reverse hits so best dphis/scores come first in the hit-index list.
1038       m_XHitSize[itrack] = pqueue_size;
1039       while (pqueue_size) {
1040         --pqueue_size;
1041         m_XHitArr.At(itrack, pqueue_size, 0) = pqueue.top().hit_index;
1042         dprintf("   %d: %f %d", pqueue_size, pqueue.top().score, pqueue.top().hit_index);
1043         pqueue.pop();
1044       }
1045       dprintf("\n");
1046 
1047     }  //itrack
1048   }
1049 
1050   //==============================================================================
1051   // AddBestHit - Best Hit Track Finding
1052   //==============================================================================
1053 
1054   void MkFinder::addBestHit(const LayerOfHits &layer_of_hits, const int N_proc, const FindingFoos &fnd_foos) {
1055     // debug = true;
1056 
1057     MatriplexHitPacker mhp(layer_of_hits.hitArray());
1058 
1059     float minChi2[NN];
1060     int bestHit[NN];
1061     int maxSize = 0;
1062 
1063     // Determine maximum number of hits for tracks in the collection.
1064     for (int it = 0; it < NN; ++it) {
1065       if (it < N_proc) {
1066         if (m_XHitSize[it] > 0) {
1067           maxSize = std::max(maxSize, m_XHitSize[it]);
1068         }
1069       }
1070 
1071       bestHit[it] = -1;
1072       minChi2[it] = getHitSelDynamicChi2Cut(it, iP);
1073     }
1074 
1075     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
1076       //fixme what if size is zero???
1077 
1078       mhp.reset();
1079 
1080       //#pragma omp simd doesn't vectorize with current compilers
1081       for (int itrack = 0; itrack < NN; ++itrack) {
1082         if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1083           mhp.addInputAt(itrack, layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)));
1084         }
1085       }
1086 
1087       mhp.pack(m_msErr, m_msPar);
1088 
1089       //now compute the chi2 of track state vs hit
1090       MPlexQF outChi2;
1091       MPlexLV tmpPropPar;
1092       clearFailFlag();
1093       (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
1094                                      m_Par[iP],
1095                                      m_Chg,
1096                                      m_msErr,
1097                                      m_msPar,
1098                                      outChi2,
1099                                      tmpPropPar,
1100                                      m_FailFlag,
1101                                      N_proc,
1102                                      m_prop_config->finding_intra_layer_pflags,
1103                                      m_prop_config->finding_requires_propagation_to_hit_pos);
1104 
1105       //update best hit in case chi2<minChi2
1106 #pragma omp simd
1107       for (int itrack = 0; itrack < NN; ++itrack) {
1108         if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1109           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1110           dprint("chi2=" << chi2 << " minChi2[itrack]=" << minChi2[itrack]);
1111           if (chi2 < minChi2[itrack]) {
1112             minChi2[itrack] = chi2;
1113             bestHit[itrack] = m_XHitArr.At(itrack, hit_cnt, 0);
1114           }
1115         }
1116       }
1117     }  // end loop over hits
1118 
1119     //#pragma omp simd
1120     for (int itrack = 0; itrack < NN; ++itrack) {
1121       if (itrack >= N_proc) {
1122         continue;
1123       }
1124 
1125       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
1126         // Why am I doing this?
1127         m_msErr.setDiagonal3x3(itrack, 666);
1128         m_msPar(itrack, 0, 0) = m_Par[iP](itrack, 0, 0);
1129         m_msPar(itrack, 1, 0) = m_Par[iP](itrack, 1, 0);
1130         m_msPar(itrack, 2, 0) = m_Par[iP](itrack, 2, 0);
1131 
1132         // XXXX If not in gap, should get back the old track params. But they are gone ...
1133         // Would actually have to do it right after SelectHitIndices where updated params are still ok.
1134         // Here they got screwed during hit matching.
1135         // So, I'd store them there (into propagated params) and retrieve them here.
1136         // Or we decide not to care ...
1137 
1138         continue;
1139       }
1140 
1141       //fixme decide what to do in case no hit found
1142       if (bestHit[itrack] >= 0) {
1143         const Hit &hit = layer_of_hits.refHit(bestHit[itrack]);
1144         const float chi2 = minChi2[itrack];
1145 
1146         dprint("ADD BEST HIT FOR TRACK #"
1147                << itrack << std::endl
1148                << "prop x=" << m_Par[iP].constAt(itrack, 0, 0) << " y=" << m_Par[iP].constAt(itrack, 1, 0) << std::endl
1149                << "copy in hit #" << bestHit[itrack] << " x=" << hit.position()[0] << " y=" << hit.position()[1]);
1150 
1151         m_msErr.copyIn(itrack, hit.errArray());
1152         m_msPar.copyIn(itrack, hit.posArray());
1153         m_Chi2(itrack, 0, 0) += chi2;
1154 
1155         add_hit(itrack, bestHit[itrack], layer_of_hits.layer_id());
1156       } else {
1157         int fake_hit_idx = Hit::kHitMissIdx;
1158 
1159         if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
1160           // YYYYYY Config::store_missed_layers
1161           fake_hit_idx = Hit::kHitEdgeIdx;
1162         } else if (num_all_minus_one_hits(itrack)) {
1163           fake_hit_idx = Hit::kHitStopIdx;
1164         }
1165 
1166         dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx)
1167                                           << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0)));
1168 
1169         m_msErr.setDiagonal3x3(itrack, 666);
1170         m_msPar(itrack, 0, 0) = m_Par[iP](itrack, 0, 0);
1171         m_msPar(itrack, 1, 0) = m_Par[iP](itrack, 1, 0);
1172         m_msPar(itrack, 2, 0) = m_Par[iP](itrack, 2, 0);
1173         // Don't update chi2
1174 
1175         add_hit(itrack, fake_hit_idx, layer_of_hits.layer_id());
1176       }
1177     }
1178 
1179     // Update the track parameters with this hit. (Note that some calculations
1180     // are already done when computing chi2. Not sure it's worth caching them?)
1181 
1182     dprint("update parameters");
1183     clearFailFlag();
1184     (*fnd_foos.m_update_param_foo)(m_Err[iP],
1185                                    m_Par[iP],
1186                                    m_Chg,
1187                                    m_msErr,
1188                                    m_msPar,
1189                                    m_Err[iC],
1190                                    m_Par[iC],
1191                                    m_FailFlag,
1192                                    N_proc,
1193                                    m_prop_config->finding_intra_layer_pflags,
1194                                    m_prop_config->finding_requires_propagation_to_hit_pos);
1195 
1196     dprint("m_Par[iP](0,0,0)=" << m_Par[iP](0, 0, 0) << " m_Par[iC](0,0,0)=" << m_Par[iC](0, 0, 0));
1197   }
1198 
1199   //=======================================================
1200   // isStripQCompatible : check if prop is consistent with the barrel/endcap strip length
1201   //=======================================================
1202   bool isStripQCompatible(
1203       int itrack, bool isBarrel, const MPlexLS &pErr, const MPlexLV &pPar, const MPlexHS &msErr, const MPlexHV &msPar) {
1204     //check module compatibility via long strip side = L/sqrt(12)
1205     if (isBarrel) {  //check z direction only
1206       const float res = std::abs(msPar.constAt(itrack, 2, 0) - pPar.constAt(itrack, 2, 0));
1207       const float hitHL = sqrt(msErr.constAt(itrack, 2, 2) * 3.f);  //half-length
1208       const float qErr = sqrt(pErr.constAt(itrack, 2, 2));
1209       dprint("qCompat " << hitHL << " + " << 3.f * qErr << " vs " << res);
1210       return hitHL + std::max(3.f * qErr, 0.5f) > res;
1211     } else {  //project on xy, assuming the strip Length >> Width
1212       const float res[2]{msPar.constAt(itrack, 0, 0) - pPar.constAt(itrack, 0, 0),
1213                          msPar.constAt(itrack, 1, 0) - pPar.constAt(itrack, 1, 0)};
1214       const float hitT2 = msErr.constAt(itrack, 0, 0) + msErr.constAt(itrack, 1, 1);
1215       const float hitT2inv = 1.f / hitT2;
1216       const float proj[3] = {msErr.constAt(itrack, 0, 0) * hitT2inv,
1217                              msErr.constAt(itrack, 0, 1) * hitT2inv,
1218                              msErr.constAt(itrack, 1, 1) * hitT2inv};
1219       const float qErr =
1220           sqrt(std::abs(pErr.constAt(itrack, 0, 0) * proj[0] + 2.f * pErr.constAt(itrack, 0, 1) * proj[1] +
1221                         pErr.constAt(itrack, 1, 1) * proj[2]));  //take abs to avoid non-pos-def cases
1222       const float resProj =
1223           sqrt(res[0] * proj[0] * res[0] + 2.f * res[1] * proj[1] * res[0] + res[1] * proj[2] * res[1]);
1224       dprint("qCompat " << sqrt(hitT2 * 3.f) << " + " << 3.f * qErr << " vs " << resProj);
1225       return sqrt(hitT2 * 3.f) + std::max(3.f * qErr, 0.5f) > resProj;
1226     }
1227   }
1228 
1229   //=======================================================
1230   // passStripChargePCMfromTrack : apply the slope correction to charge per cm and cut using hit err matrix
1231   //         the raw pcm = charge/L_normal
1232   //         the corrected qCorr = charge/L_path = charge/(L_normal*p/p_zLocal) = pcm*p_zLocal/p
1233   //=======================================================
1234   bool passStripChargePCMfromTrack(
1235       int itrack, bool isBarrel, unsigned int pcm, unsigned int pcmMin, const MPlexLV &pPar, const MPlexHS &msErr) {
1236     //skip the overflow case
1237     if (pcm >= Hit::maxChargePerCM())
1238       return true;
1239 
1240     float qSF;
1241     if (isBarrel) {  //project in x,y, assuming zero-error direction is in this plane
1242       const float hitT2 = msErr.constAt(itrack, 0, 0) + msErr.constAt(itrack, 1, 1);
1243       const float hitT2inv = 1.f / hitT2;
1244       const float proj[3] = {msErr.constAt(itrack, 0, 0) * hitT2inv,
1245                              msErr.constAt(itrack, 0, 1) * hitT2inv,
1246                              msErr.constAt(itrack, 1, 1) * hitT2inv};
1247       const bool detXY_OK =
1248           std::abs(proj[0] * proj[2] - proj[1] * proj[1]) < 0.1f;  //check that zero-direction is close
1249       const float cosP = cos(pPar.constAt(itrack, 4, 0));
1250       const float sinP = sin(pPar.constAt(itrack, 4, 0));
1251       const float sinT = std::abs(sin(pPar.constAt(itrack, 5, 0)));
1252       //qSF = sqrt[(px,py)*(1-proj)*(px,py)]/p = sinT*sqrt[(cosP,sinP)*(1-proj)*(cosP,sinP)].
1253       qSF = detXY_OK ? sinT * std::sqrt(std::abs(1.f + cosP * cosP * proj[0] + sinP * sinP * proj[2] -
1254                                                  2.f * cosP * sinP * proj[1]))
1255                      : 1.f;
1256     } else {  //project on z
1257       // p_zLocal/p = p_z/p = cosT
1258       qSF = std::abs(cos(pPar.constAt(itrack, 5, 0)));
1259     }
1260 
1261     const float qCorr = pcm * qSF;
1262     dprint("pcm " << pcm << " * " << qSF << " = " << qCorr << " vs " << pcmMin);
1263     return qCorr > pcmMin;
1264   }
1265 
1266   //==============================================================================
1267   // FindCandidates - Standard Track Finding
1268   //==============================================================================
1269 
1270   void MkFinder::findCandidates(const LayerOfHits &layer_of_hits,
1271                                 std::vector<std::vector<TrackCand>> &tmp_candidates,
1272                                 const int offset,
1273                                 const int N_proc,
1274                                 const FindingFoos &fnd_foos) {
1275     // bool debug = true;
1276 
1277     MatriplexHitPacker mhp(layer_of_hits.hitArray());
1278 
1279     int maxSize = 0;
1280 
1281     // Determine maximum number of hits for tracks in the collection.
1282     for (int it = 0; it < NN; ++it) {
1283       if (it < N_proc) {
1284         if (m_XHitSize[it] > 0) {
1285           maxSize = std::max(maxSize, m_XHitSize[it]);
1286         }
1287       }
1288     }
1289 
1290     dprintf("FindCandidates max hits to process=%d\n", maxSize);
1291 
1292     int nHitsAdded[NN]{};
1293     bool isTooLargeCluster[NN]{false};
1294 
1295     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
1296       mhp.reset();
1297 
1298       int charge_pcm[NN];
1299 
1300       //#pragma omp simd doesn't vectorize with current compilers
1301       for (int itrack = 0; itrack < NN; ++itrack) {
1302         if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1303           const auto &hit = layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0));
1304           mhp.addInputAt(itrack, hit);
1305           charge_pcm[itrack] = hit.chargePerCM();
1306         }
1307       }
1308 
1309       mhp.pack(m_msErr, m_msPar);
1310 
1311       //now compute the chi2 of track state vs hit
1312       MPlexQF outChi2;
1313       MPlexLV propPar;
1314       clearFailFlag();
1315 
1316       if constexpr (Config::usePropToPlane) {
1317         // Maybe could use 2 matriplex packers ... ModuleInfo has 3 * SVector3 and uint
1318         MPlexHV norm, dir;
1319         packModuleNormDir(layer_of_hits, hit_cnt, norm, dir, N_proc);
1320         kalmanPropagateAndComputeChi2Plane(m_Err[iP],
1321                                            m_Par[iP],
1322                                            m_Chg,
1323                                            m_msErr,
1324                                            m_msPar,
1325                                            norm,
1326                                            dir,
1327                                            outChi2,
1328                                            propPar,
1329                                            m_FailFlag,
1330                                            N_proc,
1331                                            m_prop_config->finding_intra_layer_pflags,
1332                                            m_prop_config->finding_requires_propagation_to_hit_pos);
1333       } else {
1334         (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
1335                                        m_Par[iP],
1336                                        m_Chg,
1337                                        m_msErr,
1338                                        m_msPar,
1339                                        outChi2,
1340                                        propPar,
1341                                        m_FailFlag,
1342                                        N_proc,
1343                                        m_prop_config->finding_intra_layer_pflags,
1344                                        m_prop_config->finding_requires_propagation_to_hit_pos);
1345       }
1346 
1347       // Now update the track parameters with this hit (note that some
1348       // calculations are already done when computing chi2, to be optimized).
1349       // 1. This is not needed for candidates the hit is not added to, but it's
1350       // vectorized so doing it serially below should take the same time.
1351       // 2. Still it's a waste of time in case the hit is not added to any of the
1352       // candidates, so check beforehand that at least one cand needs update.
1353       bool oneCandPassCut = false;
1354       for (int itrack = 0; itrack < NN; ++itrack) {
1355         if (itrack >= N_proc) {
1356           continue;
1357         }
1358         float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
1359 
1360         if (hit_cnt < m_XHitSize[itrack]) {
1361           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1362           dprint("chi2=" << chi2);
1363           if (chi2 < max_c2) {
1364             bool isCompatible = true;
1365             if (!layer_of_hits.is_pixel()) {
1366               //check module compatibility via long strip side = L/sqrt(12)
1367               isCompatible =
1368                   isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
1369 
1370               //rescale strip charge to track parameters and reapply the cut
1371               isCompatible &= passStripChargePCMfromTrack(
1372                   itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
1373             }
1374             // Select only SiStrip hits with cluster size < maxClusterSize
1375             if (!layer_of_hits.is_pixel()) {
1376               if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
1377                   m_iteration_params->maxClusterSize) {
1378                 isTooLargeCluster[itrack] = true;
1379                 isCompatible = false;
1380               }
1381             }
1382 
1383             if (isCompatible) {
1384               oneCandPassCut = true;
1385               break;
1386             }
1387           }
1388         }
1389       }
1390 
1391       if (oneCandPassCut) {
1392         MPlexQI tmpChg = m_Chg;
1393         clearFailFlag();
1394         (*fnd_foos.m_update_param_foo)(m_Err[iP],
1395                                        m_Par[iP],
1396                                        tmpChg,
1397                                        m_msErr,
1398                                        m_msPar,
1399                                        m_Err[iC],
1400                                        m_Par[iC],
1401                                        m_FailFlag,
1402                                        N_proc,
1403                                        m_prop_config->finding_intra_layer_pflags,
1404                                        m_prop_config->finding_requires_propagation_to_hit_pos);
1405 
1406         dprint("update parameters" << std::endl
1407                                    << "propagated track parameters x=" << m_Par[iP].constAt(0, 0, 0)
1408                                    << " y=" << m_Par[iP].constAt(0, 1, 0) << std::endl
1409                                    << "               hit position x=" << m_msPar.constAt(0, 0, 0)
1410                                    << " y=" << m_msPar.constAt(0, 1, 0) << std::endl
1411                                    << "   updated track parameters x=" << m_Par[iC].constAt(0, 0, 0)
1412                                    << " y=" << m_Par[iC].constAt(0, 1, 0));
1413 
1414         //create candidate with hit in case chi2 < max_c2
1415         //fixme: please vectorize me... (not sure it's possible in this case)
1416         for (int itrack = 0; itrack < NN; ++itrack) {
1417           if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1418             const float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
1419             const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1420             dprint("chi2=" << chi2);
1421             if (chi2 < max_c2) {
1422               bool isCompatible = true;
1423               if (!layer_of_hits.is_pixel()) {
1424                 //check module compatibility via long strip side = L/sqrt(12)
1425                 isCompatible =
1426                     isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
1427 
1428                 //rescale strip charge to track parameters and reapply the cut
1429                 isCompatible &= passStripChargePCMfromTrack(
1430                     itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
1431               }
1432               // Select only SiStrip hits with cluster size < maxClusterSize
1433               if (!layer_of_hits.is_pixel()) {
1434                 if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
1435                     m_iteration_params->maxClusterSize)
1436                   isCompatible = false;
1437               }
1438 
1439               if (isCompatible) {
1440                 bool hitExists = false;
1441                 int maxHits = m_NFoundHits(itrack, 0, 0);
1442                 if (layer_of_hits.is_pixel()) {
1443                   for (int i = 0; i <= maxHits; ++i) {
1444                     if (i > 2)
1445                       break;
1446                     if (m_HoTArrs[itrack][i].layer == layer_of_hits.layer_id()) {
1447                       hitExists = true;
1448                       break;
1449                     }
1450                   }
1451                 }
1452                 if (hitExists)
1453                   continue;
1454 
1455                 nHitsAdded[itrack]++;
1456                 dprint("chi2 cut passed, creating new candidate");
1457                 // Create a new candidate and fill the tmp_candidates output vector.
1458                 // QQQ only instantiate if it will pass, be better than N_best
1459 
1460                 const int hit_idx = m_XHitArr.At(itrack, hit_cnt, 0);
1461 
1462                 TrackCand newcand;
1463                 copy_out(newcand, itrack, iC);
1464                 newcand.setCharge(tmpChg(itrack, 0, 0));
1465                 newcand.addHitIdx(hit_idx, layer_of_hits.layer_id(), chi2);
1466                 newcand.setScore(getScoreCand(m_steering_params->m_track_scorer,
1467                                               newcand,
1468                                               true /*penalizeTailMissHits*/,
1469                                               true /*inFindCandidates*/));
1470                 newcand.setOriginIndex(m_CandIdx(itrack, 0, 0));
1471 
1472                 // To apply a fixed cut instead of dynamic cut for overlap: m_iteration_params->chi2CutOverlap
1473                 if (chi2 < max_c2) {
1474                   CombCandidate &ccand = *newcand.combCandidate();
1475                   ccand[m_CandIdx(itrack, 0, 0)].considerHitForOverlap(
1476                       hit_idx, layer_of_hits.refHit(hit_idx).detIDinLayer(), chi2);
1477                 }
1478 
1479                 dprint("updated track parameters x=" << newcand.parameters()[0] << " y=" << newcand.parameters()[1]
1480                                                      << " z=" << newcand.parameters()[2]
1481                                                      << " pt=" << 1. / newcand.parameters()[3]);
1482 
1483                 tmp_candidates[m_SeedIdx(itrack, 0, 0) - offset].emplace_back(newcand);
1484               }
1485             }
1486           }
1487         }
1488       }  //end if (oneCandPassCut)
1489 
1490     }  //end loop over hits
1491 
1492     //now add invalid hit
1493     //fixme: please vectorize me...
1494     for (int itrack = 0; itrack < NN; ++itrack) {
1495       // Cands that miss the layer are stashed away in MkBuilder(), before propagation,
1496       // and then merged back afterwards.
1497       if (itrack >= N_proc || m_XWsrResult[itrack].m_wsr == WSR_Outside) {
1498         continue;
1499       }
1500 
1501       int fake_hit_idx = ((num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand) &&
1502                           (m_NTailMinusOneHits(itrack, 0, 0) < m_iteration_params->maxConsecHoles))
1503                              ? Hit::kHitMissIdx
1504                              : Hit::kHitStopIdx;
1505 
1506       if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
1507         // YYYYYY m_iteration_params->store_missed_layers
1508         fake_hit_idx = Hit::kHitEdgeIdx;
1509       }
1510       //now add fake hit for tracks that passsed through inactive modules
1511       else if (m_XWsrResult[itrack].m_in_gap == true && nHitsAdded[itrack] == 0) {
1512         fake_hit_idx = Hit::kHitInGapIdx;
1513       }
1514       //now add fake hit for cases where hit cluster size is larger than maxClusterSize
1515       else if (isTooLargeCluster[itrack] == true && nHitsAdded[itrack] == 0) {
1516         fake_hit_idx = Hit::kHitMaxClusterIdx;
1517       }
1518 
1519       dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx)
1520                                         << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0)));
1521 
1522       // QQQ as above, only create and add if score better
1523       TrackCand newcand;
1524       copy_out(newcand, itrack, iP);
1525       newcand.addHitIdx(fake_hit_idx, layer_of_hits.layer_id(), 0.);
1526       newcand.setScore(getScoreCand(
1527           m_steering_params->m_track_scorer, newcand, true /*penalizeTailMissHits*/, true /*inFindCandidates*/));
1528       // Only relevant when we actually add a hit
1529       // newcand.setOriginIndex(m_CandIdx(itrack, 0, 0));
1530       tmp_candidates[m_SeedIdx(itrack, 0, 0) - offset].emplace_back(newcand);
1531     }
1532   }
1533 
1534   //==============================================================================
1535   // FindCandidatesCloneEngine - Clone Engine Track Finding
1536   //==============================================================================
1537 
1538   void MkFinder::findCandidatesCloneEngine(const LayerOfHits &layer_of_hits,
1539                                            CandCloner &cloner,
1540                                            const int offset,
1541                                            const int N_proc,
1542                                            const FindingFoos &fnd_foos) {
1543     // bool debug = true;
1544 
1545     MatriplexHitPacker mhp(layer_of_hits.hitArray());
1546 
1547     int maxSize = 0;
1548 
1549     // Determine maximum number of hits for tracks in the collection.
1550 #pragma omp simd
1551     for (int it = 0; it < NN; ++it) {
1552       if (it < N_proc) {
1553         if (m_XHitSize[it] > 0) {
1554           maxSize = std::max(maxSize, m_XHitSize[it]);
1555         }
1556       }
1557     }
1558 
1559     dprintf("FindCandidatesCloneEngine max hits to process=%d\n", maxSize);
1560 
1561     int nHitsAdded[NN]{};
1562     bool isTooLargeCluster[NN]{false};
1563 
1564     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
1565       mhp.reset();
1566 
1567       int charge_pcm[NN];
1568 
1569       //#pragma omp simd doesn't vectorize with current compilers
1570       for (int itrack = 0; itrack < NN; ++itrack) {
1571         if (itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1572           const auto &hit = layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0));
1573           mhp.addInputAt(itrack, hit);
1574           charge_pcm[itrack] = hit.chargePerCM();
1575         }
1576       }
1577 
1578       mhp.pack(m_msErr, m_msPar);
1579 
1580       //now compute the chi2 of track state vs hit
1581       MPlexQF outChi2;
1582       MPlexLV propPar;
1583       clearFailFlag();
1584 
1585       if constexpr (Config::usePropToPlane) {
1586         // Maybe could use 2 matriplex packers ... ModuleInfo has 3 * SVector3 and uint
1587         MPlexHV norm, dir;
1588         packModuleNormDir(layer_of_hits, hit_cnt, norm, dir, N_proc);
1589         kalmanPropagateAndComputeChi2Plane(m_Err[iP],
1590                                            m_Par[iP],
1591                                            m_Chg,
1592                                            m_msErr,
1593                                            m_msPar,
1594                                            norm,
1595                                            dir,
1596                                            outChi2,
1597                                            propPar,
1598                                            m_FailFlag,
1599                                            N_proc,
1600                                            m_prop_config->finding_intra_layer_pflags,
1601                                            m_prop_config->finding_requires_propagation_to_hit_pos);
1602       } else {
1603         (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
1604                                        m_Par[iP],
1605                                        m_Chg,
1606                                        m_msErr,
1607                                        m_msPar,
1608                                        outChi2,
1609                                        propPar,
1610                                        m_FailFlag,
1611                                        N_proc,
1612                                        m_prop_config->finding_intra_layer_pflags,
1613                                        m_prop_config->finding_requires_propagation_to_hit_pos);
1614       }
1615 
1616       //#pragma omp simd  // DOES NOT VECTORIZE AS IT IS NOW
1617       for (int itrack = 0; itrack < NN; ++itrack) {
1618         // We can be in failed state from the initial propagation before selectHitIndices
1619         // and there hit_count for track is set to -1 and WSR state to Failed, handled below.
1620         // Or we might have hit it here in propagate-to-hit.
1621         // PROP-FAIL-ENABLE FailFlag check to be enabled when propagation failure
1622         // detection is properly implemented in propagate-to-R/Z.
1623         if (/*!m_FailFlag[itrack] &&*/ itrack < N_proc && hit_cnt < m_XHitSize[itrack]) {
1624           // make sure the hit was in the compatiblity window for the candidate
1625           const float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
1626           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1627           // XXX-NUM-ERR assert(chi2 >= 0);
1628 
1629           dprint("chi2=" << chi2 << " for trkIdx=" << itrack << " hitIdx=" << m_XHitArr.At(itrack, hit_cnt, 0));
1630           if (chi2 < max_c2) {
1631             bool isCompatible = true;
1632             if (!layer_of_hits.is_pixel()) {
1633               //check module compatibility via long strip side = L/sqrt(12)
1634               isCompatible =
1635                   isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
1636 
1637               //rescale strip charge to track parameters and reapply the cut
1638               isCompatible &= passStripChargePCMfromTrack(
1639                   itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
1640             }
1641 
1642             // Select only SiStrip hits with cluster size < maxClusterSize
1643             if (!layer_of_hits.is_pixel()) {
1644               if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
1645                   m_iteration_params->maxClusterSize) {
1646                 isTooLargeCluster[itrack] = true;
1647                 isCompatible = false;
1648               }
1649             }
1650 
1651             if (isCompatible) {
1652               CombCandidate &ccand = cloner.combCandWithOriginalIndex(m_SeedIdx(itrack, 0, 0));
1653               bool hitExists = false;
1654               int maxHits = m_NFoundHits(itrack, 0, 0);
1655               if (layer_of_hits.is_pixel()) {
1656                 for (int i = 0; i <= maxHits; ++i) {
1657                   if (i > 2)
1658                     break;
1659                   if (ccand.hot(i).layer == layer_of_hits.layer_id()) {
1660                     hitExists = true;
1661                     break;
1662                   }
1663                 }
1664               }
1665               if (hitExists)
1666                 continue;
1667 
1668               nHitsAdded[itrack]++;
1669               const int hit_idx = m_XHitArr.At(itrack, hit_cnt, 0);
1670 
1671               // Register hit for overlap consideration, if chi2 cut is passed
1672               // To apply a fixed cut instead of dynamic cut for overlap: m_iteration_params->chi2CutOverlap
1673               if (chi2 < max_c2) {
1674                 ccand[m_CandIdx(itrack, 0, 0)].considerHitForOverlap(
1675                     hit_idx, layer_of_hits.refHit(hit_idx).detIDinLayer(), chi2);
1676               }
1677 
1678               IdxChi2List tmpList;
1679               tmpList.trkIdx = m_CandIdx(itrack, 0, 0);
1680               tmpList.hitIdx = hit_idx;
1681               tmpList.module = layer_of_hits.refHit(hit_idx).detIDinLayer();
1682               tmpList.nhits = m_NFoundHits(itrack, 0, 0) + 1;
1683               tmpList.ntailholes = 0;
1684               tmpList.noverlaps = m_NOverlapHits(itrack, 0, 0);
1685               tmpList.nholes = num_all_minus_one_hits(itrack);
1686               tmpList.pt = std::abs(1.0f / m_Par[iP].At(itrack, 3, 0));
1687               tmpList.chi2 = m_Chi2(itrack, 0, 0) + chi2;
1688               tmpList.chi2_hit = chi2;
1689               tmpList.score = getScoreStruct(m_steering_params->m_track_scorer, tmpList);
1690               cloner.add_cand(m_SeedIdx(itrack, 0, 0) - offset, tmpList);
1691 
1692               dprint("  adding hit with hit_cnt=" << hit_cnt << " for trkIdx=" << tmpList.trkIdx
1693                                                   << " orig Seed=" << m_Label(itrack, 0, 0));
1694             }
1695           }
1696         }
1697       }
1698 
1699     }  //end loop over hits
1700 
1701     //now add invalid hit
1702     for (int itrack = 0; itrack < NN; ++itrack) {
1703       dprint("num_all_minus_one_hits(" << itrack << ")=" << num_all_minus_one_hits(itrack));
1704 
1705       // Cands that miss the layer are stashed away in MkBuilder(), before propagation,
1706       // and then merged back afterwards.
1707       if (itrack >= N_proc || m_XWsrResult[itrack].m_wsr == WSR_Outside) {
1708         continue;
1709       }
1710 
1711       // int fake_hit_idx = num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand ? -1 : -2;
1712       int fake_hit_idx = ((num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand) &&
1713                           (m_NTailMinusOneHits(itrack, 0, 0) < m_iteration_params->maxConsecHoles))
1714                              ? Hit::kHitMissIdx
1715                              : Hit::kHitStopIdx;
1716 
1717       if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
1718         fake_hit_idx = Hit::kHitEdgeIdx;
1719       }
1720       //now add fake hit for tracks that passsed through inactive modules
1721       else if (m_XWsrResult[itrack].m_in_gap == true && nHitsAdded[itrack] == 0) {
1722         fake_hit_idx = Hit::kHitInGapIdx;
1723       }
1724       //now add fake hit for cases where hit cluster size is larger than maxClusterSize
1725       else if (isTooLargeCluster[itrack] == true && nHitsAdded[itrack] == 0) {
1726         fake_hit_idx = Hit::kHitMaxClusterIdx;
1727       }
1728 
1729       // PROP-FAIL-ENABLE The following to be enabled when propagation failure
1730       // detection is properly implemented in propagate-to-R/Z.
1731       // // Override for failed propagation, this trumps all other cases.
1732       // if (m_XWsrResult[itrack].m_wsr == WSR_Failed) {
1733       //   fake_hit_idx = Hit::kHitStopIdx;
1734       // }
1735 
1736       IdxChi2List tmpList;
1737       tmpList.trkIdx = m_CandIdx(itrack, 0, 0);
1738       tmpList.hitIdx = fake_hit_idx;
1739       tmpList.module = -1;
1740       tmpList.nhits = m_NFoundHits(itrack, 0, 0);
1741       tmpList.ntailholes = (fake_hit_idx == Hit::kHitMissIdx ? m_NTailMinusOneHits(itrack, 0, 0) + 1
1742                                                              : m_NTailMinusOneHits(itrack, 0, 0));
1743       tmpList.noverlaps = m_NOverlapHits(itrack, 0, 0);
1744       tmpList.nholes = num_inside_minus_one_hits(itrack);
1745       tmpList.pt = std::abs(1.0f / m_Par[iP].At(itrack, 3, 0));
1746       tmpList.chi2 = m_Chi2(itrack, 0, 0);
1747       tmpList.chi2_hit = 0;
1748       tmpList.score = getScoreStruct(m_steering_params->m_track_scorer, tmpList);
1749       cloner.add_cand(m_SeedIdx(itrack, 0, 0) - offset, tmpList);
1750       dprint("adding invalid hit " << fake_hit_idx);
1751     }
1752   }
1753 
1754   //==============================================================================
1755   // UpdateWithLoadedHit
1756   //==============================================================================
1757 
1758   void MkFinder::updateWithLoadedHit(int N_proc, const LayerOfHits &layer_of_hits, const FindingFoos &fnd_foos) {
1759     // See comment in MkBuilder::find_tracks_in_layer() about intra / inter flags used here
1760     // for propagation to the hit.
1761     clearFailFlag();
1762     if constexpr (Config::usePropToPlane) {
1763       MPlexHV norm, dir;
1764       packModuleNormDir(layer_of_hits, 0, norm, dir, N_proc);
1765       kalmanPropagateAndUpdatePlane(m_Err[iP],
1766                                     m_Par[iP],
1767                                     m_Chg,
1768                                     m_msErr,
1769                                     m_msPar,
1770                                     norm,
1771                                     dir,
1772                                     m_Err[iC],
1773                                     m_Par[iC],
1774                                     m_FailFlag,
1775                                     N_proc,
1776                                     m_prop_config->finding_inter_layer_pflags,
1777                                     m_prop_config->finding_requires_propagation_to_hit_pos);
1778     } else {
1779       (*fnd_foos.m_update_param_foo)(m_Err[iP],
1780                                      m_Par[iP],
1781                                      m_Chg,
1782                                      m_msErr,
1783                                      m_msPar,
1784                                      m_Err[iC],
1785                                      m_Par[iC],
1786                                      m_FailFlag,
1787                                      N_proc,
1788                                      m_prop_config->finding_inter_layer_pflags,
1789                                      m_prop_config->finding_requires_propagation_to_hit_pos);
1790     }
1791 
1792     // PROP-FAIL-ENABLE The following to be enabled when propagation failure
1793     // detection is properly implemented in propagate-to-R/Z.
1794     // for (int i = 0; i < N_proc; ++i) {
1795     //   if (m_FailFlag[i]) {
1796     //     dprintf("MkFinder::updateWithLoadedHit fail in update, recovering.\n");
1797     //     m_Err[iC].copySlot(i, m_Err[iP]);
1798     //     m_Par[iC].copySlot(i, m_Par[iP]);
1799     //   }
1800     // }
1801   }
1802 
1803   void MkFinder::chi2OfLoadedHit(int N_proc, const FindingFoos &fnd_foos) {
1804     // We expect input in iC slots from above function.
1805     // See comment in MkBuilder::find_tracks_in_layer() about intra / inter flags used here
1806     // for propagation to the hit.
1807     clearFailFlag();
1808     (*fnd_foos.m_compute_chi2_foo)(m_Err[iC],
1809                                    m_Par[iC],
1810                                    m_Chg,
1811                                    m_msErr,
1812                                    m_msPar,
1813                                    m_Chi2,
1814                                    m_Par[iP],
1815                                    m_FailFlag,
1816                                    N_proc,
1817                                    m_prop_config->finding_inter_layer_pflags,
1818                                    m_prop_config->finding_requires_propagation_to_hit_pos);
1819 
1820     // PROP-FAIL-ENABLE .... removed here
1821   }
1822 
1823   //==============================================================================
1824   // CopyOutParErr
1825   //==============================================================================
1826 
1827   void MkFinder::copyOutParErr(std::vector<CombCandidate> &seed_cand_vec, int N_proc, bool outputProp) const {
1828     const int iO = outputProp ? iP : iC;
1829 
1830     for (int i = 0; i < NN; ++i) {
1831       if (i < N_proc) {
1832         TrackCand &cand = seed_cand_vec[m_SeedIdx(i, 0, 0)][m_CandIdx(i, 0, 0)];
1833 
1834         // Set the track state to the updated parameters
1835         m_Err[iO].copyOut(i, cand.errors_nc().Array());
1836         m_Par[iO].copyOut(i, cand.parameters_nc().Array());
1837         cand.setCharge(m_Chg(i, 0, 0));
1838 
1839         dprint((outputProp ? "propagated" : "updated")
1840                << " track parameters x=" << cand.parameters()[0] << " y=" << cand.parameters()[1]
1841                << " z=" << cand.parameters()[2] << " pt=" << 1. / cand.parameters()[3] << " posEta=" << cand.posEta());
1842       }
1843     }
1844   }
1845 
1846   //==============================================================================
1847   // Backward Fit hack
1848   //==============================================================================
1849 
1850   void MkFinder::bkFitInputTracks(TrackVec &cands, int beg, int end) {
1851     // Uses HitOnTrack vector from Track directly + a local cursor array to current hit.
1852 
1853     MatriplexTrackPacker mtp(&cands[beg]);
1854 
1855     int itrack = 0;
1856 
1857     for (int i = beg; i < end; ++i, ++itrack) {
1858       const Track &trk = cands[i];
1859 
1860       m_Chg(itrack, 0, 0) = trk.charge();
1861       m_CurHit[itrack] = trk.nTotalHits() - 1;
1862       m_HoTArr[itrack] = trk.getHitsOnTrackArray();
1863 
1864       mtp.addInput(trk);
1865     }
1866 
1867     m_Chi2.setVal(0);
1868 
1869     mtp.pack(m_Err[iC], m_Par[iC]);
1870 
1871     m_Err[iC].scale(100.0f);
1872   }
1873 
1874   void MkFinder::bkFitInputTracks(EventOfCombCandidates &eocss, int beg, int end) {
1875     // Could as well use HotArrays from tracks directly + a local cursor array to last hit.
1876 
1877     // XXXX - shall we assume only TrackCand-zero is needed and that we can freely
1878     // bork the HoTNode array?
1879 
1880     MatriplexTrackPacker mtp(&eocss[beg][0]);
1881 
1882     int itrack = 0;
1883 
1884     for (int i = beg; i < end; ++i, ++itrack) {
1885       const TrackCand &trk = eocss[i][0];
1886 
1887       m_Chg(itrack, 0, 0) = trk.charge();
1888       m_CurNode[itrack] = trk.lastCcIndex();
1889       m_HoTNodeArr[itrack] = trk.combCandidate()->hotsData();
1890 
1891       // XXXX Need TrackCand* to update num-hits. Unless I collect info elsewhere
1892       // and fix it in BkFitOutputTracks.
1893       m_TrkCand[itrack] = &eocss[i][0];
1894 
1895       mtp.addInput(trk);
1896     }
1897 
1898     m_Chi2.setVal(0);
1899 
1900     mtp.pack(m_Err[iC], m_Par[iC]);
1901 
1902     m_Err[iC].scale(100.0f);
1903   }
1904 
1905   //------------------------------------------------------------------------------
1906 
1907   void MkFinder::bkFitOutputTracks(TrackVec &cands, int beg, int end, bool outputProp) {
1908     // Only copy out track params / errors / chi2, all the rest is ok.
1909 
1910     const int iO = outputProp ? iP : iC;
1911 
1912     int itrack = 0;
1913     for (int i = beg; i < end; ++i, ++itrack) {
1914       Track &trk = cands[i];
1915 
1916       m_Err[iO].copyOut(itrack, trk.errors_nc().Array());
1917       m_Par[iO].copyOut(itrack, trk.parameters_nc().Array());
1918 
1919       trk.setChi2(m_Chi2(itrack, 0, 0));
1920       if (isFinite(trk.chi2())) {
1921         trk.setScore(getScoreCand(m_steering_params->m_track_scorer, trk));
1922       }
1923     }
1924   }
1925 
1926   void MkFinder::bkFitOutputTracks(EventOfCombCandidates &eocss, int beg, int end, bool outputProp) {
1927     // Only copy out track params / errors / chi2, all the rest is ok.
1928 
1929     // XXXX - where will rejected hits get removed?
1930 
1931     const int iO = outputProp ? iP : iC;
1932 
1933     int itrack = 0;
1934     for (int i = beg; i < end; ++i, ++itrack) {
1935       TrackCand &trk = eocss[i][0];
1936 
1937       m_Err[iO].copyOut(itrack, trk.errors_nc().Array());
1938       m_Par[iO].copyOut(itrack, trk.parameters_nc().Array());
1939 
1940       trk.setChi2(m_Chi2(itrack, 0, 0));
1941       if (isFinite(trk.chi2())) {
1942         trk.setScore(getScoreCand(m_steering_params->m_track_scorer, trk));
1943       }
1944     }
1945   }
1946 
1947   //------------------------------------------------------------------------------
1948 
1949 #if defined(DEBUG_BACKWARD_FIT) || defined(DEBUG_BACKWARD_FIT_BH)
1950   namespace {
1951     float e2s(float x) { return 1e4 * std::sqrt(x); }
1952   }  // namespace
1953 #endif
1954 
1955   void MkFinder::bkFitFitTracksBH(const EventOfHits &eventofhits,
1956                                   const SteeringParams &st_par,
1957                                   const int N_proc,
1958                                   bool chiDebug) {
1959     // Prototyping final backward fit.
1960     // This works with track-finding indices, before remapping.
1961     //
1962     // Layers should be collected during track finding and list all layers that have actual hits.
1963     // Then we could avoid checking which layers actually do have hits.
1964 
1965     MPlexQF tmp_chi2;
1966     float tmp_err[6] = {666, 0, 666, 0, 0, 666};
1967     float tmp_pos[3];
1968 
1969     for (auto lp_iter = st_par.m_layer_plan.rbegin(); lp_iter != st_par.m_layer_plan.rend(); ++lp_iter) {
1970       const int layer = lp_iter->m_layer;
1971 
1972       const LayerOfHits &L = eventofhits[layer];
1973       const LayerInfo &LI = *L.layer_info();
1974 
1975       int count = 0;
1976       for (int i = 0; i < N_proc; ++i) {
1977         while (m_CurHit[i] >= 0 && m_HoTArr[i][m_CurHit[i]].index < 0)
1978           --m_CurHit[i];
1979 
1980         if (m_CurHit[i] >= 0 && m_HoTArr[i][m_CurHit[i]].layer == layer) {
1981           // Skip the overlap hits -- if they exist.
1982           // 1. Overlap hit gets placed *after* the original hit in TrackCand::exportTrack()
1983           // which is *before* in the reverse iteration that we are doing here.
1984           // 2. Seed-hit merging can result in more than two hits per layer.
1985           while (m_CurHit[i] > 0 && m_HoTArr[i][m_CurHit[i] - 1].layer == layer)
1986             --m_CurHit[i];
1987 
1988           const Hit &hit = L.refHit(m_HoTArr[i][m_CurHit[i]].index);
1989           m_msErr.copyIn(i, hit.errArray());
1990           m_msPar.copyIn(i, hit.posArray());
1991           ++count;
1992           --m_CurHit[i];
1993         } else {
1994           tmp_pos[0] = m_Par[iC](i, 0, 0);
1995           tmp_pos[1] = m_Par[iC](i, 1, 0);
1996           tmp_pos[2] = m_Par[iC](i, 2, 0);
1997           m_msErr.copyIn(i, tmp_err);
1998           m_msPar.copyIn(i, tmp_pos);
1999         }
2000       }
2001 
2002       if (count == 0)
2003         continue;
2004 
2005       // ZZZ Could add missing hits here, only if there are any actual matches.
2006 
2007       if (LI.is_barrel()) {
2008         propagateTracksToHitR(m_msPar, N_proc, m_prop_config->backward_fit_pflags);
2009 
2010         kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
2011                         m_Err[iP],
2012                         m_Par[iP],
2013                         m_msErr,
2014                         m_msPar,
2015                         m_Err[iC],
2016                         m_Par[iC],
2017                         tmp_chi2,
2018                         N_proc);
2019       } else {
2020         propagateTracksToHitZ(m_msPar, N_proc, m_prop_config->backward_fit_pflags);
2021 
2022         kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params,
2023                               m_Err[iP],
2024                               m_Par[iP],
2025                               m_msErr,
2026                               m_msPar,
2027                               m_Err[iC],
2028                               m_Par[iC],
2029                               tmp_chi2,
2030                               N_proc);
2031       }
2032 
2033       //fixup invpt sign and charge
2034       for (int n = 0; n < NN; ++n) {
2035         if (n < N_proc && m_Par[iC].At(n, 3, 0) < 0) {
2036           m_Chg.At(n, 0, 0) = -m_Chg.At(n, 0, 0);
2037           m_Par[iC].At(n, 3, 0) = -m_Par[iC].At(n, 3, 0);
2038         }
2039       }
2040 
2041 #ifdef DEBUG_BACKWARD_FIT_BH
2042       // Dump per hit chi2
2043       for (int i = 0; i < N_proc; ++i) {
2044         float r_h = std::hypot(m_msPar.At(i, 0, 0), m_msPar.At(i, 1, 0));
2045         float r_t = std::hypot(m_Par[iC].At(i, 0, 0), m_Par[iC].At(i, 1, 0));
2046 
2047         // if ((std::isnan(tmp_chi2[i]) || std::isnan(r_t)))
2048         // if ( ! std::isnan(tmp_chi2[i]) && tmp_chi2[i] > 0) // && tmp_chi2[i] > 30)
2049         if (chiDebug) {
2050           int ti = iP;
2051           printf(
2052               "CHIHIT %3d %10g %10g %10g %10g %10g %11.5g %11.5g %11.5g %10g %10g %10g %10g %11.5g %11.5g %11.5g %10g "
2053               "%10g %10g %10g %10g %11.5g %11.5g\n",
2054               layer,
2055               tmp_chi2[i],
2056               m_msPar.At(i, 0, 0),
2057               m_msPar.At(i, 1, 0),
2058               m_msPar.At(i, 2, 0),
2059               r_h,  // x_h y_h z_h r_h -- hit pos
2060               e2s(m_msErr.At(i, 0, 0)),
2061               e2s(m_msErr.At(i, 1, 1)),
2062               e2s(m_msErr.At(i, 2, 2)),  // ex_h ey_h ez_h -- hit errors
2063               m_Par[ti].At(i, 0, 0),
2064               m_Par[ti].At(i, 1, 0),
2065               m_Par[ti].At(i, 2, 0),
2066               r_t,  // x_t y_t z_t r_t -- track pos
2067               e2s(m_Err[ti].At(i, 0, 0)),
2068               e2s(m_Err[ti].At(i, 1, 1)),
2069               e2s(m_Err[ti].At(i, 2, 2)),  // ex_t ey_t ez_t -- track errors
2070               1.0f / m_Par[ti].At(i, 3, 0),
2071               m_Par[ti].At(i, 4, 0),
2072               m_Par[ti].At(i, 5, 0),                                     // pt, phi, theta
2073               std::atan2(m_msPar.At(i, 1, 0), m_msPar.At(i, 0, 0)),      // phi_h
2074               std::atan2(m_Par[ti].At(i, 1, 0), m_Par[ti].At(i, 0, 0)),  // phi_t
2075               1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0),
2076                                 m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)),  // d_xy
2077               1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0))             // d_z
2078               // e2s((m_msErr.At(i,0,0) + m_msErr.At(i,1,1)) / (r_h * r_h)),     // ephi_h
2079               // e2s((m_Err[ti].At(i,0,0) + m_Err[ti].At(i,1,1)) / (r_t * r_t))  // ephi_t
2080           );
2081         }
2082       }
2083 #endif
2084 
2085       // update chi2
2086       m_Chi2.add(tmp_chi2);
2087     }
2088   }
2089 
2090   //------------------------------------------------------------------------------
2091 
2092   void MkFinder::print_par_err(int corp, int mslot) const {
2093 #ifdef DEBUG
2094     printf("Parameters:\n");
2095     for (int i = 0; i < 6; ++i) {
2096       printf("  %12.4g", m_Par[corp].constAt(mslot, i, 0));
2097     }
2098     printf("\nError matrix\n");
2099     for (int i = 0; i < 6; ++i) {
2100       for (int j = 0; j < 6; ++j) {
2101         printf("  %12.4g", m_Err[corp].constAt(mslot, i, j));
2102       }
2103       printf("\n");
2104     }
2105 #endif
2106   }
2107 
2108   void MkFinder::bkFitFitTracks(const EventOfHits &eventofhits,
2109                                 const SteeringParams &st_par,
2110                                 const int N_proc,
2111                                 bool chiDebug) {
2112     // Prototyping final backward fit.
2113     // This works with track-finding indices, before remapping.
2114     //
2115     // Layers should be collected during track finding and list all layers that have actual hits.
2116     // Then we could avoid checking which layers actually do have hits.
2117 
2118     // bool debug = true;
2119 
2120     MPlexQF tmp_chi2;
2121     MPlexQI no_mat_effs;
2122     float tmp_err[6] = {666, 0, 666, 0, 0, 666};
2123     float tmp_pos[3];
2124 
2125 #if defined(DEBUG_PROP_UPDATE)
2126     const int DSLOT = 0;
2127     printf("bkfit entry, track in slot %d\n", DSLOT);
2128     print_par_err(iC, DSLOT);
2129 #endif
2130 
2131     for (auto lp_iter = st_par.make_iterator(SteeringParams::IT_BkwFit); lp_iter.is_valid(); ++lp_iter) {
2132       const int layer = lp_iter.layer();
2133 
2134       const LayerOfHits &L = eventofhits[layer];
2135       const LayerInfo &LI = *L.layer_info();
2136 
2137 #if defined(DEBUG_BACKWARD_FIT)
2138       const Hit *last_hit_ptr[NN];
2139 #endif
2140 
2141       no_mat_effs.setVal(0);
2142       int done_count = 0;
2143       int here_count = 0;
2144       for (int i = 0; i < N_proc; ++i) {
2145         while (m_CurNode[i] >= 0 && m_HoTNodeArr[i][m_CurNode[i]].m_hot.index < 0) {
2146           m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
2147         }
2148 
2149         if (m_CurNode[i] < 0)
2150           ++done_count;
2151 
2152         if (m_CurNode[i] >= 0 && m_HoTNodeArr[i][m_CurNode[i]].m_hot.layer == layer) {
2153           // Skip the overlap hits -- if they exist.
2154           // 1. Overlap hit gets placed *after* the original hit in TrackCand::exportTrack()
2155           // which is *before* in the reverse iteration that we are doing here.
2156           // 2. Seed-hit merging can result in more than two hits per layer.
2157           // while (m_CurHit[i] > 0 && m_HoTArr[ i ][ m_CurHit[i] - 1 ].layer == layer) --m_CurHit[i];
2158           while (m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx >= 0 &&
2159                  m_HoTNodeArr[i][m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx].m_hot.layer == layer)
2160             m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
2161 
2162           const Hit &hit = L.refHit(m_HoTNodeArr[i][m_CurNode[i]].m_hot.index);
2163 
2164 #ifdef DEBUG_BACKWARD_FIT
2165           last_hit_ptr[i] = &hit;
2166 #endif
2167           m_msErr.copyIn(i, hit.errArray());
2168           m_msPar.copyIn(i, hit.posArray());
2169           ++here_count;
2170 
2171           m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
2172         } else {
2173 #ifdef DEBUG_BACKWARD_FIT
2174           last_hit_ptr[i] = nullptr;
2175 #endif
2176           no_mat_effs[i] = 1;
2177           tmp_pos[0] = m_Par[iC](i, 0, 0);
2178           tmp_pos[1] = m_Par[iC](i, 1, 0);
2179           tmp_pos[2] = m_Par[iC](i, 2, 0);
2180           m_msErr.copyIn(i, tmp_err);
2181           m_msPar.copyIn(i, tmp_pos);
2182         }
2183       }
2184 
2185       if (done_count == N_proc)
2186         break;
2187       if (here_count == 0)
2188         continue;
2189 
2190       // ZZZ Could add missing hits here, only if there are any actual matches.
2191 
2192       clearFailFlag();
2193 
2194       // PROP-FAIL-ENABLE Once always "copy input to output on fail" is removed from
2195       // propagateToR one might want to enable this for barrel or endcap or both.
2196       if (LI.is_barrel()) {
2197         propagateTracksToHitR(m_msPar, N_proc, m_prop_config->backward_fit_pflags, &no_mat_effs);
2198 
2199         kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
2200                         m_Err[iP],
2201                         m_Par[iP],
2202                         m_msErr,
2203                         m_msPar,
2204                         m_Err[iC],
2205                         m_Par[iC],
2206                         tmp_chi2,
2207                         N_proc);
2208       } else {
2209         propagateTracksToHitZ(m_msPar, N_proc, m_prop_config->backward_fit_pflags, &no_mat_effs);
2210 
2211         kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params,
2212                               m_Err[iP],
2213                               m_Par[iP],
2214                               m_msErr,
2215                               m_msPar,
2216                               m_Err[iC],
2217                               m_Par[iC],
2218                               tmp_chi2,
2219                               N_proc);
2220       }
2221 
2222 #if defined(DEBUG_PROP_UPDATE)
2223       printf("\nbkfit at layer %d, track in slot %d -- fail=%d, had hit=%d (%g, %g, %g)\n",
2224              LI.layer_id(),
2225              DSLOT,
2226              m_FailFlag[DSLOT],
2227              1 - no_mat_effs[DSLOT],
2228              m_msPar(DSLOT, 0, 0),
2229              m_msPar(DSLOT, 1, 0),
2230              m_msPar(DSLOT, 2, 0));
2231       printf("Propagated:\n");
2232       print_par_err(iP, DSLOT);
2233       printf("Updated:\n");
2234       print_par_err(iC, DSLOT);
2235 #endif
2236 
2237       // Fixup for failed propagation or invpt sign and charge.
2238       for (int i = 0; i < NN; ++i) {
2239         // PROP-FAIL-ENABLE The following to be enabled when propagation failure
2240         // detection is properly implemented in propagate-to-R/Z.
2241         // 1. The following code was only expecting barrel state to be restored.
2242         //      auto barrel_pf(m_prop_config->backward_fit_pflags);
2243         //      barrel_pf.copy_input_state_on_fail = true;
2244         // 2. There is also check on chi2, commented out to keep physics changes minimal.
2245         /*
2246         if (m_FailFlag[i] && LI.is_barrel()) {
2247           // Barrel pflags are set to include PF_copy_input_state_on_fail.
2248           // Endcap errors are immaterial here (relevant for fwd search), with prop error codes
2249           // one could do other things.
2250           // Are there also fail conditions in KalmanUpdate?
2251 #ifdef DEBUG
2252           if (debug && g_debug) {
2253             dprintf("MkFinder::bkFitFitTracks prop fail: chi2=%f, layer=%d, label=%d. Recovering.\n",
2254                     tmp_chi2[i], LI.layer_id(), m_Label[i]);
2255             print_par_err(iC, i);
2256           }
2257 #endif
2258           m_Err[iC].copySlot(i, m_Err[iP]);
2259           m_Par[iC].copySlot(i, m_Par[iP]);
2260         } else if (tmp_chi2[i] > 200 || tmp_chi2[i] < 0) {
2261 #ifdef DEBUG
2262           if (debug && g_debug) {
2263             dprintf("MkFinder::bkFitFitTracks chi2 fail: chi2=%f, layer=%d, label=%d. Recovering.\n",
2264                     tmp_chi2[i], LI.layer_id(), m_Label[i]);
2265             print_par_err(iC, i);
2266           }
2267 #endif
2268           // Go back to propagated state (at the current hit, the previous one is lost).
2269           m_Err[iC].copySlot(i, m_Err[iP]);
2270           m_Par[iC].copySlot(i, m_Par[iP]);
2271         }
2272         */
2273         // Fixup invpt sign and charge.
2274         if (i < N_proc && m_Par[iC].At(i, 3, 0) < 0) {
2275           m_Chg.At(i, 0, 0) = -m_Chg.At(i, 0, 0);
2276           m_Par[iC].At(i, 3, 0) = -m_Par[iC].At(i, 3, 0);
2277         }
2278       }
2279 
2280 #if defined(DEBUG_BACKWARD_FIT)
2281       // clang-format off
2282       bool debug = true;
2283       const char beg_cur_sep = '/'; // set to ' ' root parsable printouts
2284       for (int i = 0; i < N_proc; ++i) {
2285         if (chiDebug && last_hit_ptr[i]) {
2286           TrackCand &bb = *m_TrkCand[i];
2287           int ti = iP;
2288           float chi = tmp_chi2.At(i, 0, 0);
2289           float chi_prnt = std::isfinite(chi) ? chi : -9;
2290 
2291 #if defined(MKFIT_STANDALONE)
2292           const MCHitInfo &mchi = m_event->simHitsInfo_[last_hit_ptr[i]->mcHitID()];
2293 
2294           dprintf("BKF_OVERLAP %d %d %d %d %d %d %d "
2295                   "%f%c%f %f %f%c%f %f %f %f %d %d %d %d "
2296                   "%f %f %f %f %f\n",
2297               m_event->evtID(),
2298 #else
2299           dprintf("BKF_OVERLAP %d %d %d %d %d %d "
2300                   "%f%c%f %f %f%c%f %f %f %f %d %d %d "
2301                   "%f %f %f %f %f\n",
2302 #endif
2303               bb.label(), (int)bb.prodType(), bb.isFindable(),
2304               layer, L.is_stereo(), L.is_barrel(),
2305               bb.pT(), beg_cur_sep, 1.0f / m_Par[ti].At(i, 3, 0),
2306               bb.posEta(),
2307               bb.posPhi(), beg_cur_sep, std::atan2(m_Par[ti].At(i, 1, 0), m_Par[ti].At(i, 0, 0)),
2308               std::hypot(m_Par[ti].At(i, 0, 0), m_Par[ti].At(i, 1, 0)),
2309               m_Par[ti].At(i, 2, 0),
2310               chi_prnt,
2311               std::isnan(chi), std::isfinite(chi), chi > 0,
2312 #if defined(MKFIT_STANDALONE)
2313               mchi.mcTrackID(),
2314 #endif
2315               // The following three can get negative / prouce nans in e2s.
2316               // std::abs the args for FPE hunt.
2317               e2s(std::abs(m_Err[ti].At(i, 0, 0))),
2318               e2s(std::abs(m_Err[ti].At(i, 1, 1))),
2319               e2s(std::abs(m_Err[ti].At(i, 2, 2))),  // sx_t sy_t sz_t -- track errors
2320               1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0),
2321                                 m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)),  // d_xy
2322               1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0))             // d_z
2323           );
2324         }
2325       }
2326       // clang-format on
2327 #endif
2328 
2329       // update chi2
2330       m_Chi2.add(tmp_chi2);
2331     }
2332   }
2333 
2334   //------------------------------------------------------------------------------
2335 
2336   void MkFinder::bkFitPropTracksToPCA(const int N_proc) {
2337     propagateTracksToPCAZ(N_proc, m_prop_config->pca_prop_pflags);
2338   }
2339 
2340 }  // end namespace mkfit