Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-12 05:47:34

0001 #include "MkFitter.h"
0002 
0003 #include "KalmanUtilsMPlex.h"
0004 #include "MatriplexPackers.h"
0005 
0006 //#define DEBUG
0007 #include "Debug.h"
0008 
0009 #include <sstream>
0010 
0011 namespace mkfit {
0012 
0013   void MkFitter::checkAlignment() {
0014     printf("MkFitter alignment check:\n");
0015     Matriplex::align_check("  m_Err[0]   =", &m_Err[0].fArray[0]);
0016     Matriplex::align_check("  m_Err[1]   =", &m_Err[1].fArray[0]);
0017     Matriplex::align_check("  m_Par[0]   =", &m_Par[0].fArray[0]);
0018     Matriplex::align_check("  m_Par[1]   =", &m_Par[1].fArray[0]);
0019     Matriplex::align_check("  m_msErr[0] =", &m_msErr[0].fArray[0]);
0020     Matriplex::align_check("  m_msPar[0] =", &m_msPar[0].fArray[0]);
0021   }
0022 
0023   void MkFitter::printPt(int idx) {
0024     for (int i = 0; i < NN; ++i) {
0025       printf("%5.2f  ", std::hypot(m_Par[idx].At(i, 3, 0), m_Par[idx].At(i, 4, 0)));
0026     }
0027   }
0028 
0029   int MkFitter::countValidHits(int itrack, int end_hit) const {
0030     int result = 0;
0031     for (int hi = 0; hi < end_hit; ++hi) {
0032       if (m_HoTArr[hi](itrack, 0, 0).index >= 0)
0033         result++;
0034     }
0035     return result;
0036   }
0037 
0038   int MkFitter::countInvalidHits(int itrack, int end_hit) const {
0039     int result = 0;
0040     for (int hi = 0; hi < end_hit; ++hi) {
0041       // XXXX MT: Should also count -2 hits as invalid?
0042       if (m_HoTArr[hi](itrack, 0, 0).index == -1)
0043         result++;
0044     }
0045     return result;
0046   }
0047 
0048   //==============================================================================
0049 
0050   void MkFitter::inputTracksAndHits(const std::vector<Track>& tracks,
0051                                     const std::vector<HitVec>& layerHits,
0052                                     int beg,
0053                                     int end) {
0054     // Assign track parameters to initial state and copy hit values in.
0055 
0056     // This might not be true for the last chunk!
0057     // assert(end - beg == NN);
0058 
0059     int itrack = 0;
0060 
0061     for (int i = beg; i < end; ++i, ++itrack) {
0062       const Track& trk = tracks[i];
0063 
0064       m_Err[iC].copyIn(itrack, trk.errors().Array());
0065       m_Par[iC].copyIn(itrack, trk.parameters().Array());
0066 
0067       m_Chg(itrack, 0, 0) = trk.charge();
0068       m_Chi2(itrack, 0, 0) = trk.chi2();
0069       m_Label(itrack, 0, 0) = trk.label();
0070 
0071       // CopyIn seems fast enough, but indirections are quite slow.
0072       for (int hi = 0; hi < m_Nhits; ++hi) {
0073         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
0074 
0075         const int hidx = trk.getHitIdx(hi);
0076         if (hidx < 0)
0077           continue;
0078 
0079         const Hit& hit = layerHits[hi][hidx];
0080         m_msErr[hi].copyIn(itrack, hit.errArray());
0081         m_msPar[hi].copyIn(itrack, hit.posArray());
0082       }
0083     }
0084   }
0085 
0086   void MkFitter::inputTracksAndHits(const std::vector<Track>& tracks,
0087                                     const std::vector<LayerOfHits>& layerHits,
0088                                     int beg,
0089                                     int end) {
0090     // Assign track parameters to initial state and copy hit values in.
0091 
0092     // This might not be true for the last chunk!
0093     // assert(end - beg == NN);
0094 
0095     int itrack;
0096 
0097     for (int i = beg; i < end; ++i) {
0098       itrack = i - beg;
0099       const Track& trk = tracks[i];
0100 
0101       m_Label(itrack, 0, 0) = trk.label();
0102 
0103       m_Err[iC].copyIn(itrack, trk.errors().Array());
0104       m_Par[iC].copyIn(itrack, trk.parameters().Array());
0105 
0106       m_Chg(itrack, 0, 0) = trk.charge();
0107       m_Chi2(itrack, 0, 0) = trk.chi2();
0108 
0109       // CopyIn seems fast enough, but indirections are quite slow.
0110       for (int hi = 0; hi < m_Nhits; ++hi) {
0111         const int hidx = trk.getHitIdx(hi);
0112         const int hlyr = trk.getHitLyr(hi);
0113         const Hit& hit = layerHits[hlyr].refHit(hidx);
0114 
0115         m_msErr[hi].copyIn(itrack, hit.errArray());
0116         m_msPar[hi].copyIn(itrack, hit.posArray());
0117 
0118         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
0119       }
0120     }
0121   }
0122 
0123   void MkFitter::slurpInTracksAndHits(const std::vector<Track>& tracks,
0124                                       const std::vector<HitVec>& layerHits,
0125                                       int beg,
0126                                       int end) {
0127     // Assign track parameters to initial state and copy hit values in.
0128 
0129     // This might not be true for the last chunk!
0130     // assert(end - beg == NN);
0131 
0132     MatriplexTrackPacker mtp(&tracks[beg]);
0133 
0134     for (int i = beg; i < end; ++i) {
0135       int itrack = i - beg;
0136       const Track& trk = tracks[i];
0137 
0138       m_Label(itrack, 0, 0) = trk.label();
0139 
0140       mtp.addInput(trk);
0141 
0142       m_Chg(itrack, 0, 0) = trk.charge();
0143       m_Chi2(itrack, 0, 0) = trk.chi2();
0144     }
0145 
0146     mtp.pack(m_Err[iC], m_Par[iC]);
0147 
0148     // CopyIn seems fast enough, but indirections are quite slow.
0149     for (int hi = 0; hi < m_Nhits; ++hi) {
0150       MatriplexHitPacker mhp(layerHits[hi].data());
0151 
0152       for (int i = beg; i < end; ++i) {
0153         const int hidx = tracks[i].getHitIdx(hi);
0154         const Hit& hit = layerHits[hi][hidx];
0155 
0156         m_HoTArr[hi](i - beg, 0, 0) = tracks[i].getHitOnTrack(hi);
0157 
0158         mhp.addInput(hit);
0159       }
0160 
0161       mhp.pack(m_msErr[hi], m_msPar[hi]);
0162     }
0163   }
0164 
0165   void MkFitter::inputTracksAndHitIdx(const std::vector<Track>& tracks, int beg, int end, bool inputProp) {
0166     // Assign track parameters to initial state and copy hit values in.
0167 
0168     // This might not be true for the last chunk!
0169     // assert(end - beg == NN);
0170 
0171     const int iI = inputProp ? iP : iC;
0172 
0173     int itrack = 0;
0174     for (int i = beg; i < end; ++i, ++itrack) {
0175       const Track& trk = tracks[i];
0176 
0177       m_Err[iI].copyIn(itrack, trk.errors().Array());
0178       m_Par[iI].copyIn(itrack, trk.parameters().Array());
0179 
0180       m_Chg(itrack, 0, 0) = trk.charge();
0181       m_Chi2(itrack, 0, 0) = trk.chi2();
0182       m_Label(itrack, 0, 0) = trk.label();
0183 
0184       for (int hi = 0; hi < m_Nhits; ++hi) {
0185         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
0186       }
0187     }
0188   }
0189 
0190   void MkFitter::inputTracksAndHitIdx(const std::vector<std::vector<Track> >& tracks,
0191                                       const std::vector<std::pair<int, int> >& idxs,
0192                                       int beg,
0193                                       int end,
0194                                       bool inputProp) {
0195     // Assign track parameters to initial state and copy hit values in.
0196 
0197     // This might not be true for the last chunk!
0198     // assert(end - beg == NN);
0199 
0200     const int iI = inputProp ? iP : iC;
0201 
0202     int itrack = 0;
0203     for (int i = beg; i < end; ++i, ++itrack) {
0204       const Track& trk = tracks[idxs[i].first][idxs[i].second];
0205 
0206       m_Label(itrack, 0, 0) = trk.label();
0207       m_SeedIdx(itrack, 0, 0) = idxs[i].first;
0208       m_CandIdx(itrack, 0, 0) = idxs[i].second;
0209 
0210       m_Err[iI].copyIn(itrack, trk.errors().Array());
0211       m_Par[iI].copyIn(itrack, trk.parameters().Array());
0212 
0213       m_Chg(itrack, 0, 0) = trk.charge();
0214       m_Chi2(itrack, 0, 0) = trk.chi2();
0215 
0216       for (int hi = 0; hi < m_Nhits; ++hi) {
0217         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
0218       }
0219     }
0220   }
0221 
0222   void MkFitter::inputSeedsTracksAndHits(const std::vector<Track>& seeds,
0223                                          const std::vector<Track>& tracks,
0224                                          const std::vector<HitVec>& layerHits,
0225                                          int beg,
0226                                          int end) {
0227     // Assign track parameters to initial state and copy hit values in.
0228 
0229     // This might not be true for the last chunk!
0230     // assert(end - beg == NN);
0231 
0232     int itrack;
0233     for (int i = beg; i < end; ++i) {
0234       itrack = i - beg;
0235 
0236       const Track& see = seeds[i];
0237 
0238       m_Label(itrack, 0, 0) = see.label();
0239       if (see.label() < 0)
0240         continue;
0241 
0242       m_Err[iC].copyIn(itrack, see.errors().Array());
0243       m_Par[iC].copyIn(itrack, see.parameters().Array());
0244 
0245       m_Chg(itrack, 0, 0) = see.charge();
0246       m_Chi2(itrack, 0, 0) = see.chi2();
0247 
0248       const Track& trk = tracks[see.label()];
0249 
0250       // CopyIn seems fast enough, but indirections are quite slow.
0251       for (int hi = 0; hi < m_Nhits; ++hi) {
0252         m_HoTArr[hi](itrack, 0, 0) = trk.getHitOnTrack(hi);
0253 
0254         const int hidx = trk.getHitIdx(hi);
0255         if (hidx < 0)
0256           continue;  //fixme, check if this is harmless
0257 
0258         const Hit& hit = layerHits[hi][hidx];
0259         m_msErr[hi].copyIn(itrack, hit.errArray());
0260         m_msPar[hi].copyIn(itrack, hit.posArray());
0261       }
0262     }
0263   }
0264 
0265   //------------------------------------------------------------------------------
0266   // Fitting with interleaved hit loading
0267   //------------------------------------------------------------------------------
0268 
0269   void MkFitter::inputTracksForFit(const std::vector<Track>& tracks, int beg, int end) {
0270     // Loads track parameters and hit indices.
0271 
0272     // XXXXMT4K has Config::nLayers: How many hits do we read in?
0273     // Check for max? Expect an argument?
0274     // What to do with invalid hits? Skip?
0275 
0276     // XXXX MT Here the same idx array WAS used for slurping in of tracks and
0277     // hots. With this, two index arrays are built, one within each packer.
0278 
0279     MatriplexTrackPacker mtp(&tracks[beg]);
0280     MatriplexHoTPacker mhotp(tracks[beg].getHitsOnTrackArray());
0281 
0282     int itrack = 0;
0283 
0284     for (int i = beg; i < end; ++i, ++itrack) {
0285       const Track& trk = tracks[i];
0286 
0287       m_Chg(itrack, 0, 0) = trk.charge();
0288       m_Chi2(itrack, 0, 0) = trk.chi2();
0289       m_Label(itrack, 0, 0) = trk.label();
0290 
0291       mtp.addInput(trk);
0292 
0293       mhotp.addInput(*trk.getHitsOnTrackArray());
0294     }
0295 
0296     mtp.pack(m_Err[iC], m_Par[iC]);
0297     for (int ll = 0; ll < Config::nLayers; ++ll) {
0298       mhotp.pack(m_HoTArr[ll], ll);
0299     }
0300   }
0301 
0302   void MkFitter::fitTracksWithInterSlurp(const std::vector<HitVec>& layersohits, const int N_proc) {
0303     // XXXX This has potential issues hits coming from different layers!
0304     // Expected to only work reliably with barrel (consecutive layers from 0 -> m_Nhits)
0305     // and with hits present on every layer for every track.
0306 
0307     // Loops over layers and:
0308     // a) slurps in hit parameters;
0309     // b) propagates and updates tracks
0310 
0311     for (int ii = 0; ii < m_Nhits; ++ii) {
0312       // XXXX Assuming hit index corresponds to layer!
0313       MatriplexHitPacker mhp(layersohits[ii].data());
0314 
0315       for (int i = 0; i < N_proc; ++i) {
0316         const int hidx = m_HoTArr[ii](i, 0, 0).index;
0317         const int hlyr = m_HoTArr[ii](i, 0, 0).layer;
0318 
0319         // XXXXMT4K What to do with hidx < 0 ????
0320         // This could solve the unbalanced fit.
0321         // Or, if the hidx is the "universal" missing hit, it could just work.
0322         // Say, hidx = 0 ... grr ... but then we don't know it is missing.
0323 
0324         if (hidx < 0 || hlyr < 0) {
0325           mhp.addNullInput();
0326         } else {
0327           mhp.addInput(layersohits[hlyr][hidx]);
0328         }
0329       }
0330 
0331       mhp.pack(m_msErr[0], m_msPar[0]);
0332 
0333       propagateTracksToHitR(m_msPar[0], N_proc, PropagationConfig::get_default().forward_fit_pflags);
0334 
0335       kalmanUpdate(m_Err[iP], m_Par[iP], m_msErr[0], m_msPar[0], m_Err[iC], m_Par[iC], N_proc);
0336     }
0337   }
0338 
0339   //==============================================================================
0340   // Fitting functions
0341   //==============================================================================
0342 
0343   void MkFitter::outputTracks(std::vector<Track>& tracks, int beg, int end, int iCP) const {
0344     // Copies last track parameters (updated) into Track objects.
0345     // The tracks vector should be resized to allow direct copying.
0346 
0347     int itrack = 0;
0348     for (int i = beg; i < end; ++i, ++itrack) {
0349       m_Err[iCP].copyOut(itrack, tracks[i].errors_nc().Array());
0350       m_Par[iCP].copyOut(itrack, tracks[i].parameters_nc().Array());
0351 
0352       tracks[i].setCharge(m_Chg(itrack, 0, 0));
0353 
0354       // XXXXX chi2 is not set (also not in SMatrix fit, it seems)
0355       tracks[i].setChi2(m_Chi2(itrack, 0, 0));
0356       tracks[i].setLabel(m_Label(itrack, 0, 0));
0357     }
0358   }
0359 
0360   void MkFitter::outputFittedTracksAndHitIdx(std::vector<Track>& tracks, int beg, int end, bool outputProp) const {
0361     // Copies last track parameters (updated) into Track objects and up to m_Nhits.
0362     // The tracks vector should be resized to allow direct copying.
0363 
0364     const int iO = outputProp ? iP : iC;
0365 
0366     int itrack = 0;
0367     for (int i = beg; i < end; ++i, ++itrack) {
0368       m_Err[iO].copyOut(itrack, tracks[i].errors_nc().Array());
0369       m_Par[iO].copyOut(itrack, tracks[i].parameters_nc().Array());
0370 
0371       tracks[i].setCharge(m_Chg(itrack, 0, 0));
0372       tracks[i].setChi2(m_Chi2(itrack, 0, 0));
0373       tracks[i].setLabel(m_Label(itrack, 0, 0));
0374 
0375       // QQQQ Could do resize and std::copy, as in MkFinder::copy_out(), but
0376       // we do not know the correct N_found_hits.
0377       tracks[i].resetHits();
0378       tracks[i].reserveHits(m_Nhits);
0379       for (int hi = 0; hi < m_Nhits; ++hi) {
0380         tracks[i].addHitIdx(m_HoTArr[hi](itrack, 0, 0), 0.);
0381       }
0382     }
0383   }
0384 
0385 }  // end namespace mkfit