MkFitCore/src/MkFinder.h

0001 #ifndef RecoTracker_MkFitCore_src_MkFinder_h
0002 #define RecoTracker_MkFitCore_src_MkFinder_h
0003
0004 #include "MkBase.h"
0005 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0006 #include "RecoTracker/MkFitCore/interface/Track.h"
0007
0008 #include "RecoTracker/MkFitCore/interface/HitStructures.h"
0009 #include "RecoTracker/MkFitCore/interface/TrackStructures.h"
0010
0011 // Define to get printouts about track and hit chi2.
0012 // See also MkBuilder::backwardFit().
0013
0014 //#define DEBUG_BACKWARD_FIT_BH
0015 //#define DEBUG_BACKWARD_FIT
0016
0017 namespace mkfit {
0018
0019   class CandCloner;
0020   class CombCandidate;
0021   class LayerOfHits;
0022   class FindingFoos;
0023   class IterationParams;
0024   class IterationLayerConfig;
0025   class SteeringParams;
0026
0027   class Event;
0028
0029   struct UpdateIndices {
0030     int seed_idx;
0031     int cand_idx;
0032     int hit_idx;
0033     int ovlp_idx;
0034
0035     UpdateIndices(int si, int ci, int hi, int oi) : seed_idx(si), cand_idx(ci), hit_idx(hi), ovlp_idx(oi) {}
0036   };
0037
0038   class MkFinder : public MkBase {
0039     friend class MkBuilder;
0040
0041   public:
0042     static constexpr int MPlexHitIdxMax = 16;
0043
0044     using MPlexHitIdx = Matriplex::Matriplex<int, MPlexHitIdxMax, 1, NN>;
0045     using MPlexQHoT = Matriplex::Matriplex<HitOnTrack, 1, 1, NN>;
0046
0047     //----------------------------------------------------------------------------
0048
0049     MkFinder() {}
0050
0051     void setup(const PropagationConfig &pc,
0052                const IterationConfig &ic,
0053                const IterationParams &ip,
0054                const IterationLayerConfig &ilc,
0055                const SteeringParams &sp,
0056                const std::vector<bool> *ihm,
0057                const Event *ev,
0058                int region,
0059                bool infwd);
0060     void setup_bkfit(const PropagationConfig &pc, const SteeringParams &sp, const Event *ev);
0061     void release();
0062
0063     void begin_layer(const LayerOfHits &layer_of_hits);
0064     void end_layer();
0065
0066     //----------------------------------------------------------------------------
0067
0068     void inputTracksAndHitIdx(const std::vector<Track> &tracks, int beg, int end, bool inputProp);
0069
0070     void inputTracksAndHitIdx(const std::vector<Track> &tracks,
0071                               const std::vector<int> &idxs,
0072                               int beg,
0073                               int end,
0074                               bool inputProp,
0075                               int mp_offset);
0076
0077     void inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0078                               const std::vector<std::pair<int, int>> &idxs,
0079                               int beg,
0080                               int end,
0081                               bool inputProp);
0082
0083     void inputTracksAndHits(const std::vector<CombCandidate> &tracks,
0084                             const LayerOfHits &layer_of_hits,
0085                             const std::vector<UpdateIndices> &idxs,
0086                             int beg,
0087                             int end,
0088                             bool inputProp);
0089     void inputOverlapHits(const LayerOfHits &layer_of_hits, const std::vector<UpdateIndices> &idxs, int beg, int end);
0090
0091     void inputTracksAndHitIdx(const std::vector<CombCandidate> &tracks,
0092                               const std::vector<std::pair<int, IdxChi2List>> &idxs,
0093                               int beg,
0094                               int end,
0095                               bool inputProp);
0096
0097     void outputTracksAndHitIdx(std::vector<Track> &tracks, int beg, int end, bool outputProp) const;
0098
0099     void outputTracksAndHitIdx(
0100         std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool outputProp) const;
0101
0102     void outputTrackAndHitIdx(Track &track, int itrack, bool outputProp) const {
0103       const int iO = outputProp ? iP : iC;
0104       copy_out(track, itrack, iO);
0105     }
0106
0107     void outputNonStoppedTracksAndHitIdx(
0108         std::vector<Track> &tracks, const std::vector<int> &idxs, int beg, int end, bool outputProp) const {
0109       const int iO = outputProp ? iP : iC;
0110       for (int i = beg, imp = 0; i < end; ++i, ++imp) {
0111         if (!m_Stopped[imp])
0112           copy_out(tracks[idxs[i]], imp, iO);
0113       }
0114     }
0115
0116     HitOnTrack bestHitLastHoT(int itrack) const { return m_HoTArrs[itrack][m_NHits(itrack, 0, 0) - 1]; }
0117
0118     void packModuleNormDirPnt(
0119         const LayerOfHits &layer_of_hits, int hit_cnt, MPlexHV &norm, MPlexHV &dir, MPlexHV &pnt, int N_proc) const;
0120
0121     //----------------------------------------------------------------------------
0122
0123     void getHitSelDynamicWindows(
0124         const float invpt, const float theta, float &min_dq, float &max_dq, float &min_dphi, float &max_dphi);
0125
0126     float getHitSelDynamicChi2Cut(const int itrk, const int ipar);
0127
0128     void selectHitIndices(const LayerOfHits &layer_of_hits, const int N_proc, bool fill_binsearch_only = false);
0129     void selectHitIndicesV2(const LayerOfHits &layer_of_hits, const int N_proc);
0130
0131     void addBestHit(const LayerOfHits &layer_of_hits, const int N_proc, const FindingFoos &fnd_foos);
0132
0133     //----------------------------------------------------------------------------
0134
0135     void findCandidates(const LayerOfHits &layer_of_hits,
0136                         std::vector<std::vector<TrackCand>> &tmp_candidates,
0137                         const int offset,
0138                         const int N_proc,
0139                         const FindingFoos &fnd_foos);
0140
0141     //----------------------------------------------------------------------------
0142
0143     void findCandidatesCloneEngine(const LayerOfHits &layer_of_hits,
0144                                    CandCloner &cloner,
0145                                    const int offset,
0146                                    const int N_proc,
0147                                    const FindingFoos &fnd_foos);
0148
0149     void updateWithLoadedHit(int N_proc, const LayerOfHits &layer_of_hits, const FindingFoos &fnd_foos);
0150
0151     void chi2OfLoadedHit(int N_proc, const FindingFoos &fnd_foos);
0152
0153     void copyOutParErr(std::vector<CombCandidate> &seed_cand_vec, int N_proc, bool outputProp) const;
0154
0155     //----------------------------------------------------------------------------
0156     // Backward fit
0157
0158     void bkFitInputTracks(TrackVec &cands, int beg, int end);
0159     void bkFitOutputTracks(TrackVec &cands, int beg, int end, bool outputProp);
0160
0161     void bkFitInputTracks(EventOfCombCandidates &eocss, int beg, int end);
0162     void bkFitOutputTracks(EventOfCombCandidates &eocss, int beg, int end, bool outputProp);
0163
0164     void bkFitFitTracksBH(const EventOfHits &eventofhits,
0165                           const SteeringParams &st_par,
0166                           const int N_proc,
0167                           bool chiDebug = false);
0168
0169     void bkFitFitTracks(const EventOfHits &eventofhits,
0170                         const SteeringParams &st_par,
0171                         const int N_proc,
0172                         bool chiDebug = false);
0173
0174     void bkFitPropTracksToPCA(const int N_proc);
0175
0176     //----------------------------------------------------------------------------
0177
0178   private:
0179     void copy_in(const Track &trk, const int mslot, const int tslot) {
0180       m_Err[tslot].copyIn(mslot, trk.errors().Array());
0181       m_Par[tslot].copyIn(mslot, trk.parameters().Array());
0182
0183       m_Chg(mslot, 0, 0) = trk.charge();
0184       m_Chi2(mslot, 0, 0) = trk.chi2();
0185       m_Label(mslot, 0, 0) = trk.label();
0186
0187       m_NHits(mslot, 0, 0) = trk.nTotalHits();
0188       m_NFoundHits(mslot, 0, 0) = trk.nFoundHits();
0189
0190       m_NInsideMinusOneHits(mslot, 0, 0) = trk.nInsideMinusOneHits();
0191       m_NTailMinusOneHits(mslot, 0, 0) = trk.nTailMinusOneHits();
0192
0193       std::copy(trk.beginHitsOnTrack(), trk.endHitsOnTrack(), m_HoTArrs[mslot]);
0194     }
0195
0196     void copy_out(Track &trk, const int mslot, const int tslot) const {
0197       m_Err[tslot].copyOut(mslot, trk.errors_nc().Array());
0198       m_Par[tslot].copyOut(mslot, trk.parameters_nc().Array());
0199
0200       trk.setCharge(m_Chg(mslot, 0, 0));
0201       trk.setChi2(m_Chi2(mslot, 0, 0));
0202       trk.setLabel(m_Label(mslot, 0, 0));
0203
0204       trk.resizeHits(m_NHits(mslot, 0, 0), m_NFoundHits(mslot, 0, 0));
0205       std::copy(m_HoTArrs[mslot], &m_HoTArrs[mslot][m_NHits(mslot, 0, 0)], trk.beginHitsOnTrack_nc());
0206     }
0207
0208     void copy_in(const TrackCand &trk, const int mslot, const int tslot) {
0209       m_Err[tslot].copyIn(mslot, trk.errors().Array());
0210       m_Par[tslot].copyIn(mslot, trk.parameters().Array());
0211
0212       m_Chg(mslot, 0, 0) = trk.charge();
0213       m_Chi2(mslot, 0, 0) = trk.chi2();
0214       m_Label(mslot, 0, 0) = trk.label();
0215
0216       m_LastHitCcIndex(mslot, 0, 0) = trk.lastCcIndex();
0217       m_NFoundHits(mslot, 0, 0) = trk.nFoundHits();
0218       m_NMissingHits(mslot, 0, 0) = trk.nMissingHits();
0219       m_NOverlapHits(mslot, 0, 0) = trk.nOverlapHits();
0220
0221       m_NInsideMinusOneHits(mslot, 0, 0) = trk.nInsideMinusOneHits();
0222       m_NTailMinusOneHits(mslot, 0, 0) = trk.nTailMinusOneHits();
0223
0224       m_CombCand[mslot] = trk.combCandidate();
0225       m_TrkStatus[mslot] = trk.getStatus();
0226     }
0227
0228     void copy_out(TrackCand &trk, const int mslot, const int tslot) const {
0229       m_Err[tslot].copyOut(mslot, trk.errors_nc().Array());
0230       m_Par[tslot].copyOut(mslot, trk.parameters_nc().Array());
0231
0232       trk.setCharge(m_Chg(mslot, 0, 0));
0233       trk.setChi2(m_Chi2(mslot, 0, 0));
0234       trk.setLabel(m_Label(mslot, 0, 0));
0235
0236       trk.setLastCcIndex(m_LastHitCcIndex(mslot, 0, 0));
0237       trk.setNFoundHits(m_NFoundHits(mslot, 0, 0));
0238       trk.setNMissingHits(m_NMissingHits(mslot, 0, 0));
0239       trk.setNOverlapHits(m_NOverlapHits(mslot, 0, 0));
0240
0241       trk.setNInsideMinusOneHits(m_NInsideMinusOneHits(mslot, 0, 0));
0242       trk.setNTailMinusOneHits(m_NTailMinusOneHits(mslot, 0, 0));
0243
0244       trk.setCombCandidate(m_CombCand[mslot]);
0245       trk.setStatus(m_TrkStatus[mslot]);
0246     }
0247
0248     void add_hit(const int mslot, int index, int layer) {
0249       // Only used by BestHit.
0250       // m_NInsideMinusOneHits and m_NTailMinusOneHits are maintained here but are
0251       // not used and are not copied out (as Track does not have these members).
0252
0253       int &n_tot_hits = m_NHits(mslot, 0, 0);
0254       int &n_fnd_hits = m_NFoundHits(mslot, 0, 0);
0255
0256       if (n_tot_hits < Config::nMaxTrkHits) {
0257         m_HoTArrs[mslot][n_tot_hits++] = {index, layer};
0258         if (index >= 0) {
0259           ++n_fnd_hits;
0260           m_NInsideMinusOneHits(mslot, 0, 0) += m_NTailMinusOneHits(mslot, 0, 0);
0261           m_NTailMinusOneHits(mslot, 0, 0) = 0;
0262         } else if (index == -1) {
0263           ++m_NTailMinusOneHits(mslot, 0, 0);
0264         }
0265       } else {
0266         // printf("WARNING MkFinder::add_hit hit-on-track limit reached for label=%d\n", label_);
0267
0268         const int LH = Config::nMaxTrkHits - 1;
0269
0270         if (index >= 0) {
0271           if (m_HoTArrs[mslot][LH].index < 0)
0272             ++n_fnd_hits;
0273           m_HoTArrs[mslot][LH] = {index, layer};
0274         } else if (index == -2) {
0275           if (m_HoTArrs[mslot][LH].index >= 0)
0276             --n_fnd_hits;
0277           m_HoTArrs[mslot][LH] = {index, layer};
0278         }
0279       }
0280     }
0281
0282     int num_all_minus_one_hits(const int mslot) const {
0283       return m_NInsideMinusOneHits(mslot, 0, 0) + m_NTailMinusOneHits(mslot, 0, 0);
0284     }
0285
0286     int num_inside_minus_one_hits(const int mslot) const { return m_NInsideMinusOneHits(mslot, 0, 0); }
0287
0288     void print_par_err(int corp, int mslot) const;
0289
0290     //----------------------------------------------------------------------------
0291
0292     MPlexQF m_Chi2;
0293     MPlexQI m_Label;  // seed index in global seed vector (for MC truth match)
0294
0295     MPlexQI m_NHits;
0296     MPlexQI m_NFoundHits;
0297
0298     HitOnTrack m_HoTArrs[NN][Config::nMaxTrkHits];
0299
0300     MPlexQI m_SeedIdx;        // seed index in local thread (for bookkeeping at thread level)
0301     MPlexQI m_CandIdx;        // candidate index for the given seed (for bookkeeping of clone engine)
0302     MPlexQI m_SeedOriginIdx;  // seed index in MkBuilder seed input vector
0303
0304     MPlexQI m_Stopped;  // Flag for BestHit that a track has been stopped (and copied out already)
0305
0306     // Additions / substitutions for TrackCand copy_in/out()
0307     // One could really access the original TrackCand for all of those, especially the ones that
0308     // are STD only. This then requires access back to that TrackCand memory.
0309     // So maybe one should just have flags for CopyIn methods (or several versions). Yay, etc.
0310     MPlexQI m_NMissingHits;             // sub: m_NHits, sort of, STD only
0311     MPlexQI m_NOverlapHits;             // add: num of overlaps registered in HitOnTrack, STD only
0312     MPlexQI m_NInsideMinusOneHits;      // sub: before we copied all hit idcs and had a loop counting them only
0313     MPlexQI m_NTailMinusOneHits;        // sub: before we copied all hit idcs and had a loop counting them only
0314     MPlexQI m_LastHitCcIndex;           // add: index of last hit in m_CombCand hit tree, STD only
0315     TrackBase::Status m_TrkStatus[NN];  // STD only, status bits
0316     CombCandidate *m_CombCand[NN];
0317     // const TrackCand *m_TrkCand[NN]; // hmmh, could get all data through this guy ... but scattered
0318     // storing it in now for bkfit debug printouts
0319     TrackCand *m_TrkCand[NN];
0320
0321     // Hit indices into LayerOfHits to explore.
0322     WSR_Result m_XWsrResult[NN];  // Could also merge it with m_XHitSize. Or use smaller arrays.
0323     MPlexQI m_XHitSize;
0324     MPlexHitIdx m_XHitArr;
0325
0326     // Hit errors / parameters for hit matching, update.
0327     MPlexHS m_msErr{0.0f};
0328     MPlexHV m_msPar{0.0f};
0329
0330     // An idea: Do propagation to hit in FindTracksXYZZ functions.
0331     // Have some state / functions here that make this short to write.
0332     // This would simplify KalmanUtils (remove the propagate functions).
0333     // Track errors / parameters propagated to current hit.
0334     // MPlexLS    candErrAtCurrHit;
0335     // MPlexLV    candParAtCurrHit;
0336
0337     const PropagationConfig *m_prop_config = nullptr;
0338     const IterationConfig *m_iteration_config = nullptr;
0339     const IterationParams *m_iteration_params = nullptr;
0340     const IterationLayerConfig *m_iteration_layer_config = nullptr;
0341     const SteeringParams *m_steering_params = nullptr;
0342     const std::vector<bool> *m_iteration_hit_mask = nullptr;
0343     const Event *m_event = nullptr;
0344     int m_current_region = -1;
0345     bool m_in_fwd = true;
0346
0347     // Backward fit
0348     int m_CurHit[NN];
0349     const HitOnTrack *m_HoTArr[NN];
0350     int m_CurNode[NN];
0351     const HoTNode *m_HoTNodeArr[NN];
0352   };
0353
0354 }  // end namespace mkfit
0355 #endif