Draft PR for prototype final fit

RecoTracker/MkFitCMS/standalone/buildtestMPlex.cc

@@ -562,6 +562,12 @@ namespace mkfit {
 
         builder.export_tracks(ev.fitTracks_);
       }
+
+      ////main call to fit - this is where you can measure the time as for findTracksCloneEngine
+      //builder.fittracks();
+      ////export tracks as "fit" tracks instead of the ones coming from backward fit (should remove the export above)
+      //builder.export_tracks(ev.fitTracks_);
+
       ev.resetCurrentSeedTracks();
 
       builder.end_event();

RecoTracker/MkFitCore/interface/MkBuilder.h

@@ -116,6 +116,12 @@ namespace mkfit {
     void backwardFit();
     void fit_cands(MkFinder *mkfndr, int start_cand, int end_cand, int region);
 
+    //refit
+
+    void fittracks();
+    void fit_tracks(MkFinder *mkfndr, int nFoundTracks, std::vector<int> inds, int start_trk, int end_trk);
+    void check_tracks(std::vector<int> inds, int start_trk, int end_trk);
+
   private:
     void fit_one_seed_set(TrackVec &simtracks, int itrack, int end, MkFitter *mkfttr, const bool is_brl[]);
 

RecoTracker/MkFitCore/src/KalmanUtilsMPlex.cc

@@ -1255,7 +1255,7 @@ namespace mkfit {
     if (propToHit) {
       MPlexLS propErr;
       MPlexLV propPar;
-      propagateHelixToPlaneMPlex(psErr, psPar, Chg, msPar, plNrm, propErr, propPar, outFailFlag, N_proc, propFlags);
+      propagateHelixToPlaneMPlex(psErr, psPar, Chg, plPnt, plNrm, propErr, propPar, outFailFlag, N_proc, propFlags);
 
       kalmanOperationPlaneLocal(KFO_Update_Params | KFO_Local_Cov,
                                 propErr,
@@ -1295,6 +1295,68 @@ namespace mkfit {
 
   //------------------------------------------------------------------------------
 
+  void kalmanPropagateAndUpdateAndChi2Plane(const MPlexLS& psErr,
+                                            const MPlexLV& psPar,
+                                            MPlexQI& Chg,
+                                            const MPlexHS& msErr,
+                                            const MPlexHV& msPar,
+                                            const MPlexHV& plNrm,
+                                            const MPlexHV& plDir,
+                                            const MPlexHV& plPnt,
+                                            MPlexLS& outErr,
+                                            MPlexLV& outPar,
+                                            MPlexQI& outFailFlag,
+                                            MPlexQF& outChi2,
+                                            const int N_proc,
+                                            const PropagationFlags& propFlags,
+                                            const bool propToHit,
+                                            const MPlexQI* noMatEffPtr) {
+    if (propToHit) {
+      MPlexLS propErr;
+      MPlexLV propPar;
+
+      propagateHelixToPlaneMPlex(
+          psErr, psPar, Chg, plPnt, plNrm, propErr, propPar, outFailFlag, N_proc, propFlags, noMatEffPtr);
+
+      kalmanOperationPlaneLocal(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
+                                propErr,
+                                propPar,
+                                Chg,
+                                msErr,
+                                msPar,
+                                plNrm,
+                                plDir,
+                                plPnt,
+                                outErr,
+                                outPar,
+                                outChi2,
+                                N_proc);
+
+    } else {
+      kalmanOperationPlaneLocal(KFO_Calculate_Chi2 | KFO_Update_Params | KFO_Local_Cov,
+                                psErr,
+                                psPar,
+                                Chg,
+                                msErr,
+                                msPar,
+                                plNrm,
+                                plDir,
+                                plPnt,
+                                outErr,
+                                outPar,
+                                outChi2,
+                                N_proc);
+    }
+    for (int n = 0; n < NN; ++n) {
+      if (outPar.At(n, 3, 0) < 0) {
+        Chg.At(n, 0, 0) = -Chg.At(n, 0, 0);
+        outPar.At(n, 3, 0) = -outPar.At(n, 3, 0);
+      }
+    }
+  }
+
+  //------------------------------------------------------------------------------
+
   void kalmanComputeChi2Plane(const MPlexLS& psErr,
                               const MPlexLV& psPar,
                               const MPlexQI& inChg,
@@ -1337,7 +1399,7 @@ namespace mkfit {
     propPar = psPar;
     if (propToHit) {
       MPlexLS propErr;
-      propagateHelixToPlaneMPlex(psErr, psPar, inChg, msPar, plNrm, propErr, propPar, outFailFlag, N_proc, propFlags);
+      propagateHelixToPlaneMPlex(psErr, psPar, inChg, plPnt, plNrm, propErr, propPar, outFailFlag, N_proc, propFlags);
 
       kalmanOperationPlaneLocal(KFO_Calculate_Chi2,
                                 propErr,

RecoTracker/MkFitCore/src/KalmanUtilsMPlex.h

@@ -143,6 +143,23 @@ namespace mkfit {
                                      const PropagationFlags& propFlags,
                                      const bool propToHit);
 
+  void kalmanPropagateAndUpdateAndChi2Plane(const MPlexLS& psErr,
+                                            const MPlexLV& psPar,
+                                            MPlexQI& Chg,
+                                            const MPlexHS& msErr,
+                                            const MPlexHV& msPar,
+                                            const MPlexHV& plNrm,
+                                            const MPlexHV& plDir,
+                                            const MPlexHV& plPnt,
+                                            MPlexLS& outErr,
+                                            MPlexLV& outPar,
+                                            MPlexQI& outFailFlag,
+                                            MPlexQF& outChi2,
+                                            const int N_proc,
+                                            const PropagationFlags& propFlags,
+                                            const bool propToHit,
+                                            const MPlexQI* noMatEffPtr = nullptr);
+
   void kalmanComputeChi2Plane(const MPlexLS& psErr,
                               const MPlexLV& psPar,
                               const MPlexQI& inChg,

RecoTracker/MkFitCore/src/MkBuilder.cc

@@ -19,6 +19,7 @@
 #endif
 
 //#define DEBUG
+//#define DEBUG_FIT
 #include "Debug.h"
 //#define DEBUG_FINAL_FIT
 
@@ -1404,4 +1405,118 @@ namespace mkfit {
     mkfndr->release();
   }
 
+  //==============================================================================
+  // ReFit
+  //==============================================================================
+
+  void MkBuilder::fit_tracks(MkFinder *mkfndr, int nFoundHits, std::vector<int> inds, int start_trk, int end_trk) {
+    //could be wrapped into some setup_fit
+    const TrackerInfo &ti = m_job->m_trk_info;
+    PropagationFlags my_flags = PropagationFlags(PF_use_param_b_field | PF_apply_material);
+    my_flags.tracker_info = &ti;
+    mkfndr->refit_flags = &my_flags;
+
+    mkfndr->m_event = m_event;
+
+    for (int icand = start_trk; icand < end_trk; icand += NN) {
+      const int end = std::min(icand + NN, end_trk);
+
+      // input candidate tracks
+      mkfndr->fwdFitInputTracks(m_tracks, inds, icand, end);
+
+      //prepare indices
+      std::vector<std::vector<int>> indices_R2 = mkfndr->reFitIndices(m_job->m_event_of_hits, end - icand, nFoundHits);
+
+      // fit the tracks from the input in fwd direction
+      mkfndr->fwdFitFitTracks(m_job->m_event_of_hits, end - icand, nFoundHits, indices_R2);
+
+      // fwdFitOutput
+      mkfndr->reFitOutputTracks(m_tracks, inds, icand, end, nFoundHits);
+
+      // input candidate tracks
+      mkfndr->bkReFitInputTracks(m_tracks, inds, icand, end);
+
+      // fit the tracks from the input in bkw direction
+      mkfndr->bkReFitFitTracks(m_job->m_event_of_hits, end - icand, nFoundHits, indices_R2);
+
+      // fwdFitOutput
+      mkfndr->reFitOutputTracks(m_tracks, inds, icand, end, nFoundHits, true);
+    }
+
+    mkfndr->release();
+  }
+
+  void MkBuilder::check_tracks(std::vector<int> inds, int start_trk, int end_trk) {
+    for (int icand = start_trk; icand < end_trk; icand += NN) {
+      const int end = std::min(icand + NN, end_trk);
+
+      std::cout << "BEGIN CHECKs" << std::endl;
+      std::cout << " strt " << start_trk << " end " << end_trk << " end " << end_trk - start_trk << " NN " << NN
+                << std::endl;
+
+      for (int i = start_trk; i < end; ++i) {
+        const Track &trk = m_tracks[inds[i]];
+
+        std::cout << "trk pt " << trk.pT() << " trk eta " << trk.momEta() << " trk phi " << trk.momPhi() << std::endl;
+        std::cout << "trk nTotalHits " << trk.nTotalHits() << " trk nFoundHits " << trk.nFoundHits() << std::endl;
+        std::cout << "END CHECK" << std::endl;
+      }
+    }
+  }
+
+  void MkBuilder::fittracks() {
+#ifdef DEBUG_FIT
+    std::cout << "here are N tracks " << m_tracks.size() << std::endl;
+#endif
+    int N = 0;
+
+    std::map<int, std::vector<int>> mapFoundHits;
+    for (auto &t : m_tracks) {
+#ifdef DEBUG_FIT
+      std::cout << "______________ " << std::endl;
+      std::cout << "track N " << N << " nhits " << t.nFoundHits() << std::endl;
+      std::cout << "trk pt " << t.pT() << " trk eta " << t.momEta() << std::endl;
+      std::cout << "trk nTotalHits " << t.nTotalHits() << " trk nFoundHits " << t.nFoundHits() << std::endl;
+      std::cout << "trk last l " << t.getLastFoundHitLyr() << " trk last idx " << t.getLastFoundHitIdx() << std::endl;
+      for (int i = 0; i < t.nTotalHits(); i++) {
+        std::cout << "index " << t.getHitIdx(i) << " layer " << t.getHitLyr(i) << std::endl;
+      }
+#endif
+      int foundh = t.nFoundHits();
+      if (mapFoundHits.count(foundh)) {
+        mapFoundHits[foundh].push_back(N);
+      } else
+        mapFoundHits[foundh] = {N};
+      N++;
+    }
+#ifdef DEBUG_FIT
+    for (auto &m : mapFoundHits) {
+      std::cout << m.first << " the key " << std::endl;
+      for (auto i : m.second) {
+        std::cout << i << " index ";
+      }
+      std::cout << "\n";
+    }
+#endif
+    auto mkfndr = g_exe_ctx.m_finders.makeOrGet();
+    for (auto &m : mapFoundHits) {
+      int ntimes = m.second.size() / NN;
+#ifdef DEBUG_FIT
+      std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
+      std::cout << "ntimes " << ntimes << " size  " << m.second.size() << " extra " << m.second.size() - NN * ntimes
+                << std::endl;
+#endif
+      for (int i = 0; i < ntimes; i++) {
+#ifdef DEBUG_FIT
+        check_tracks(m.second, NN * i, NN * (i + 1));
+#endif
+        fit_tracks(mkfndr.get(), m.first, m.second, NN * i, NN * (i + 1));
+      }
+#ifdef DEBUG_FIT
+      check_tracks(m.second, NN * ntimes, m.second.size());
+#endif
+      fit_tracks(mkfndr.get(), m.first, m.second, NN * ntimes, m.second.size());
+    }
+  }
+
 }  // end namespace mkfit