Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-23 02:42:47

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