Project CMSSW displayed by LXR CMSSW_15_1_X_2025-02-16-2300/​RecoTracker/​MkFitCore/​src/​MkFitter.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_15_1_X_2025-02-16-2300 CMSSW_15_1_X_2025-02-14-2300 CMSSW_15_1_X_2025-02-13-2300 CMSSW_15_1_X_2025-02-12-2300 CMSSW_15_1_X_2025-02-11-2300 CMSSW_15_1_X_2025-02-10-2300 Revert   Change

File indexing completed on 2025-01-18 03:42:20

 #include "MkFitter.h"
 
 #include "KalmanUtilsMPlex.h"
 #include "MatriplexPackers.h"
 
 //#define DEBUG
 #include "Debug.h"
 
 #include <sstream>
 
 namespace mkfit {
 
   void MkFitter::checkAlignment() {
     printf("MkFitter alignment check:\n");
     Matriplex::align_check("  m_Err[0]   =", &m_Err[0].fArray[0]);
     Matriplex::align_check("  m_Err[1]   =", &m_Err[1].fArray[0]);
     Matriplex::align_check("  m_Par[0]   =", &m_Par[0].fArray[0]);
     Matriplex::align_check("  m_Par[1]   =", &m_Par[1].fArray[0]);
     Matriplex::align_check("  m_msErr[0] =", &m_msErr[0].fArray[0]);
     Matriplex::align_check("  m_msPar[0] =", &m_msPar[0].fArray[0]);
   }
 
   void MkFitter::printPt(int idx) {
     for (int i = 0; i < NN; ++i) {
       printf("%5.2f  ", getHypot(m_Par[idx].At(i, 3, 0), m_Par[idx].At(i, 4, 0)));
     }
   }
 
   int MkFitter::countValidHits(int itrack, int end_hit) const {
     int result = 0;
     for (int hi = 0; hi < end_hit; ++hi) {
       if (m_HoTArr[hi](itrack, 0, 0).index >= 0)
         result++;
     }
     return result;
   }
 
   int MkFitter::countInvalidHits(int itrack, int end_hit) const {
     int result = 0;
     for (int hi = 0; hi < end_hit; ++hi) {
       // XXXX MT: Should also count -2 hits as invalid?
       if (m_HoTArr[hi](itrack, 0, 0).index == -1)
         result++;
     }
     return result;
   }
 
   //==============================================================================
 
   void MkFitter::inputTracksAndHits(const std::vector<Track>& tracks,
                                     const std::vector<HitVec>& layerHits,
                                     int beg,
                                     int end) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     int itrack = 0;
 
     for (int i = beg; i < end; ++i, ++itrack) {
       const Track& trk = tracks[i];
 
       m_Err[iC].copyIn(itrack, trk.errors().Array());
       m_Par[iC].copyIn(itrack, trk.parameters().Array());
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
       m_Label(itrack, 0, 0) = trk.label();
 
       // CopyIn seems fast enough, but indirections are quite slow.
       for (int hi = 0; hi < m_Nhits; ++hi) {
         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
 
         const int hidx = trk.getHitIdx(hi);
         if (hidx < 0)
           continue;
 
         const Hit& hit = layerHits[hi][hidx];
         m_msErr[hi].copyIn(itrack, hit.errArray());
         m_msPar[hi].copyIn(itrack, hit.posArray());
       }
     }
   }
 
   void MkFitter::inputTracksAndHits(const std::vector<Track>& tracks,
                                     const std::vector<LayerOfHits>& layerHits,
                                     int beg,
                                     int end) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     int itrack;
 
     for (int i = beg; i < end; ++i) {
       itrack = i - beg;
       const Track& trk = tracks[i];
 
       m_Label(itrack, 0, 0) = trk.label();
 
       m_Err[iC].copyIn(itrack, trk.errors().Array());
       m_Par[iC].copyIn(itrack, trk.parameters().Array());
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
 
       // CopyIn seems fast enough, but indirections are quite slow.
       for (int hi = 0; hi < m_Nhits; ++hi) {
         const int hidx = trk.getHitIdx(hi);
         const int hlyr = trk.getHitLyr(hi);
         const Hit& hit = layerHits[hlyr].refHit(hidx);
 
         m_msErr[hi].copyIn(itrack, hit.errArray());
         m_msPar[hi].copyIn(itrack, hit.posArray());
 
         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
       }
     }
   }
 
   void MkFitter::slurpInTracksAndHits(const std::vector<Track>& tracks,
                                       const std::vector<HitVec>& layerHits,
                                       int beg,
                                       int end) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     MatriplexTrackPacker mtp(&tracks[beg]);
 
     for (int i = beg; i < end; ++i) {
       int itrack = i - beg;
       const Track& trk = tracks[i];
 
       m_Label(itrack, 0, 0) = trk.label();
 
       mtp.addInput(trk);
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
     }
 
     mtp.pack(m_Err[iC], m_Par[iC]);
 
     // CopyIn seems fast enough, but indirections are quite slow.
     for (int hi = 0; hi < m_Nhits; ++hi) {
       MatriplexHitPacker mhp(layerHits[hi].data());
 
       for (int i = beg; i < end; ++i) {
         const int hidx = tracks[i].getHitIdx(hi);
         const Hit& hit = layerHits[hi][hidx];
 
         m_HoTArr[hi](i - beg, 0, 0) = tracks[i].getHitOnTrack(hi);
 
         mhp.addInput(hit);
       }
 
       mhp.pack(m_msErr[hi], m_msPar[hi]);
     }
   }
 
   void MkFitter::inputTracksAndHitIdx(const std::vector<Track>& tracks, int beg, int end, bool inputProp) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     const int iI = inputProp ? iP : iC;
 
     int itrack = 0;
     for (int i = beg; i < end; ++i, ++itrack) {
       const Track& trk = tracks[i];
 
       m_Err[iI].copyIn(itrack, trk.errors().Array());
       m_Par[iI].copyIn(itrack, trk.parameters().Array());
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
       m_Label(itrack, 0, 0) = trk.label();
 
       for (int hi = 0; hi < m_Nhits; ++hi) {
         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
       }
     }
   }
 
   void MkFitter::inputTracksAndHitIdx(const std::vector<std::vector<Track> >& tracks,
                                       const std::vector<std::pair<int, int> >& idxs,
                                       int beg,
                                       int end,
                                       bool inputProp) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     const int iI = inputProp ? iP : iC;
 
     int itrack = 0;
     for (int i = beg; i < end; ++i, ++itrack) {
       const Track& trk = tracks[idxs[i].first][idxs[i].second];
 
       m_Label(itrack, 0, 0) = trk.label();
       m_SeedIdx(itrack, 0, 0) = idxs[i].first;
       m_CandIdx(itrack, 0, 0) = idxs[i].second;
 
       m_Err[iI].copyIn(itrack, trk.errors().Array());
       m_Par[iI].copyIn(itrack, trk.parameters().Array());
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
 
       for (int hi = 0; hi < m_Nhits; ++hi) {
         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
       }
     }
   }
 
   void MkFitter::inputSeedsTracksAndHits(const std::vector<Track>& seeds,
                                          const std::vector<Track>& tracks,
                                          const std::vector<HitVec>& layerHits,
                                          int beg,
                                          int end) {
     // Assign track parameters to initial state and copy hit values in.
 
     // This might not be true for the last chunk!
     // assert(end - beg == NN);
 
     int itrack;
     for (int i = beg; i < end; ++i) {
       itrack = i - beg;
 
       const Track& see = seeds[i];
 
       m_Label(itrack, 0, 0) = see.label();
       if (see.label() < 0)
         continue;
 
       m_Err[iC].copyIn(itrack, see.errors().Array());
       m_Par[iC].copyIn(itrack, see.parameters().Array());
 
       m_Chg(itrack, 0, 0) = see.charge();
       m_Chi2(itrack, 0, 0) = see.chi2();
 
       const Track& trk = tracks[see.label()];
 
       // CopyIn seems fast enough, but indirections are quite slow.
       for (int hi = 0; hi < m_Nhits; ++hi) {
         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
 
         const int hidx = trk.getHitIdx(hi);
         if (hidx < 0)
           continue;  //fixme, check if this is harmless
 
         const Hit& hit = layerHits[hi][hidx];
         m_msErr[hi].copyIn(itrack, hit.errArray());
         m_msPar[hi].copyIn(itrack, hit.posArray());
       }
     }
   }
 
   //------------------------------------------------------------------------------
   // Fitting with interleaved hit loading
   //------------------------------------------------------------------------------
 
   void MkFitter::inputTracksForFit(const std::vector<Track>& tracks, int beg, int end) {
     // Loads track parameters and hit indices.
 
     // XXXXMT4K has Config::nLayers: How many hits do we read in?
     // Check for max? Expect an argument?
     // What to do with invalid hits? Skip?
 
     // XXXX MT Here the same idx array WAS used for slurping in of tracks and
     // hots. With this, two index arrays are built, one within each packer.
 
     MatriplexTrackPacker mtp(&tracks[beg]);
     MatriplexHoTPacker mhotp(tracks[beg].getHitsOnTrackArray());
 
     int itrack = 0;
 
     for (int i = beg; i < end; ++i, ++itrack) {
       const Track& trk = tracks[i];
 
       m_Chg(itrack, 0, 0) = trk.charge();
       m_Chi2(itrack, 0, 0) = trk.chi2();
       m_Label(itrack, 0, 0) = trk.label();
 
       mtp.addInput(trk);
 
       mhotp.addInput(*trk.getHitsOnTrackArray());
     }
 
     mtp.pack(m_Err[iC], m_Par[iC]);
     for (int ll = 0; ll < Config::nLayers; ++ll) {
       mhotp.pack(m_HoTArr[ll], ll);
     }
   }
 
   void MkFitter::fitTracksWithInterSlurp(const std::vector<HitVec>& layersohits,
                                          const PropagationFlags& pflags,
                                          const int N_proc) {
     // XXXX This has potential issues hits coming from different layers!
     // Expected to only work reliably with barrel (consecutive layers from 0 -> m_Nhits)
     // and with hits present on every layer for every track.
 
     // Loops over layers and:
     // a) slurps in hit parameters;
     // b) propagates and updates tracks
 
     for (int ii = 0; ii < m_Nhits; ++ii) {
       // XXXX Assuming hit index corresponds to layer!
       MatriplexHitPacker mhp(layersohits[ii].data());
 
       for (int i = 0; i < N_proc; ++i) {
         const int hidx = m_HoTArr[ii](i, 0, 0).index;
         const int hlyr = m_HoTArr[ii](i, 0, 0).layer;
 
         // XXXXMT4K What to do with hidx < 0 ????
         // This could solve the unbalanced fit.
         // Or, if the hidx is the "universal" missing hit, it could just work.
         // Say, hidx = 0 ... grr ... but then we don't know it is missing.
 
         if (hidx < 0 || hlyr < 0) {
           mhp.addNullInput();
         } else {
           mhp.addInput(layersohits[hlyr][hidx]);
         }
       }
 
       mhp.pack(m_msErr[0], m_msPar[0]);
 
       propagateTracksToHitR(m_msPar[0], N_proc, pflags);
 
       kalmanUpdate(m_Err[iP], m_Par[iP], m_msErr[0], m_msPar[0], m_Err[iC], m_Par[iC], N_proc);
     }
   }
 
   //==============================================================================
   // Fitting functions
   //==============================================================================
 
   void MkFitter::outputTracks(std::vector<Track>& tracks, int beg, int end, int iCP) const {
     // Copies last track parameters (updated) into Track objects.
     // The tracks vector should be resized to allow direct copying.
 
     int itrack = 0;
     for (int i = beg; i < end; ++i, ++itrack) {
       m_Err[iCP].copyOut(itrack, tracks[i].errors_nc().Array());
       m_Par[iCP].copyOut(itrack, tracks[i].parameters_nc().Array());
 
       tracks[i].setCharge(m_Chg(itrack, 0, 0));
 
       // XXXXX chi2 is not set (also not in SMatrix fit, it seems)
       tracks[i].setChi2(m_Chi2(itrack, 0, 0));
       tracks[i].setLabel(m_Label(itrack, 0, 0));
     }
   }
 
   void MkFitter::outputFittedTracksAndHitIdx(std::vector<Track>& tracks, int beg, int end, bool outputProp) const {
     // Copies last track parameters (updated) into Track objects and up to m_Nhits.
     // The tracks vector should be resized to allow direct copying.
 
     const int iO = outputProp ? iP : iC;
 
     int itrack = 0;
     for (int i = beg; i < end; ++i, ++itrack) {
       m_Err[iO].copyOut(itrack, tracks[i].errors_nc().Array());
       m_Par[iO].copyOut(itrack, tracks[i].parameters_nc().Array());
 
       tracks[i].setCharge(m_Chg(itrack, 0, 0));
       tracks[i].setChi2(m_Chi2(itrack, 0, 0));
       tracks[i].setLabel(m_Label(itrack, 0, 0));
 
       // QQQQ Could do resize and std::copy, as in MkFinder::copy_out(), but
       // we do not know the correct N_found_hits.
       tracks[i].resetHits();
       tracks[i].reserveHits(m_Nhits);
       for (int hi = 0; hi < m_Nhits; ++hi) {
         tracks[i].addHitIdx(m_HoTArr[hi](itrack, 0, 0), 0.);
       }
     }
   }
 
 }  // end namespace mkfit

This page was automatically generated by the 2.2.1 LXR engine. The LXR team