Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-17 23:57:50

0001 #include "MkFinder.h"
0002 
0003 #include "CandCloner.h"
0004 #include "RecoTracker/MkFitCore/interface/IterationConfig.h"
0005 #include "FindingFoos.h"
0006 
0007 #include "KalmanUtilsMPlex.h"
0008 
0009 #include "MatriplexPackers.h"
0010 
0011 //#define DEBUG
0012 #include "Debug.h"
0013 
0014 #if (defined(DUMPHITWINDOW) || defined(DEBUG_BACKWARD_FIT)) && defined(MKFIT_STANDALONE)
0015 #include "RecoTracker/MkFitCore/standalone/Event.h"
0016 #endif
0017 
0018 #ifndef MKFIT_STANDALONE
0019 #include "FWCore/Utilities/interface/isFinite.h"
0020 #endif
0021 
0022 #include <algorithm>
0023 
0024 namespace {
0025   bool isFinite(float x) {
0026 #ifndef MKFIT_STANDALONE
0027     return edm::isFinite(x);
0028 #else
0029     return std::isfinite(x);
0030 #endif
0031   }
0032 }  // namespace
0033 
0034 namespace mkfit {
0035 
0036   void MkFinder::setup(const PropagationConfig &pc,
0037                        const IterationParams &ip,
0038                        const IterationLayerConfig &ilc,
0039                        const std::vector<bool> *ihm) {
0040     m_prop_config = &pc;
0041     m_iteration_params = &ip;
0042     m_iteration_layer_config = &ilc;
0043     m_iteration_hit_mask = ihm;
0044   }
0045 
0046   void MkFinder::setup_bkfit(const PropagationConfig &pc) { m_prop_config = &pc; }
0047 
0048   void MkFinder::release() {
0049     m_prop_config = nullptr;
0050     m_iteration_params = nullptr;
0051     m_iteration_layer_config = nullptr;
0052     m_iteration_hit_mask = nullptr;
0053   }
0054 
0055   //==============================================================================
0056   // Input / Output TracksAndHitIdx
0057   //==============================================================================
0058 
0059   void MkFinder::inputTracksAndHitIdx(const std::vector<Track> &tracks, int beg, int end, bool inputProp) {
0060     // Assign track parameters to initial state and copy hit values in.
0061 
0062     // This might not be true for the last chunk!
0063     // assert(end - beg == NN);
0064 
0065     const int iI = inputProp ? iP : iC;
0066 
0067     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0068       copy_in(tracks[i], imp, iI);
0069     }
0070   }
0071 
0072   void MkFinder::inputTracksAndHitIdx(
0073       const std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool inputProp, int mp_offset) {
0074     // Assign track parameters to initial state and copy hit values in.
0075 
0076     // This might not be true for the last chunk!
0077     // assert(end - beg == NN);
0078 
0079     const int iI = inputProp ? iP : iC;
0080 
0081     for (int i = beg, imp = mp_offset; i < end; ++i, ++imp) {
0082       copy_in(tracks[idxs[i]], imp, iI);
0083     }
0084   }
0085 
0086   void MkFinder::inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0087                                       const std::vector<std::pair<int, int>> &idxs,
0088                                       int beg,
0089                                       int end,
0090                                       bool inputProp) {
0091     // Assign track parameters to initial state and copy hit values in.
0092 
0093     // This might not be true for the last chunk!
0094     // assert(end - beg == NN);
0095 
0096     const int iI = inputProp ? iP : iC;
0097 
0098     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0099       const TrackCand &trk = tracks[idxs[i].first][idxs[i].second];
0100 
0101       copy_in(trk, imp, iI);
0102 
0103 #ifdef DUMPHITWINDOW
0104       m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].first].seed_algo();
0105       m_SeedLabel(imp, 0, 0) = tracks[idxs[i].first].seed_label();
0106 #endif
0107 
0108       m_SeedIdx(imp, 0, 0) = idxs[i].first;
0109       m_CandIdx(imp, 0, 0) = idxs[i].second;
0110     }
0111   }
0112 
0113   void MkFinder::inputTracksAndHits(const std::vector<CombCandidate> &tracks,
0114                                     const LayerOfHits &layer_of_hits,
0115                                     const std::vector<UpdateIndices> &idxs,
0116                                     int beg,
0117                                     int end,
0118                                     bool inputProp) {
0119     // Assign track parameters to initial state and copy hit values in.
0120 
0121     // This might not be true for the last chunk!
0122     // assert(end - beg == NN);
0123 
0124     const int iI = inputProp ? iP : iC;
0125 
0126     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0127       const TrackCand &trk = tracks[idxs[i].seed_idx][idxs[i].cand_idx];
0128 
0129       copy_in(trk, imp, iI);
0130 
0131 #ifdef DUMPHITWINDOW
0132       m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].seed_idx].seed_algo();
0133       m_SeedLabel(imp, 0, 0) = tracks[idxs[i].seed_idx.seed_label();
0134 #endif
0135 
0136       m_SeedIdx(imp, 0, 0) = idxs[i].seed_idx;
0137       m_CandIdx(imp, 0, 0) = idxs[i].cand_idx;
0138 
0139       const Hit &hit = layer_of_hits.refHit(idxs[i].hit_idx);
0140       m_msErr.copyIn(imp, hit.errArray());
0141       m_msPar.copyIn(imp, hit.posArray());
0142     }
0143   }
0144 
0145   void MkFinder::inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0146                                       const std::vector<std::pair<int, IdxChi2List>> &idxs,
0147                                       int beg,
0148                                       int end,
0149                                       bool inputProp) {
0150     // Assign track parameters to initial state and copy hit values in.
0151 
0152     // This might not be true for the last chunk!
0153     // assert(end - beg == NN);
0154 
0155     const int iI = inputProp ? iP : iC;
0156 
0157     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0158       const TrackCand &trk = tracks[idxs[i].first][idxs[i].second.trkIdx];
0159 
0160       copy_in(trk, imp, iI);
0161 
0162 #ifdef DUMPHITWINDOW
0163       m_SeedAlgo(imp, 0, 0) = tracks[idxs[i].first].seed_algo();
0164       m_SeedLabel(imp, 0, 0) = tracks[idxs[i].first].seed_label();
0165 #endif
0166 
0167       m_SeedIdx(imp, 0, 0) = idxs[i].first;
0168       m_CandIdx(imp, 0, 0) = idxs[i].second.trkIdx;
0169     }
0170   }
0171 
0172   void MkFinder::outputTracksAndHitIdx(std::vector<Track> &tracks, int beg, int end, bool outputProp) const {
0173     // Copies requested track parameters into Track objects.
0174     // The tracks vector should be resized to allow direct copying.
0175 
0176     const int iO = outputProp ? iP : iC;
0177 
0178     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0179       copy_out(tracks[i], imp, iO);
0180     }
0181   }
0182 
0183   void MkFinder::outputTracksAndHitIdx(
0184       std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool outputProp) const {
0185     // Copies requested track parameters into Track objects.
0186     // The tracks vector should be resized to allow direct copying.
0187 
0188     const int iO = outputProp ? iP : iC;
0189 
0190     for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0191       copy_out(tracks[idxs[i]], imp, iO);
0192     }
0193   }
0194 
0195   //==============================================================================
0196   // getHitSelDynamicWindows
0197   //==============================================================================
0198   // From HitSelectionWindows.h: track-related config on hit selection windows
0199 
0200   void MkFinder::getHitSelDynamicWindows(
0201       const float invpt, const float theta, float &min_dq, float &max_dq, float &min_dphi, float &max_dphi) {
0202     const IterationLayerConfig &ILC = *m_iteration_layer_config;
0203 
0204     float max_invpt = invpt;
0205     if (invpt > 10.0)
0206       max_invpt = 10.0;  // => pT>0.1 GeV
0207 
0208     // dq hit selection window
0209     float this_dq = (ILC.c_dq_0) * max_invpt + (ILC.c_dq_1) * theta + (ILC.c_dq_2);
0210     // In case layer is missing (e.g., seeding layers, or too low stats for training), leave original limits
0211     if ((ILC.c_dq_sf) * this_dq > 0.f) {
0212       min_dq = (ILC.c_dq_sf) * this_dq;
0213       max_dq = 2.0f * min_dq;
0214     }
0215 
0216     // dphi hit selection window
0217     float this_dphi = (ILC.c_dp_0) * max_invpt + (ILC.c_dp_1) * theta + (ILC.c_dp_2);
0218     // In case layer is missing (e.g., seeding layers, or too low stats for training), leave original limits
0219     if ((ILC.c_dp_sf) * this_dphi > min_dphi) {
0220       min_dphi = (ILC.c_dp_sf) * this_dphi;
0221       max_dphi = 2.0f * min_dphi;
0222     }
0223 
0224     //// For future optimization: for layer & iteration dependend hit chi2 cut
0225     //float this_c2 = (ILC.c_c2_0)*invpt+(ILC.c_c2_1)*theta+(ILC.c_c2_2);
0226     //// In case layer is missing (e.g., seeding layers, or too low stats for training), leave original limits
0227     //if(this_c2>0.f){
0228     //  max_c2 = (ILC.c_c2_sf)*this_c2;
0229     //}
0230   }
0231 
0232   //==============================================================================
0233   // getHitSelDynamicChi2Cut
0234   //==============================================================================
0235   // From HitSelectionWindows.h: track-related config on hit selection windows
0236 
0237   inline float MkFinder::getHitSelDynamicChi2Cut(const int itrk, const int ipar) {
0238     const IterationLayerConfig &ILC = *m_iteration_layer_config;
0239 
0240     const float minChi2Cut = m_iteration_params->chi2Cut_min;
0241     const float invpt = m_Par[ipar].At(itrk, 3, 0);
0242     const float theta = std::abs(m_Par[ipar].At(itrk, 5, 0) - Const::PIOver2);
0243 
0244     float max_invpt = invpt;
0245     if (invpt > 10.0)
0246       max_invpt = 10.0;
0247 
0248     float this_c2 = ILC.c_c2_0 * max_invpt + ILC.c_c2_1 * theta + ILC.c_c2_2;
0249     // In case layer is missing (e.g., seeding layers, or too low stats for training), leave original limits
0250     if ((ILC.c_c2_sf) * this_c2 > minChi2Cut)
0251       return ILC.c_c2_sf * this_c2;
0252     else
0253       return minChi2Cut;
0254   }
0255 
0256   //==============================================================================
0257   // SelectHitIndices
0258   //==============================================================================
0259 
0260   void MkFinder::selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc) {
0261     // bool debug = true;
0262     using bidx_t = LayerOfHits::bin_index_t;
0263     using bcnt_t = LayerOfHits::bin_content_t;
0264     const LayerOfHits &L = layer_of_hits;
0265     const IterationLayerConfig &ILC = *m_iteration_layer_config;
0266 
0267     const int iI = iP;
0268     const float nSigmaPhi = 3;
0269     const float nSigmaZ = 3;
0270     const float nSigmaR = 3;
0271 
0272     dprintf("LayerOfHits::SelectHitIndices %s layer=%d N_proc=%d\n",
0273             L.is_barrel() ? "barrel" : "endcap",
0274             L.layer_id(),
0275             N_proc);
0276 
0277     float dqv[NN], dphiv[NN], qv[NN], phiv[NN];
0278     bidx_t qb1v[NN], qb2v[NN], qbv[NN], pb1v[NN], pb2v[NN];
0279 
0280     const auto assignbins = [&](int itrack,
0281                                 float q,
0282                                 float dq,
0283                                 float phi,
0284                                 float dphi,
0285                                 float min_dq,
0286                                 float max_dq,
0287                                 float min_dphi,
0288                                 float max_dphi) {
0289       dphi = std::clamp(std::abs(dphi), min_dphi, max_dphi);
0290       dq = std::clamp(dq, min_dq, max_dq);
0291       //
0292       qv[itrack] = q;
0293       phiv[itrack] = phi;
0294       dphiv[itrack] = dphi;
0295       dqv[itrack] = dq;
0296       //
0297       qbv[itrack] = L.qBinChecked(q);
0298       qb1v[itrack] = L.qBinChecked(q - dq);
0299       qb2v[itrack] = L.qBinChecked(q + dq) + 1;
0300       pb1v[itrack] = L.phiBinChecked(phi - dphi);
0301       pb2v[itrack] = L.phiMaskApply(L.phiBin(phi + dphi) + 1);
0302     };
0303 
0304     const auto calcdphi2 = [&](int itrack, float dphidx, float dphidy) {
0305       return dphidx * dphidx * m_Err[iI].constAt(itrack, 0, 0) + dphidy * dphidy * m_Err[iI].constAt(itrack, 1, 1) +
0306              2 * dphidx * dphidy * m_Err[iI].constAt(itrack, 0, 1);
0307     };
0308 
0309     const auto calcdphi = [&](float dphi2, float min_dphi) {
0310       return std::max(nSigmaPhi * std::sqrt(std::abs(dphi2)), min_dphi);
0311     };
0312 
0313     if (L.is_barrel()) {
0314       // Pull out the part of the loop that vectorizes
0315 #pragma omp simd
0316       for (int itrack = 0; itrack < NN; ++itrack) {
0317         m_XHitSize[itrack] = 0;
0318 
0319         float min_dq = ILC.min_dq();
0320         float max_dq = ILC.max_dq();
0321         float min_dphi = ILC.min_dphi();
0322         float max_dphi = ILC.max_dphi();
0323 
0324         const float invpt = m_Par[iI].At(itrack, 3, 0);
0325         const float theta = std::fabs(m_Par[iI].At(itrack, 5, 0) - Const::PIOver2);
0326         getHitSelDynamicWindows(invpt, theta, min_dq, max_dq, min_dphi, max_dphi);
0327 
0328         const float x = m_Par[iI].constAt(itrack, 0, 0);
0329         const float y = m_Par[iI].constAt(itrack, 1, 0);
0330         const float r2 = x * x + y * y;
0331         const float dphidx = -y / r2, dphidy = x / r2;
0332         const float dphi2 = calcdphi2(itrack, dphidx, dphidy);
0333 #ifdef HARD_CHECK
0334         assert(dphi2 >= 0);
0335 #endif
0336 
0337         const float phi = getPhi(x, y);
0338         float dphi = calcdphi(dphi2, min_dphi);
0339 
0340         const float z = m_Par[iI].constAt(itrack, 2, 0);
0341         const float dz = std::abs(nSigmaZ * std::sqrt(m_Err[iI].constAt(itrack, 2, 2)));
0342         const float edgeCorr = std::abs(0.5f * (L.layer_info()->rout() - L.layer_info()->rin()) /
0343                                         std::tan(m_Par[iI].constAt(itrack, 5, 0)));
0344         // XXX-NUM-ERR above, m_Err(2,2) gets negative!
0345 
0346         m_XWsrResult[itrack] = L.is_within_z_sensitive_region(z, std::sqrt(dz * dz + edgeCorr * edgeCorr));
0347         assignbins(itrack, z, dz, phi, dphi, min_dq, max_dq, min_dphi, max_dphi);
0348       }
0349     } else  // endcap
0350     {
0351       //layer half-thikness for dphi spread calculation; only for very restrictive iters
0352       const float layerD = std::abs(L.layer_info()->zmax() - L.layer_info()->zmin()) * 0.5f *
0353                            (m_iteration_params->maxConsecHoles == 0 || m_iteration_params->maxHolesPerCand == 0);
0354       // Pull out the part of the loop that vectorizes
0355 #pragma omp simd
0356       for (int itrack = 0; itrack < NN; ++itrack) {
0357         m_XHitSize[itrack] = 0;
0358 
0359         float min_dq = ILC.min_dq();
0360         float max_dq = ILC.max_dq();
0361         float min_dphi = ILC.min_dphi();
0362         float max_dphi = ILC.max_dphi();
0363 
0364         const float invpt = m_Par[iI].At(itrack, 3, 0);
0365         const float theta = std::fabs(m_Par[iI].At(itrack, 5, 0) - Const::PIOver2);
0366         getHitSelDynamicWindows(invpt, theta, min_dq, max_dq, min_dphi, max_dphi);
0367 
0368         const float x = m_Par[iI].constAt(itrack, 0, 0);
0369         const float y = m_Par[iI].constAt(itrack, 1, 0);
0370         const float r2 = x * x + y * y;
0371         const float r2Inv = 1.f / r2;
0372         const float dphidx = -y * r2Inv, dphidy = x * r2Inv;
0373         const float phi = getPhi(x, y);
0374         const float dphi2 =
0375             calcdphi2(itrack, dphidx, dphidy)
0376             //range from finite layer thickness
0377             + std::pow(layerD * std::tan(m_Par[iI].At(itrack, 5, 0)) * std::sin(m_Par[iI].At(itrack, 4, 0) - phi), 2) *
0378                   r2Inv;
0379 #ifdef HARD_CHECK
0380         assert(dphi2 >= 0);
0381 #endif
0382 
0383         float dphi = calcdphi(dphi2, min_dphi);
0384 
0385         const float r = std::sqrt(r2);
0386         const float dr = nSigmaR * std::sqrt(std::abs(x * x * m_Err[iI].constAt(itrack, 0, 0) +
0387                                                       y * y * m_Err[iI].constAt(itrack, 1, 1) +
0388                                                       2 * x * y * m_Err[iI].constAt(itrack, 0, 1)) /
0389                                              r2);
0390         const float edgeCorr = std::abs(0.5f * (L.layer_info()->zmax() - L.layer_info()->zmin()) *
0391                                         std::tan(m_Par[iI].constAt(itrack, 5, 0)));
0392 
0393         m_XWsrResult[itrack] = L.is_within_r_sensitive_region(r, std::sqrt(dr * dr + edgeCorr * edgeCorr));
0394         assignbins(itrack, r, dr, phi, dphi, min_dq, max_dq, min_dphi, max_dphi);
0395       }
0396     }
0397 
0398     // Vectorizing this makes it run slower!
0399     //#pragma omp simd
0400     for (int itrack = 0; itrack < N_proc; ++itrack) {
0401       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
0402         m_XHitSize[itrack] = -1;
0403         continue;
0404       }
0405 
0406       const bidx_t qb = qbv[itrack];
0407       const bidx_t qb1 = qb1v[itrack];
0408       const bidx_t qb2 = qb2v[itrack];
0409       const bidx_t pb1 = pb1v[itrack];
0410       const bidx_t pb2 = pb2v[itrack];
0411 
0412       // Used only by usePhiQArrays
0413       const float q = qv[itrack];
0414       const float phi = phiv[itrack];
0415       const float dphi = dphiv[itrack];
0416       const float dq = dqv[itrack];
0417       // clang-format off
0418       dprintf("  %2d/%2d: %6.3f %6.3f %6.6f %7.5f %3u %3u %4u %4u\n",
0419               L.layer_id(), itrack, q, phi, dq, dphi,
0420               qb1, qb2, pb1, pb2);
0421       // clang-format on
0422       // MT: One could iterate in "spiral" order, to pick hits close to the center.
0423       // http://stackoverflow.com/questions/398299/looping-in-a-spiral
0424       // This would then work best with relatively small bin sizes.
0425       // Or, set them up so I can always take 3x3 array around the intersection.
0426 
0427 #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE)
0428       int thisseedmcid = -999999;
0429       {
0430         int seedlabel = m_SeedLabel(itrack, 0, 0);
0431         TrackVec &seedtracks = m_event->seedTracks_;
0432         int thisidx = -999999;
0433         for (int i = 0; i < int(seedtracks.size()); ++i) {
0434           auto &thisseed = seedtracks[i];
0435           if (thisseed.label() == seedlabel) {
0436             thisidx = i;
0437             break;
0438           }
0439         }
0440         if (thisidx > -999999) {
0441           auto &seedtrack = m_event->seedTracks_[thisidx];
0442           std::vector<int> thismcHitIDs;
0443           seedtrack.mcHitIDsVec(m_event->layerHits_, m_event->simHitsInfo_, thismcHitIDs);
0444           if (std::adjacent_find(thismcHitIDs.begin(), thismcHitIDs.end(), std::not_equal_to<>()) ==
0445               thismcHitIDs.end()) {
0446             thisseedmcid = thismcHitIDs.at(0);
0447           }
0448         }
0449       }
0450 #endif
0451 
0452       for (bidx_t qi = qb1; qi != qb2; ++qi) {
0453         for (bidx_t pi = pb1; pi != pb2; pi = L.phiMaskApply(pi + 1)) {
0454           // Limit to central Q-bin
0455           if (qi == qb && L.isBinDead(pi, qi) == true) {
0456             dprint("dead module for track in layer=" << L.layer_id() << " qb=" << qi << " pi=" << pi << " q=" << q
0457                                                      << " phi=" << phi);
0458             m_XWsrResult[itrack].m_in_gap = true;
0459           }
0460 
0461           // MT: The following line is the biggest hog (4% total run time).
0462           // This comes from cache misses, I presume.
0463           // It might make sense to make first loop to extract bin indices
0464           // and issue prefetches at the same time.
0465           // Then enter vectorized loop to actually collect the hits in proper order.
0466 
0467           //SK: ~20x1024 bin sizes give mostly 1 hit per bin. Commented out for 128 bins or less
0468           // #pragma nounroll
0469           auto pbi = L.phiQBinContent(pi, qi);
0470           for (bcnt_t hi = pbi.begin(); hi < pbi.end(); ++hi) {
0471             // MT: Access into m_hit_zs and m_hit_phis is 1% run-time each.
0472 
0473             unsigned int hi_orig = L.getOriginalHitIndex(hi);
0474 
0475             if (m_iteration_hit_mask && (*m_iteration_hit_mask)[hi_orig]) {
0476               dprintf(
0477                   "Yay, denying masked hit on layer %u, hi %u, orig idx %u\n", L.layer_info()->layer_id(), hi, hi_orig);
0478               continue;
0479             }
0480 
0481             if (Config::usePhiQArrays) {
0482               if (m_XHitSize[itrack] >= MPlexHitIdxMax)
0483                 break;
0484 
0485               const float ddq = std::abs(q - L.hit_q(hi));
0486               const float ddphi = cdist(std::abs(phi - L.hit_phi(hi)));
0487 
0488 #if defined(DUMPHITWINDOW) && defined(MKFIT_STANDALONE)
0489               {
0490                 const MCHitInfo &mchinfo = m_event->simHitsInfo_[L.refHit(hi).mcHitID()];
0491                 int mchid = mchinfo.mcTrackID();
0492                 int st_isfindable = 0;
0493                 int st_label = -999999;
0494                 int st_prodtype = 0;
0495                 int st_nhits = -1;
0496                 int st_charge = 0;
0497                 float st_r = -999.;
0498                 float st_z = -999.;
0499                 float st_pt = -999.;
0500                 float st_eta = -999.;
0501                 float st_phi = -999.;
0502                 if (mchid >= 0) {
0503                   Track simtrack = m_event->simTracks_[mchid];
0504                   st_isfindable = (int)simtrack.isFindable();
0505                   st_label = simtrack.label();
0506                   st_prodtype = (int)simtrack.prodType();
0507                   st_pt = simtrack.pT();
0508                   st_eta = simtrack.momEta();
0509                   st_phi = simtrack.momPhi();
0510                   st_nhits = simtrack.nTotalHits();
0511                   st_charge = simtrack.charge();
0512                   st_r = simtrack.posR();
0513                   st_z = simtrack.z();
0514                 }
0515 
0516                 const Hit &thishit = L.refHit(hi);
0517                 m_msErr.copyIn(itrack, thishit.errArray());
0518                 m_msPar.copyIn(itrack, thishit.posArray());
0519 
0520                 MPlexQF thisOutChi2;
0521                 MPlexLV tmpPropPar;
0522                 const FindingFoos &fnd_foos = FindingFoos::get_finding_foos(L.is_barrel());
0523                 (*fnd_foos.m_compute_chi2_foo)(m_Err[iI],
0524                                                m_Par[iI],
0525                                                m_Chg,
0526                                                m_msErr,
0527                                                m_msPar,
0528                                                thisOutChi2,
0529                                                tmpPropPar,
0530                                                N_proc,
0531                                                m_prop_config->finding_intra_layer_pflags,
0532                                                m_prop_config->finding_requires_propagation_to_hit_pos);
0533 
0534                 float hx = thishit.x();
0535                 float hy = thishit.y();
0536                 float hz = thishit.z();
0537                 float hr = std::hypot(hx, hy);
0538                 float hphi = std::atan2(hy, hx);
0539                 float hex = std::sqrt(thishit.exx());
0540                 if (std::isnan(hex))
0541                   hex = -999.;
0542                 float hey = std::sqrt(thishit.eyy());
0543                 if (std::isnan(hey))
0544                   hey = -999.;
0545                 float hez = std::sqrt(thishit.ezz());
0546                 if (std::isnan(hez))
0547                   hez = -999.;
0548                 float her = std::sqrt(
0549                     (hx * hx * thishit.exx() + hy * hy * thishit.eyy() + 2.0f * hx * hy * m_msErr.At(itrack, 0, 1)) /
0550                     (hr * hr));
0551                 if (std::isnan(her))
0552                   her = -999.;
0553                 float hephi = std::sqrt(thishit.ephi());
0554                 if (std::isnan(hephi))
0555                   hephi = -999.;
0556                 float hchi2 = thisOutChi2[itrack];
0557                 if (std::isnan(hchi2))
0558                   hchi2 = -999.;
0559                 float tx = m_Par[iI].At(itrack, 0, 0);
0560                 float ty = m_Par[iI].At(itrack, 1, 0);
0561                 float tz = m_Par[iI].At(itrack, 2, 0);
0562                 float tr = std::hypot(tx, ty);
0563                 float tphi = std::atan2(ty, tx);
0564                 float tchi2 = m_Chi2(itrack, 0, 0);
0565                 if (std::isnan(tchi2))
0566                   tchi2 = -999.;
0567                 float tex = std::sqrt(m_Err[iI].At(itrack, 0, 0));
0568                 if (std::isnan(tex))
0569                   tex = -999.;
0570                 float tey = std::sqrt(m_Err[iI].At(itrack, 1, 1));
0571                 if (std::isnan(tey))
0572                   tey = -999.;
0573                 float tez = std::sqrt(m_Err[iI].At(itrack, 2, 2));
0574                 if (std::isnan(tez))
0575                   tez = -999.;
0576                 float ter = std::sqrt(
0577                     (tx * tx * tex * tex + ty * ty * tey * tey + 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) /
0578                     (tr * tr));
0579                 if (std::isnan(ter))
0580                   ter = -999.;
0581                 float tephi = std::sqrt(
0582                     (ty * ty * tex * tex + tx * tx * tey * tey - 2.0f * tx * ty * m_Err[iI].At(itrack, 0, 1)) /
0583                     (tr * tr * tr * tr));
0584                 if (std::isnan(tephi))
0585                   tephi = -999.;
0586                 float ht_dxy = std::hypot(hx - tx, hy - ty);
0587                 float ht_dz = hz - tz;
0588                 float ht_dphi = cdist(std::abs(hphi - tphi));
0589 
0590                 static bool first = true;
0591                 if (first) {
0592                   printf(
0593                       "HITWINDOWSEL "
0594                       "evt_id/I:"
0595                       "lyr_id/I:lyr_isbrl/I:hit_idx/I:"
0596                       "trk_cnt/I:trk_idx/I:trk_label/I:"
0597                       "trk_pt/F:trk_eta/F:trk_mphi/F:trk_chi2/F:"
0598                       "nhits/I:"
0599                       "seed_idx/I:seed_label/I:seed_algo/I:seed_mcid/I:"
0600                       "hit_mcid/I:"
0601                       "st_isfindable/I:st_prodtype/I:st_label/I:"
0602                       "st_pt/F:st_eta/F:st_phi/F:"
0603                       "st_nhits/I:st_charge/I:st_r/F:st_z/F:"
0604                       "trk_q/F:hit_q/F:dq_trkhit/F:dq_cut/F:trk_phi/F:hit_phi/F:dphi_trkhit/F:dphi_cut/F:"
0605                       "t_x/F:t_y/F:t_r/F:t_phi/F:t_z/F:"
0606                       "t_ex/F:t_ey/F:t_er/F:t_ephi/F:t_ez/F:"
0607                       "h_x/F:h_y/F:h_r/F:h_phi/F:h_z/F:"
0608                       "h_ex/F:h_ey/F:h_er/F:h_ephi/F:h_ez/F:"
0609                       "ht_dxy/F:ht_dz/F:ht_dphi/F:"
0610                       "h_chi2/F"
0611                       "\n");
0612                   first = false;
0613                 }
0614 
0615                 if (!(std::isnan(phi)) && !(std::isnan(getEta(m_Par[iI].At(itrack, 5, 0))))) {
0616                   //|| std::isnan(ter) || std::isnan(her) || std::isnan(m_Chi2(itrack, 0, 0)) || std::isnan(hchi2)))
0617                   // clang-format off
0618                   printf("HITWINDOWSEL "
0619                          "%d "
0620                          "%d %d %d "
0621                          "%d %d %d "
0622                          "%6.3f %6.3f %6.3f %6.3f "
0623                          "%d "
0624                          "%d %d %d %d "
0625                          "%d "
0626                          "%d %d %d "
0627                          "%6.3f %6.3f %6.3f "
0628                          "%d %d %6.3f %6.3f "
0629                          "%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f "
0630                          "%6.3f %6.3f %6.3f %6.3f %6.3f "
0631                          "%6.6f %6.6f %6.6f %6.6f %6.6f "
0632                          "%6.3f %6.3f %6.3f %6.3f %6.3f "
0633                          "%6.6f %6.6f %6.6f %6.6f %6.6f "
0634                          "%6.3f %6.3f %6.3f "
0635                          "%6.3f"
0636                          "\n",
0637                          m_event->evtID(),
0638                          L.layer_id(), L.is_barrel(), L.getOriginalHitIndex(hi),
0639                          itrack, m_CandIdx(itrack, 0, 0), m_Label(itrack, 0, 0),
0640                          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),
0641                          m_NFoundHits(itrack, 0, 0),
0642                          m_SeedIdx(itrack, 0, 0), m_SeedLabel(itrack, 0, 0), m_SeedAlgo(itrack, 0, 0), thisseedmcid,
0643                          mchid,
0644                          st_isfindable, st_prodtype, st_label,
0645                          st_pt, st_eta, st_phi,
0646                          st_nhits, st_charge, st_r, st_z,
0647                          q, L.hit_q(hi), ddq, dq, phi, L.hit_phi(hi), ddphi, dphi,
0648                          tx, ty, tr, tphi, tz,
0649                          tex, tey, ter, tephi, tez,
0650                          hx, hy, hr, hphi, hz,
0651                          hex, hey, her, hephi, hez,
0652                          ht_dxy, ht_dz, ht_dphi,
0653                          hchi2);
0654                   // clang-format on
0655                 }
0656               }
0657 #endif
0658 
0659               if (ddq >= dq)
0660                 continue;
0661               if (ddphi >= dphi)
0662                 continue;
0663               // clang-format off
0664               dprintf("     SHI %3u %4u %5u  %6.3f %6.3f %6.4f %7.5f   %s\n",
0665                       qi, pi, hi, L.hit_q(hi), L.hit_phi(hi),
0666                       ddq, ddphi, (ddq < dq && ddphi < dphi) ? "PASS" : "FAIL");
0667               // clang-format on
0668 
0669               // MT: Removing extra check gives full efficiency ...
0670               //     and means our error estimations are wrong!
0671               // Avi says we should have *minimal* search windows per layer.
0672               // Also ... if bins are sufficiently small, we do not need the extra
0673               // checks, see above.
0674               m_XHitArr.At(itrack, m_XHitSize[itrack]++, 0) = hi_orig;
0675             } else {
0676               // MT: The following check alone makes more sense with spiral traversal,
0677               // we'd be taking in closest hits first.
0678 
0679               // Hmmh -- there used to be some more checks here.
0680               // Or, at least, the phi binning was much smaller and no further checks were done.
0681               assert(false && "this code has not been used in a while -- see comments in code");
0682 
0683               if (m_XHitSize[itrack] < MPlexHitIdxMax) {
0684                 m_XHitArr.At(itrack, m_XHitSize[itrack]++, 0) = hi_orig;
0685               }
0686             }
0687           }  //hi
0688         }    //pi
0689       }      //qi
0690     }        //itrack
0691   }
0692 
0693   //==============================================================================
0694   // AddBestHit - Best Hit Track Finding
0695   //==============================================================================
0696 
0697   void MkFinder::addBestHit(const LayerOfHits &layer_of_hits, const int N_proc, const FindingFoos &fnd_foos) {
0698     // debug = true;
0699 
0700     MatriplexHitPacker mhp(layer_of_hits.hitArray());
0701 
0702     float minChi2[NN];
0703     int bestHit[NN];
0704     int maxSize = 0;
0705 
0706     // Determine maximum number of hits for tracks in the collection.
0707     for (int it = 0; it < NN; ++it) {
0708       if (it < N_proc) {
0709         if (m_XHitSize[it] > 0) {
0710           maxSize = std::max(maxSize, m_XHitSize[it]);
0711         }
0712       }
0713 
0714       bestHit[it] = -1;
0715       minChi2[it] = getHitSelDynamicChi2Cut(it, iP);
0716     }
0717 
0718     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
0719       //fixme what if size is zero???
0720 
0721       mhp.reset();
0722 
0723 #pragma omp simd
0724       for (int itrack = 0; itrack < N_proc; ++itrack) {
0725         if (hit_cnt < m_XHitSize[itrack]) {
0726           mhp.addInputAt(itrack, layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)));
0727         }
0728       }
0729 
0730       mhp.pack(m_msErr, m_msPar);
0731 
0732       //now compute the chi2 of track state vs hit
0733       MPlexQF outChi2;
0734       MPlexLV tmpPropPar;
0735       (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
0736                                      m_Par[iP],
0737                                      m_Chg,
0738                                      m_msErr,
0739                                      m_msPar,
0740                                      outChi2,
0741                                      tmpPropPar,
0742                                      N_proc,
0743                                      m_prop_config->finding_intra_layer_pflags,
0744                                      m_prop_config->finding_requires_propagation_to_hit_pos);
0745 
0746       //update best hit in case chi2<minChi2
0747 #pragma omp simd
0748       for (int itrack = 0; itrack < N_proc; ++itrack) {
0749         if (hit_cnt < m_XHitSize[itrack]) {
0750           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
0751           dprint("chi2=" << chi2 << " minChi2[itrack]=" << minChi2[itrack]);
0752           if (chi2 < minChi2[itrack]) {
0753             minChi2[itrack] = chi2;
0754             bestHit[itrack] = m_XHitArr.At(itrack, hit_cnt, 0);
0755           }
0756         }
0757       }
0758     }  // end loop over hits
0759 
0760     //#pragma omp simd
0761     for (int itrack = 0; itrack < N_proc; ++itrack) {
0762       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
0763         // Why am I doing this?
0764         m_msErr.setDiagonal3x3(itrack, 666);
0765         m_msPar(itrack, 0, 0) = m_Par[iP](itrack, 0, 0);
0766         m_msPar(itrack, 1, 0) = m_Par[iP](itrack, 1, 0);
0767         m_msPar(itrack, 2, 0) = m_Par[iP](itrack, 2, 0);
0768 
0769         // XXXX If not in gap, should get back the old track params. But they are gone ...
0770         // Would actually have to do it right after SelectHitIndices where updated params are still ok.
0771         // Here they got screwed during hit matching.
0772         // So, I'd store them there (into propagated params) and retrieve them here.
0773         // Or we decide not to care ...
0774 
0775         continue;
0776       }
0777 
0778       //fixme decide what to do in case no hit found
0779       if (bestHit[itrack] >= 0) {
0780         const Hit &hit = layer_of_hits.refHit(bestHit[itrack]);
0781         const float chi2 = minChi2[itrack];
0782 
0783         dprint("ADD BEST HIT FOR TRACK #"
0784                << itrack << std::endl
0785                << "prop x=" << m_Par[iP].constAt(itrack, 0, 0) << " y=" << m_Par[iP].constAt(itrack, 1, 0) << std::endl
0786                << "copy in hit #" << bestHit[itrack] << " x=" << hit.position()[0] << " y=" << hit.position()[1]);
0787 
0788         m_msErr.copyIn(itrack, hit.errArray());
0789         m_msPar.copyIn(itrack, hit.posArray());
0790         m_Chi2(itrack, 0, 0) += chi2;
0791 
0792         add_hit(itrack, bestHit[itrack], layer_of_hits.layer_id());
0793       } else {
0794         int fake_hit_idx = Hit::kHitMissIdx;
0795 
0796         if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
0797           // YYYYYY Config::store_missed_layers
0798           fake_hit_idx = Hit::kHitEdgeIdx;
0799         } else if (num_all_minus_one_hits(itrack)) {
0800           fake_hit_idx = Hit::kHitStopIdx;
0801         }
0802 
0803         dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx)
0804                                           << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0)));
0805 
0806         m_msErr.setDiagonal3x3(itrack, 666);
0807         m_msPar(itrack, 0, 0) = m_Par[iP](itrack, 0, 0);
0808         m_msPar(itrack, 1, 0) = m_Par[iP](itrack, 1, 0);
0809         m_msPar(itrack, 2, 0) = m_Par[iP](itrack, 2, 0);
0810         // Don't update chi2
0811 
0812         add_hit(itrack, fake_hit_idx, layer_of_hits.layer_id());
0813       }
0814     }
0815 
0816     // Update the track parameters with this hit. (Note that some calculations
0817     // are already done when computing chi2. Not sure it's worth caching them?)
0818 
0819     dprint("update parameters");
0820     (*fnd_foos.m_update_param_foo)(m_Err[iP],
0821                                    m_Par[iP],
0822                                    m_Chg,
0823                                    m_msErr,
0824                                    m_msPar,
0825                                    m_Err[iC],
0826                                    m_Par[iC],
0827                                    N_proc,
0828                                    m_prop_config->finding_intra_layer_pflags,
0829                                    m_prop_config->finding_requires_propagation_to_hit_pos);
0830 
0831     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));
0832   }
0833 
0834   //=======================================================
0835   // isStripQCompatible : check if prop is consistent with the barrel/endcap strip length
0836   //=======================================================
0837   bool isStripQCompatible(
0838       int itrack, bool isBarrel, const MPlexLS &pErr, const MPlexLV &pPar, const MPlexHS &msErr, const MPlexHV &msPar) {
0839     //check module compatibility via long strip side = L/sqrt(12)
0840     if (isBarrel) {  //check z direction only
0841       const float res = std::abs(msPar.constAt(itrack, 2, 0) - pPar.constAt(itrack, 2, 0));
0842       const float hitHL = sqrt(msErr.constAt(itrack, 2, 2) * 3.f);  //half-length
0843       const float qErr = sqrt(pErr.constAt(itrack, 2, 2));
0844       dprint("qCompat " << hitHL << " + " << 3.f * qErr << " vs " << res);
0845       return hitHL + std::max(3.f * qErr, 0.5f) > res;
0846     } else {  //project on xy, assuming the strip Length >> Width
0847       const float res[2]{msPar.constAt(itrack, 0, 0) - pPar.constAt(itrack, 0, 0),
0848                          msPar.constAt(itrack, 1, 0) - pPar.constAt(itrack, 1, 0)};
0849       const float hitT2 = msErr.constAt(itrack, 0, 0) + msErr.constAt(itrack, 1, 1);
0850       const float hitT2inv = 1.f / hitT2;
0851       const float proj[3] = {msErr.constAt(itrack, 0, 0) * hitT2inv,
0852                              msErr.constAt(itrack, 0, 1) * hitT2inv,
0853                              msErr.constAt(itrack, 1, 1) * hitT2inv};
0854       const float qErr =
0855           sqrt(std::abs(pErr.constAt(itrack, 0, 0) * proj[0] + 2.f * pErr.constAt(itrack, 0, 1) * proj[1] +
0856                         pErr.constAt(itrack, 1, 1) * proj[2]));  //take abs to avoid non-pos-def cases
0857       const float resProj =
0858           sqrt(res[0] * proj[0] * res[0] + 2.f * res[1] * proj[1] * res[0] + res[1] * proj[2] * res[1]);
0859       dprint("qCompat " << sqrt(hitT2 * 3.f) << " + " << 3.f * qErr << " vs " << resProj);
0860       return sqrt(hitT2 * 3.f) + std::max(3.f * qErr, 0.5f) > resProj;
0861     }
0862   }
0863 
0864   //=======================================================
0865   // passStripChargePCMfromTrack : apply the slope correction to charge per cm and cut using hit err matrix
0866   //         the raw pcm = charge/L_normal
0867   //         the corrected qCorr = charge/L_path = charge/(L_normal*p/p_zLocal) = pcm*p_zLocal/p
0868   //=======================================================
0869   bool passStripChargePCMfromTrack(
0870       int itrack, bool isBarrel, unsigned int pcm, unsigned int pcmMin, const MPlexLV &pPar, const MPlexHS &msErr) {
0871     //skip the overflow case
0872     if (pcm >= Hit::maxChargePerCM())
0873       return true;
0874 
0875     float qSF;
0876     if (isBarrel) {  //project in x,y, assuming zero-error direction is in this plane
0877       const float hitT2 = msErr.constAt(itrack, 0, 0) + msErr.constAt(itrack, 1, 1);
0878       const float hitT2inv = 1.f / hitT2;
0879       const float proj[3] = {msErr.constAt(itrack, 0, 0) * hitT2inv,
0880                              msErr.constAt(itrack, 0, 1) * hitT2inv,
0881                              msErr.constAt(itrack, 1, 1) * hitT2inv};
0882       const bool detXY_OK =
0883           std::abs(proj[0] * proj[2] - proj[1] * proj[1]) < 0.1f;  //check that zero-direction is close
0884       const float cosP = cos(pPar.constAt(itrack, 4, 0));
0885       const float sinP = sin(pPar.constAt(itrack, 4, 0));
0886       const float sinT = std::abs(sin(pPar.constAt(itrack, 5, 0)));
0887       //qSF = sqrt[(px,py)*(1-proj)*(px,py)]/p = sinT*sqrt[(cosP,sinP)*(1-proj)*(cosP,sinP)].
0888       qSF = detXY_OK ? sinT * std::sqrt(std::abs(1.f + cosP * cosP * proj[0] + sinP * sinP * proj[2] -
0889                                                  2.f * cosP * sinP * proj[1]))
0890                      : 1.f;
0891     } else {  //project on z
0892       // p_zLocal/p = p_z/p = cosT
0893       qSF = std::abs(cos(pPar.constAt(itrack, 5, 0)));
0894     }
0895 
0896     const float qCorr = pcm * qSF;
0897     dprint("pcm " << pcm << " * " << qSF << " = " << qCorr << " vs " << pcmMin);
0898     return qCorr > pcmMin;
0899   }
0900 
0901   //==============================================================================
0902   // FindCandidates - Standard Track Finding
0903   //==============================================================================
0904 
0905   void MkFinder::findCandidates(const LayerOfHits &layer_of_hits,
0906                                 std::vector<std::vector<TrackCand>> &tmp_candidates,
0907                                 const int offset,
0908                                 const int N_proc,
0909                                 const FindingFoos &fnd_foos) {
0910     // bool debug = true;
0911 
0912     MatriplexHitPacker mhp(layer_of_hits.hitArray());
0913 
0914     int maxSize = 0;
0915 
0916     // Determine maximum number of hits for tracks in the collection.
0917     for (int it = 0; it < NN; ++it) {
0918       if (it < N_proc) {
0919         if (m_XHitSize[it] > 0) {
0920           maxSize = std::max(maxSize, m_XHitSize[it]);
0921         }
0922       }
0923     }
0924 
0925     dprintf("FindCandidates max hits to process=%d\n", maxSize);
0926 
0927     int nHitsAdded[NN]{};
0928     bool isTooLargeCluster[NN]{false};
0929 
0930     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
0931       mhp.reset();
0932 
0933       int charge_pcm[NN];
0934 
0935 #pragma omp simd
0936       for (int itrack = 0; itrack < N_proc; ++itrack) {
0937         if (hit_cnt < m_XHitSize[itrack]) {
0938           const auto &hit = layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0));
0939           mhp.addInputAt(itrack, hit);
0940           charge_pcm[itrack] = hit.chargePerCM();
0941         }
0942       }
0943 
0944       mhp.pack(m_msErr, m_msPar);
0945 
0946       //now compute the chi2 of track state vs hit
0947       MPlexQF outChi2;
0948       MPlexLV propPar;
0949       (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
0950                                      m_Par[iP],
0951                                      m_Chg,
0952                                      m_msErr,
0953                                      m_msPar,
0954                                      outChi2,
0955                                      propPar,
0956                                      N_proc,
0957                                      m_prop_config->finding_intra_layer_pflags,
0958                                      m_prop_config->finding_requires_propagation_to_hit_pos);
0959 
0960       // Now update the track parameters with this hit (note that some
0961       // calculations are already done when computing chi2, to be optimized).
0962       // 1. This is not needed for candidates the hit is not added to, but it's
0963       // vectorized so doing it serially below should take the same time.
0964       // 2. Still it's a waste of time in case the hit is not added to any of the
0965       // candidates, so check beforehand that at least one cand needs update.
0966       bool oneCandPassCut = false;
0967       for (int itrack = 0; itrack < N_proc; ++itrack) {
0968         float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
0969 
0970         if (hit_cnt < m_XHitSize[itrack]) {
0971           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
0972           dprint("chi2=" << chi2);
0973           if (chi2 < max_c2) {
0974             bool isCompatible = true;
0975             if (!layer_of_hits.is_pixel()) {
0976               //check module compatibility via long strip side = L/sqrt(12)
0977               isCompatible =
0978                   isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
0979 
0980               //rescale strip charge to track parameters and reapply the cut
0981               isCompatible &= passStripChargePCMfromTrack(
0982                   itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
0983             }
0984             // Select only SiStrip hits with cluster size < maxClusterSize
0985             if (!layer_of_hits.is_pixel()) {
0986               if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
0987                   m_iteration_params->maxClusterSize) {
0988                 isTooLargeCluster[itrack] = true;
0989                 isCompatible = false;
0990               }
0991             }
0992 
0993             if (isCompatible) {
0994               oneCandPassCut = true;
0995               break;
0996             }
0997           }
0998         }
0999       }
1000 
1001       if (oneCandPassCut) {
1002         MPlexQI tmpChg = m_Chg;
1003         (*fnd_foos.m_update_param_foo)(m_Err[iP],
1004                                        m_Par[iP],
1005                                        tmpChg,
1006                                        m_msErr,
1007                                        m_msPar,
1008                                        m_Err[iC],
1009                                        m_Par[iC],
1010                                        N_proc,
1011                                        m_prop_config->finding_intra_layer_pflags,
1012                                        m_prop_config->finding_requires_propagation_to_hit_pos);
1013 
1014         dprint("update parameters" << std::endl
1015                                    << "propagated track parameters x=" << m_Par[iP].constAt(0, 0, 0)
1016                                    << " y=" << m_Par[iP].constAt(0, 1, 0) << std::endl
1017                                    << "               hit position x=" << m_msPar.constAt(0, 0, 0)
1018                                    << " y=" << m_msPar.constAt(0, 1, 0) << std::endl
1019                                    << "   updated track parameters x=" << m_Par[iC].constAt(0, 0, 0)
1020                                    << " y=" << m_Par[iC].constAt(0, 1, 0));
1021 
1022         //create candidate with hit in case chi2 < max_c2
1023         //fixme: please vectorize me... (not sure it's possible in this case)
1024         for (int itrack = 0; itrack < N_proc; ++itrack) {
1025           float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
1026 
1027           if (hit_cnt < m_XHitSize[itrack]) {
1028             const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1029             dprint("chi2=" << chi2);
1030             if (chi2 < max_c2) {
1031               bool isCompatible = true;
1032               if (!layer_of_hits.is_pixel()) {
1033                 //check module compatibility via long strip side = L/sqrt(12)
1034                 isCompatible =
1035                     isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
1036 
1037                 //rescale strip charge to track parameters and reapply the cut
1038                 isCompatible &= passStripChargePCMfromTrack(
1039                     itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
1040               }
1041               // Select only SiStrip hits with cluster size < maxClusterSize
1042               if (!layer_of_hits.is_pixel()) {
1043                 if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
1044                     m_iteration_params->maxClusterSize)
1045                   isCompatible = false;
1046               }
1047 
1048               if (isCompatible) {
1049                 bool hitExists = false;
1050                 int maxHits = m_NFoundHits(itrack, 0, 0);
1051                 if (layer_of_hits.is_pixel()) {
1052                   for (int i = 0; i <= maxHits; ++i) {
1053                     if (i > 2)
1054                       break;
1055                     if (m_HoTArrs[itrack][i].layer == layer_of_hits.layer_id()) {
1056                       hitExists = true;
1057                       break;
1058                     }
1059                   }
1060                 }
1061                 if (hitExists)
1062                   continue;
1063 
1064                 nHitsAdded[itrack]++;
1065                 dprint("chi2 cut passed, creating new candidate");
1066                 // Create a new candidate and fill the tmp_candidates output vector.
1067                 // QQQ only instantiate if it will pass, be better than N_best
1068 
1069                 const int hit_idx = m_XHitArr.At(itrack, hit_cnt, 0);
1070 
1071                 TrackCand newcand;
1072                 copy_out(newcand, itrack, iC);
1073                 newcand.setCharge(tmpChg(itrack, 0, 0));
1074                 newcand.addHitIdx(hit_idx, layer_of_hits.layer_id(), chi2);
1075                 newcand.setScore(getScoreCand(newcand, true /*penalizeTailMissHits*/, true /*inFindCandidates*/));
1076                 newcand.setOriginIndex(m_CandIdx(itrack, 0, 0));
1077 
1078                 // To apply a fixed cut instead of dynamic cut for overlap: m_iteration_params->chi2CutOverlap
1079                 if (chi2 < max_c2) {
1080                   CombCandidate &ccand = *newcand.combCandidate();
1081                   ccand[m_CandIdx(itrack, 0, 0)].considerHitForOverlap(
1082                       hit_idx, layer_of_hits.refHit(hit_idx).detIDinLayer(), chi2);
1083                 }
1084 
1085                 dprint("updated track parameters x=" << newcand.parameters()[0] << " y=" << newcand.parameters()[1]
1086                                                      << " z=" << newcand.parameters()[2]
1087                                                      << " pt=" << 1. / newcand.parameters()[3]);
1088 
1089                 tmp_candidates[m_SeedIdx(itrack, 0, 0) - offset].emplace_back(newcand);
1090               }
1091             }
1092           }
1093         }
1094       }  //end if (oneCandPassCut)
1095 
1096     }  //end loop over hits
1097 
1098     //now add invalid hit
1099     //fixme: please vectorize me...
1100     for (int itrack = 0; itrack < N_proc; ++itrack) {
1101       // Cands that miss the layer are stashed away in MkBuilder(), before propagation,
1102       // and then merged back afterwards.
1103       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
1104         continue;
1105       }
1106 
1107       int fake_hit_idx = ((num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand) &&
1108                           (m_NTailMinusOneHits(itrack, 0, 0) < m_iteration_params->maxConsecHoles))
1109                              ? Hit::kHitMissIdx
1110                              : Hit::kHitStopIdx;
1111 
1112       if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
1113         // YYYYYY m_iteration_params->store_missed_layers
1114         fake_hit_idx = Hit::kHitEdgeIdx;
1115       }
1116       //now add fake hit for tracks that passsed through inactive modules
1117       else if (m_XWsrResult[itrack].m_in_gap == true && nHitsAdded[itrack] == 0) {
1118         fake_hit_idx = Hit::kHitInGapIdx;
1119       }
1120       //now add fake hit for cases where hit cluster size is larger than maxClusterSize
1121       else if (isTooLargeCluster[itrack] == true && nHitsAdded[itrack] == 0) {
1122         fake_hit_idx = Hit::kHitMaxClusterIdx;
1123       }
1124 
1125       dprint("ADD FAKE HIT FOR TRACK #" << itrack << " withinBounds=" << (fake_hit_idx != Hit::kHitEdgeIdx)
1126                                         << " r=" << std::hypot(m_Par[iP](itrack, 0, 0), m_Par[iP](itrack, 1, 0)));
1127 
1128       // QQQ as above, only create and add if score better
1129       TrackCand newcand;
1130       copy_out(newcand, itrack, iP);
1131       newcand.addHitIdx(fake_hit_idx, layer_of_hits.layer_id(), 0.);
1132       newcand.setScore(getScoreCand(newcand, true /*penalizeTailMissHits*/, true /*inFindCandidates*/));
1133       // Only relevant when we actually add a hit
1134       // newcand.setOriginIndex(m_CandIdx(itrack, 0, 0));
1135       tmp_candidates[m_SeedIdx(itrack, 0, 0) - offset].emplace_back(newcand);
1136     }
1137   }
1138 
1139   //==============================================================================
1140   // FindCandidatesCloneEngine - Clone Engine Track Finding
1141   //==============================================================================
1142 
1143   void MkFinder::findCandidatesCloneEngine(const LayerOfHits &layer_of_hits,
1144                                            CandCloner &cloner,
1145                                            const int offset,
1146                                            const int N_proc,
1147                                            const FindingFoos &fnd_foos) {
1148     // bool debug = true;
1149 
1150     MatriplexHitPacker mhp(layer_of_hits.hitArray());
1151 
1152     int maxSize = 0;
1153 
1154     // Determine maximum number of hits for tracks in the collection.
1155 #pragma omp simd
1156     for (int it = 0; it < NN; ++it) {
1157       if (it < N_proc) {
1158         if (m_XHitSize[it] > 0) {
1159           maxSize = std::max(maxSize, m_XHitSize[it]);
1160         }
1161       }
1162     }
1163 
1164     dprintf("FindCandidatesCloneEngine max hits to process=%d\n", maxSize);
1165 
1166     int nHitsAdded[NN]{};
1167     bool isTooLargeCluster[NN]{false};
1168 
1169     for (int hit_cnt = 0; hit_cnt < maxSize; ++hit_cnt) {
1170       mhp.reset();
1171 
1172       int charge_pcm[NN];
1173 
1174 #pragma omp simd
1175       for (int itrack = 0; itrack < N_proc; ++itrack) {
1176         if (hit_cnt < m_XHitSize[itrack]) {
1177           const auto &hit = layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0));
1178           mhp.addInputAt(itrack, hit);
1179           charge_pcm[itrack] = hit.chargePerCM();
1180         }
1181       }
1182 
1183       mhp.pack(m_msErr, m_msPar);
1184 
1185       //now compute the chi2 of track state vs hit
1186       MPlexQF outChi2;
1187       MPlexLV propPar;
1188       (*fnd_foos.m_compute_chi2_foo)(m_Err[iP],
1189                                      m_Par[iP],
1190                                      m_Chg,
1191                                      m_msErr,
1192                                      m_msPar,
1193                                      outChi2,
1194                                      propPar,
1195                                      N_proc,
1196                                      m_prop_config->finding_intra_layer_pflags,
1197                                      m_prop_config->finding_requires_propagation_to_hit_pos);
1198 
1199 #pragma omp simd  // DOES NOT VECTORIZE AS IT IS NOW
1200       for (int itrack = 0; itrack < N_proc; ++itrack) {
1201         // make sure the hit was in the compatiblity window for the candidate
1202 
1203         float max_c2 = getHitSelDynamicChi2Cut(itrack, iP);
1204 
1205         if (hit_cnt < m_XHitSize[itrack]) {
1206           // XXX-NUM-ERR assert(chi2 >= 0);
1207           const float chi2 = std::abs(outChi2[itrack]);  //fixme negative chi2 sometimes...
1208 
1209           dprint("chi2=" << chi2 << " for trkIdx=" << itrack << " hitIdx=" << m_XHitArr.At(itrack, hit_cnt, 0));
1210           if (chi2 < max_c2) {
1211             bool isCompatible = true;
1212             if (!layer_of_hits.is_pixel()) {
1213               //check module compatibility via long strip side = L/sqrt(12)
1214               isCompatible =
1215                   isStripQCompatible(itrack, layer_of_hits.is_barrel(), m_Err[iP], propPar, m_msErr, m_msPar);
1216 
1217               //rescale strip charge to track parameters and reapply the cut
1218               isCompatible &= passStripChargePCMfromTrack(
1219                   itrack, layer_of_hits.is_barrel(), charge_pcm[itrack], Hit::minChargePerCM(), propPar, m_msErr);
1220             }
1221 
1222             // Select only SiStrip hits with cluster size < maxClusterSize
1223             if (!layer_of_hits.is_pixel()) {
1224               if (layer_of_hits.refHit(m_XHitArr.At(itrack, hit_cnt, 0)).spanRows() >=
1225                   m_iteration_params->maxClusterSize) {
1226                 isTooLargeCluster[itrack] = true;
1227                 isCompatible = false;
1228               }
1229             }
1230 
1231             if (isCompatible) {
1232               CombCandidate &ccand = cloner.combCandWithOriginalIndex(m_SeedIdx(itrack, 0, 0));
1233               bool hitExists = false;
1234               int maxHits = m_NFoundHits(itrack, 0, 0);
1235               if (layer_of_hits.is_pixel()) {
1236                 for (int i = 0; i <= maxHits; ++i) {
1237                   if (i > 2)
1238                     break;
1239                   if (ccand.hot(i).layer == layer_of_hits.layer_id()) {
1240                     hitExists = true;
1241                     break;
1242                   }
1243                 }
1244               }
1245               if (hitExists)
1246                 continue;
1247 
1248               nHitsAdded[itrack]++;
1249               const int hit_idx = m_XHitArr.At(itrack, hit_cnt, 0);
1250 
1251               // Register hit for overlap consideration, if chi2 cut is passed
1252               // To apply a fixed cut instead of dynamic cut for overlap: m_iteration_params->chi2CutOverlap
1253               if (chi2 < max_c2) {
1254                 ccand[m_CandIdx(itrack, 0, 0)].considerHitForOverlap(
1255                     hit_idx, layer_of_hits.refHit(hit_idx).detIDinLayer(), chi2);
1256               }
1257 
1258               IdxChi2List tmpList;
1259               tmpList.trkIdx = m_CandIdx(itrack, 0, 0);
1260               tmpList.hitIdx = hit_idx;
1261               tmpList.module = layer_of_hits.refHit(hit_idx).detIDinLayer();
1262               tmpList.nhits = m_NFoundHits(itrack, 0, 0) + 1;
1263               tmpList.ntailholes = 0;
1264               tmpList.noverlaps = m_NOverlapHits(itrack, 0, 0);
1265               tmpList.nholes = num_all_minus_one_hits(itrack);
1266               tmpList.pt = std::abs(1.0f / m_Par[iP].At(itrack, 3, 0));
1267               tmpList.chi2 = m_Chi2(itrack, 0, 0) + chi2;
1268               tmpList.chi2_hit = chi2;
1269               tmpList.score = getScoreStruct(tmpList);
1270               cloner.add_cand(m_SeedIdx(itrack, 0, 0) - offset, tmpList);
1271 
1272               dprint("  adding hit with hit_cnt=" << hit_cnt << " for trkIdx=" << tmpList.trkIdx
1273                                                   << " orig Seed=" << m_Label(itrack, 0, 0));
1274             }
1275           }
1276         }
1277       }
1278 
1279     }  //end loop over hits
1280 
1281     //now add invalid hit
1282     for (int itrack = 0; itrack < N_proc; ++itrack) {
1283       dprint("num_all_minus_one_hits(" << itrack << ")=" << num_all_minus_one_hits(itrack));
1284 
1285       // Cands that miss the layer are stashed away in MkBuilder(), before propagation,
1286       // and then merged back afterwards.
1287       if (m_XWsrResult[itrack].m_wsr == WSR_Outside) {
1288         continue;
1289       }
1290 
1291       // int fake_hit_idx = num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand ? -1 : -2;
1292       int fake_hit_idx = ((num_all_minus_one_hits(itrack) < m_iteration_params->maxHolesPerCand) &&
1293                           (m_NTailMinusOneHits(itrack, 0, 0) < m_iteration_params->maxConsecHoles))
1294                              ? Hit::kHitMissIdx
1295                              : Hit::kHitStopIdx;
1296 
1297       if (m_XWsrResult[itrack].m_wsr == WSR_Edge) {
1298         fake_hit_idx = Hit::kHitEdgeIdx;
1299       }
1300       //now add fake hit for tracks that passsed through inactive modules
1301       else if (m_XWsrResult[itrack].m_in_gap == true && nHitsAdded[itrack] == 0) {
1302         fake_hit_idx = Hit::kHitInGapIdx;
1303       }
1304       //now add fake hit for cases where hit cluster size is larger than maxClusterSize
1305       else if (isTooLargeCluster[itrack] == true && nHitsAdded[itrack] == 0) {
1306         fake_hit_idx = Hit::kHitMaxClusterIdx;
1307       }
1308 
1309       IdxChi2List tmpList;
1310       tmpList.trkIdx = m_CandIdx(itrack, 0, 0);
1311       tmpList.hitIdx = fake_hit_idx;
1312       tmpList.module = -1;
1313       tmpList.nhits = m_NFoundHits(itrack, 0, 0);
1314       tmpList.ntailholes = (fake_hit_idx == Hit::kHitMissIdx ? m_NTailMinusOneHits(itrack, 0, 0) + 1
1315                                                              : m_NTailMinusOneHits(itrack, 0, 0));
1316       tmpList.noverlaps = m_NOverlapHits(itrack, 0, 0);
1317       tmpList.nholes = num_inside_minus_one_hits(itrack);
1318       tmpList.pt = std::abs(1.0f / m_Par[iP].At(itrack, 3, 0));
1319       tmpList.chi2 = m_Chi2(itrack, 0, 0);
1320       tmpList.chi2_hit = 0;
1321       tmpList.score = getScoreStruct(tmpList);
1322       cloner.add_cand(m_SeedIdx(itrack, 0, 0) - offset, tmpList);
1323       dprint("adding invalid hit " << fake_hit_idx);
1324     }
1325   }
1326 
1327   //==============================================================================
1328   // UpdateWithLastHit
1329   //==============================================================================
1330 
1331   void MkFinder::updateWithLoadedHit(int N_proc, const FindingFoos &fnd_foos) {
1332     // See comment in MkBuilder::find_tracks_in_layer() about intra / inter flags used here
1333     // for propagation to the hit.
1334     (*fnd_foos.m_update_param_foo)(m_Err[iP],
1335                                    m_Par[iP],
1336                                    m_Chg,
1337                                    m_msErr,
1338                                    m_msPar,
1339                                    m_Err[iC],
1340                                    m_Par[iC],
1341                                    N_proc,
1342                                    m_prop_config->finding_inter_layer_pflags,
1343                                    m_prop_config->finding_requires_propagation_to_hit_pos);
1344   }
1345 
1346   //==============================================================================
1347   // CopyOutParErr
1348   //==============================================================================
1349 
1350   void MkFinder::copyOutParErr(std::vector<CombCandidate> &seed_cand_vec, int N_proc, bool outputProp) const {
1351     const int iO = outputProp ? iP : iC;
1352 
1353     for (int i = 0; i < N_proc; ++i) {
1354       TrackCand &cand = seed_cand_vec[m_SeedIdx(i, 0, 0)][m_CandIdx(i, 0, 0)];
1355 
1356       // Set the track state to the updated parameters
1357       m_Err[iO].copyOut(i, cand.errors_nc().Array());
1358       m_Par[iO].copyOut(i, cand.parameters_nc().Array());
1359       cand.setCharge(m_Chg(i, 0, 0));
1360 
1361       dprint((outputProp ? "propagated" : "updated")
1362              << " track parameters x=" << cand.parameters()[0] << " y=" << cand.parameters()[1]
1363              << " z=" << cand.parameters()[2] << " pt=" << 1. / cand.parameters()[3] << " posEta=" << cand.posEta());
1364     }
1365   }
1366 
1367   //==============================================================================
1368   // Backward Fit hack
1369   //==============================================================================
1370 
1371   void MkFinder::bkFitInputTracks(TrackVec &cands, int beg, int end) {
1372     // Uses HitOnTrack vector from Track directly + a local cursor array to current hit.
1373 
1374     MatriplexTrackPacker mtp(&cands[beg]);
1375 
1376     int itrack = 0;
1377 
1378     for (int i = beg; i < end; ++i, ++itrack) {
1379       const Track &trk = cands[i];
1380 
1381       m_Chg(itrack, 0, 0) = trk.charge();
1382       m_CurHit[itrack] = trk.nTotalHits() - 1;
1383       m_HoTArr[itrack] = trk.getHitsOnTrackArray();
1384 
1385       mtp.addInput(trk);
1386     }
1387 
1388     m_Chi2.setVal(0);
1389 
1390     mtp.pack(m_Err[iC], m_Par[iC]);
1391 
1392     m_Err[iC].scale(100.0f);
1393   }
1394 
1395   void MkFinder::bkFitInputTracks(EventOfCombCandidates &eocss, int beg, int end) {
1396     // Could as well use HotArrays from tracks directly + a local cursor array to last hit.
1397 
1398     // XXXX - shall we assume only TrackCand-zero is needed and that we can freely
1399     // bork the HoTNode array?
1400 
1401     MatriplexTrackPacker mtp(&eocss[beg][0]);
1402 
1403     int itrack = 0;
1404 
1405     for (int i = beg; i < end; ++i, ++itrack) {
1406       const TrackCand &trk = eocss[i][0];
1407 
1408       m_Chg(itrack, 0, 0) = trk.charge();
1409       m_CurNode[itrack] = trk.lastCcIndex();
1410       m_HoTNodeArr[itrack] = trk.combCandidate()->hotsData();
1411 
1412       // XXXX Need TrackCand* to update num-hits. Unless I collect info elsewhere
1413       // and fix it in BkFitOutputTracks.
1414       m_TrkCand[itrack] = &eocss[i][0];
1415 
1416       mtp.addInput(trk);
1417     }
1418 
1419     m_Chi2.setVal(0);
1420 
1421     mtp.pack(m_Err[iC], m_Par[iC]);
1422 
1423     m_Err[iC].scale(100.0f);
1424   }
1425 
1426   //------------------------------------------------------------------------------
1427 
1428   void MkFinder::bkFitOutputTracks(TrackVec &cands, int beg, int end, bool outputProp) {
1429     // Only copy out track params / errors / chi2, all the rest is ok.
1430 
1431     const int iO = outputProp ? iP : iC;
1432 
1433     int itrack = 0;
1434     for (int i = beg; i < end; ++i, ++itrack) {
1435       Track &trk = cands[i];
1436 
1437       m_Err[iO].copyOut(itrack, trk.errors_nc().Array());
1438       m_Par[iO].copyOut(itrack, trk.parameters_nc().Array());
1439 
1440       trk.setChi2(m_Chi2(itrack, 0, 0));
1441       if (isFinite(trk.chi2())) {
1442         trk.setScore(getScoreCand(trk));
1443       }
1444     }
1445   }
1446 
1447   void MkFinder::bkFitOutputTracks(EventOfCombCandidates &eocss, int beg, int end, bool outputProp) {
1448     // Only copy out track params / errors / chi2, all the rest is ok.
1449 
1450     // XXXX - where will rejected hits get removed?
1451 
1452     const int iO = outputProp ? iP : iC;
1453 
1454     int itrack = 0;
1455     for (int i = beg; i < end; ++i, ++itrack) {
1456       TrackCand &trk = eocss[i][0];
1457 
1458       m_Err[iO].copyOut(itrack, trk.errors_nc().Array());
1459       m_Par[iO].copyOut(itrack, trk.parameters_nc().Array());
1460 
1461       trk.setChi2(m_Chi2(itrack, 0, 0));
1462       if (isFinite(trk.chi2())) {
1463         trk.setScore(getScoreCand(trk));
1464       }
1465     }
1466   }
1467 
1468   //------------------------------------------------------------------------------
1469 
1470 #if defined(DEBUG_BACKWARD_FIT) || defined(DEBUG_BACKWARD_FIT_BH)
1471   namespace {
1472     float e2s(float x) { return 1e4 * std::sqrt(x); }
1473   }  // namespace
1474 #endif
1475 
1476   void MkFinder::bkFitFitTracksBH(const EventOfHits &eventofhits,
1477                                   const SteeringParams &st_par,
1478                                   const int N_proc,
1479                                   bool chiDebug) {
1480     // Prototyping final backward fit.
1481     // This works with track-finding indices, before remapping.
1482     //
1483     // Layers should be collected during track finding and list all layers that have actual hits.
1484     // Then we could avoid checking which layers actually do have hits.
1485 
1486     MPlexQF tmp_chi2;
1487     float tmp_err[6] = {666, 0, 666, 0, 0, 666};
1488     float tmp_pos[3];
1489 
1490     for (auto lp_iter = st_par.m_layer_plan.rbegin(); lp_iter != st_par.m_layer_plan.rend(); ++lp_iter) {
1491       const int layer = lp_iter->m_layer;
1492 
1493       const LayerOfHits &L = eventofhits[layer];
1494       const LayerInfo &LI = *L.layer_info();
1495 
1496       int count = 0;
1497       for (int i = 0; i < N_proc; ++i) {
1498         while (m_CurHit[i] >= 0 && m_HoTArr[i][m_CurHit[i]].index < 0)
1499           --m_CurHit[i];
1500 
1501         if (m_CurHit[i] >= 0 && m_HoTArr[i][m_CurHit[i]].layer == layer) {
1502           // Skip the overlap hits -- if they exist.
1503           // 1. Overlap hit gets placed *after* the original hit in TrackCand::exportTrack()
1504           // which is *before* in the reverse iteration that we are doing here.
1505           // 2. Seed-hit merging can result in more than two hits per layer.
1506           while (m_CurHit[i] > 0 && m_HoTArr[i][m_CurHit[i] - 1].layer == layer)
1507             --m_CurHit[i];
1508 
1509           const Hit &hit = L.refHit(m_HoTArr[i][m_CurHit[i]].index);
1510           m_msErr.copyIn(i, hit.errArray());
1511           m_msPar.copyIn(i, hit.posArray());
1512           ++count;
1513           --m_CurHit[i];
1514         } else {
1515           tmp_pos[0] = m_Par[iC](i, 0, 0);
1516           tmp_pos[1] = m_Par[iC](i, 1, 0);
1517           tmp_pos[2] = m_Par[iC](i, 2, 0);
1518           m_msErr.copyIn(i, tmp_err);
1519           m_msPar.copyIn(i, tmp_pos);
1520         }
1521       }
1522 
1523       if (count == 0)
1524         continue;
1525 
1526       // ZZZ Could add missing hits here, only if there are any actual matches.
1527 
1528       if (LI.is_barrel()) {
1529         propagateTracksToHitR(m_msPar, N_proc, m_prop_config->backward_fit_pflags);
1530 
1531         kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
1532                         m_Err[iP],
1533                         m_Par[iP],
1534                         m_msErr,
1535                         m_msPar,
1536                         m_Err[iC],
1537                         m_Par[iC],
1538                         tmp_chi2,
1539                         N_proc);
1540       } else {
1541         propagateTracksToHitZ(m_msPar, N_proc, m_prop_config->backward_fit_pflags);
1542 
1543         kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params,
1544                               m_Err[iP],
1545                               m_Par[iP],
1546                               m_msErr,
1547                               m_msPar,
1548                               m_Err[iC],
1549                               m_Par[iC],
1550                               tmp_chi2,
1551                               N_proc);
1552       }
1553 
1554       //fixup invpt sign and charge
1555       for (int n = 0; n < N_proc; ++n) {
1556         if (m_Par[iC].At(n, 3, 0) < 0) {
1557           m_Chg.At(n, 0, 0) = -m_Chg.At(n, 0, 0);
1558           m_Par[iC].At(n, 3, 0) = -m_Par[iC].At(n, 3, 0);
1559         }
1560       }
1561 
1562 #ifdef DEBUG_BACKWARD_FIT_BH
1563       // Dump per hit chi2
1564       for (int i = 0; i < N_proc; ++i) {
1565         float r_h = std::hypot(m_msPar.At(i, 0, 0), m_msPar.At(i, 1, 0));
1566         float r_t = std::hypot(m_Par[iC].At(i, 0, 0), m_Par[iC].At(i, 1, 0));
1567 
1568         // if ((std::isnan(tmp_chi2[i]) || std::isnan(r_t)))
1569         // if ( ! std::isnan(tmp_chi2[i]) && tmp_chi2[i] > 0) // && tmp_chi2[i] > 30)
1570         if (chiDebug) {
1571           int ti = iP;
1572           printf(
1573               "CHIHIT %3d %10g %10g %10g %10g %10g %11.5g %11.5g %11.5g %10g %10g %10g %10g %11.5g %11.5g %11.5g %10g "
1574               "%10g %10g %10g %10g %11.5g %11.5g\n",
1575               layer,
1576               tmp_chi2[i],
1577               m_msPar.At(i, 0, 0),
1578               m_msPar.At(i, 1, 0),
1579               m_msPar.At(i, 2, 0),
1580               r_h,  // x_h y_h z_h r_h -- hit pos
1581               e2s(m_msErr.At(i, 0, 0)),
1582               e2s(m_msErr.At(i, 1, 1)),
1583               e2s(m_msErr.At(i, 2, 2)),  // ex_h ey_h ez_h -- hit errors
1584               m_Par[ti].At(i, 0, 0),
1585               m_Par[ti].At(i, 1, 0),
1586               m_Par[ti].At(i, 2, 0),
1587               r_t,  // x_t y_t z_t r_t -- track pos
1588               e2s(m_Err[ti].At(i, 0, 0)),
1589               e2s(m_Err[ti].At(i, 1, 1)),
1590               e2s(m_Err[ti].At(i, 2, 2)),  // ex_t ey_t ez_t -- track errors
1591               1.0f / m_Par[ti].At(i, 3, 0),
1592               m_Par[ti].At(i, 4, 0),
1593               m_Par[ti].At(i, 5, 0),                                     // pt, phi, theta
1594               std::atan2(m_msPar.At(i, 1, 0), m_msPar.At(i, 0, 0)),      // phi_h
1595               std::atan2(m_Par[ti].At(i, 1, 0), m_Par[ti].At(i, 0, 0)),  // phi_t
1596               1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0),
1597                                 m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)),  // d_xy
1598               1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0))             // d_z
1599               // e2s((m_msErr.At(i,0,0) + m_msErr.At(i,1,1)) / (r_h * r_h)),     // ephi_h
1600               // e2s((m_Err[ti].At(i,0,0) + m_Err[ti].At(i,1,1)) / (r_t * r_t))  // ephi_t
1601           );
1602         }
1603       }
1604 #endif
1605 
1606       // update chi2
1607       m_Chi2.add(tmp_chi2);
1608     }
1609   }
1610 
1611   //------------------------------------------------------------------------------
1612 
1613   void MkFinder::bkFitFitTracks(const EventOfHits &eventofhits,
1614                                 const SteeringParams &st_par,
1615                                 const int N_proc,
1616                                 bool chiDebug) {
1617     // Prototyping final backward fit.
1618     // This works with track-finding indices, before remapping.
1619     //
1620     // Layers should be collected during track finding and list all layers that have actual hits.
1621     // Then we could avoid checking which layers actually do have hits.
1622 
1623     MPlexQF tmp_chi2;
1624     MPlexQI no_mat_effs;
1625     float tmp_err[6] = {666, 0, 666, 0, 0, 666};
1626     float tmp_pos[3];
1627 
1628     for (auto lp_iter = st_par.make_iterator(SteeringParams::IT_BkwFit); lp_iter.is_valid(); ++lp_iter) {
1629       const int layer = lp_iter.layer();
1630 
1631       const LayerOfHits &L = eventofhits[layer];
1632       const LayerInfo &LI = *L.layer_info();
1633 
1634       // XXXX
1635 #if defined(DEBUG_BACKWARD_FIT)
1636       const Hit *last_hit_ptr[NN];
1637 #endif
1638 
1639       no_mat_effs.setVal(0);
1640       int done_count = 0;
1641       int here_count = 0;
1642       for (int i = 0; i < N_proc; ++i) {
1643         while (m_CurNode[i] >= 0 && m_HoTNodeArr[i][m_CurNode[i]].m_hot.index < 0) {
1644           m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
1645         }
1646 
1647         if (m_CurNode[i] < 0)
1648           ++done_count;
1649 
1650         if (m_CurNode[i] >= 0 && m_HoTNodeArr[i][m_CurNode[i]].m_hot.layer == layer) {
1651           // Skip the overlap hits -- if they exist.
1652           // 1. Overlap hit gets placed *after* the original hit in TrackCand::exportTrack()
1653           // which is *before* in the reverse iteration that we are doing here.
1654           // 2. Seed-hit merging can result in more than two hits per layer.
1655           // while (m_CurHit[i] > 0 && m_HoTArr[ i ][ m_CurHit[i] - 1 ].layer == layer) --m_CurHit[i];
1656           while (m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx >= 0 &&
1657                  m_HoTNodeArr[i][m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx].m_hot.layer == layer)
1658             m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
1659 
1660           const Hit &hit = L.refHit(m_HoTNodeArr[i][m_CurNode[i]].m_hot.index);
1661 
1662 #ifdef DEBUG_BACKWARD_FIT
1663           last_hit_ptr[i] = &hit;
1664 #endif
1665           m_msErr.copyIn(i, hit.errArray());
1666           m_msPar.copyIn(i, hit.posArray());
1667           ++here_count;
1668 
1669           m_CurNode[i] = m_HoTNodeArr[i][m_CurNode[i]].m_prev_idx;
1670         } else {
1671 #ifdef DEBUG_BACKWARD_FIT
1672           last_hit_ptr[i] = nullptr;
1673 #endif
1674           no_mat_effs[i] = 1;
1675           tmp_pos[0] = m_Par[iC](i, 0, 0);
1676           tmp_pos[1] = m_Par[iC](i, 1, 0);
1677           tmp_pos[2] = m_Par[iC](i, 2, 0);
1678           m_msErr.copyIn(i, tmp_err);
1679           m_msPar.copyIn(i, tmp_pos);
1680         }
1681       }
1682 
1683       if (done_count == N_proc)
1684         break;
1685       if (here_count == 0)
1686         continue;
1687 
1688       // ZZZ Could add missing hits here, only if there are any actual matches.
1689 
1690       if (LI.is_barrel()) {
1691         propagateTracksToHitR(m_msPar, N_proc, m_prop_config->backward_fit_pflags, &no_mat_effs);
1692 
1693         kalmanOperation(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
1694                         m_Err[iP],
1695                         m_Par[iP],
1696                         m_msErr,
1697                         m_msPar,
1698                         m_Err[iC],
1699                         m_Par[iC],
1700                         tmp_chi2,
1701                         N_proc);
1702       } else {
1703         propagateTracksToHitZ(m_msPar, N_proc, m_prop_config->backward_fit_pflags, &no_mat_effs);
1704 
1705         kalmanOperationEndcap(KFO_Calculate_Chi2 | KFO_Update_Params,
1706                               m_Err[iP],
1707                               m_Par[iP],
1708                               m_msErr,
1709                               m_msPar,
1710                               m_Err[iC],
1711                               m_Par[iC],
1712                               tmp_chi2,
1713                               N_proc);
1714       }
1715 
1716       //fixup invpt sign and charge
1717       for (int n = 0; n < N_proc; ++n) {
1718         if (m_Par[iC].At(n, 3, 0) < 0) {
1719           m_Chg.At(n, 0, 0) = -m_Chg.At(n, 0, 0);
1720           m_Par[iC].At(n, 3, 0) = -m_Par[iC].At(n, 3, 0);
1721         }
1722       }
1723 
1724       for (int i = 0; i < N_proc; ++i) {
1725 #if defined(DEBUG_BACKWARD_FIT)
1726         if (chiDebug && last_hit_ptr[i]) {
1727           TrackCand &bb = *m_TrkCand[i];
1728           int ti = iP;
1729           float chi = tmp_chi2.At(i, 0, 0);
1730           float chi_prnt = std::isfinite(chi) ? chi : -9;
1731 
1732 #if defined(MKFIT_STANDALONE)
1733           const MCHitInfo &mchi = m_event->simHitsInfo_[last_hit_ptr[i]->mcHitID()];
1734 
1735           printf(
1736               "BKF_OVERLAP %d %d %d %d %d %d %d "
1737               "%f %f %f %f %d %d %d %d "
1738               "%f %f %f %f %f\n",
1739               m_event->evtID(),
1740 #else
1741           printf(
1742               "BKF_OVERLAP %d %d %d %d %d %d "
1743               "%f %f %f %f %d %d %d "
1744               "%f %f %f %f %f\n",
1745 #endif
1746               bb.label(),
1747               (int)bb.prodType(),
1748               bb.isFindable(),
1749               layer,
1750               L.is_stereo(),
1751               L.is_barrel(),
1752               bb.pT(),
1753               bb.posEta(),
1754               bb.posPhi(),
1755               chi_prnt,
1756               std::isnan(chi),
1757               std::isfinite(chi),
1758               chi > 0,
1759 #if defined(MKFIT_STANDALONE)
1760               mchi.mcTrackID(),
1761 #endif
1762               e2s(m_Err[ti].At(i, 0, 0)),
1763               e2s(m_Err[ti].At(i, 1, 1)),
1764               e2s(m_Err[ti].At(i, 2, 2)),  // sx_t sy_t sz_t -- track errors
1765               1e4f * std::hypot(m_msPar.At(i, 0, 0) - m_Par[ti].At(i, 0, 0),
1766                                 m_msPar.At(i, 1, 0) - m_Par[ti].At(i, 1, 0)),  // d_xy
1767               1e4f * (m_msPar.At(i, 2, 0) - m_Par[ti].At(i, 2, 0))             // d_z
1768           );
1769         }
1770 #endif
1771       }
1772 
1773       // update chi2
1774       m_Chi2.add(tmp_chi2);
1775     }
1776   }
1777 
1778   //------------------------------------------------------------------------------
1779 
1780   void MkFinder::bkFitPropTracksToPCA(const int N_proc) {
1781     propagateTracksToPCAZ(N_proc, m_prop_config->pca_prop_pflags);
1782   }
1783 
1784 }  // end namespace mkfit