Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-25 02:14:12 

0001 #ifndef RecoTracker_MkFitCore_interface_Track_h
0002 #define RecoTracker_MkFitCore_interface_Track_h
0003 
0004 #include "RecoTracker/MkFitCore/interface/Config.h"
0005 #include "RecoTracker/MkFitCore/interface/MatrixSTypes.h"
0006 #include "RecoTracker/MkFitCore/interface/FunctionTypes.h"
0007 #include "RecoTracker/MkFitCore/interface/Hit.h"
0008 #include "RecoTracker/MkFitCore/interface/IdxChi2List.h"
0009 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0010 
0011 #include <vector>
0012 #include <map>
0013 #include <limits>
0014 
0015 namespace mkfit {
0016 
0017   typedef std::pair<int, int> SimTkIDInfo;
0018   typedef std::vector<int> HitIdxVec;
0019   typedef std::map<int, std::vector<int> > HitLayerMap;
0020 
0021   inline int calculateCharge(const Hit& hit0, const Hit& hit1, const Hit& hit2) {
0022     return ((hit2.y() - hit0.y()) * (hit2.x() - hit1.x()) > (hit2.y() - hit1.y()) * (hit2.x() - hit0.x()) ? 1 : -1);
0023   }
0024 
0025   inline int calculateCharge(const float hit0_x,
0026                              const float hit0_y,
0027                              const float hit1_x,
0028                              const float hit1_y,
0029                              const float hit2_x,
0030                              const float hit2_y) {
0031     return ((hit2_y - hit0_y) * (hit2_x - hit1_x) > (hit2_y - hit1_y) * (hit2_x - hit0_x) ? 1 : -1);
0032   }
0033 
0034   //==============================================================================
0035   // TrackState
0036   //==============================================================================
0037 
0038   struct TrackState  //  possible to add same accessors as track?
0039   {
0040   public:
0041     TrackState() : valid(true) {}
0042     TrackState(int charge, const SVector3& pos, const SVector3& mom, const SMatrixSym66& err)
0043         : parameters(SVector6(pos.At(0), pos.At(1), pos.At(2), mom.At(0), mom.At(1), mom.At(2))),
0044           errors(err),
0045           charge(charge),
0046           valid(true) {}
0047     SVector3 position() const { return SVector3(parameters[0], parameters[1], parameters[2]); }
0048     SVector6 parameters;
0049     SMatrixSym66 errors;
0050     short charge;
0051     bool valid;
0052 
0053     // track state position
0054     float x() const { return parameters.At(0); }
0055     float y() const { return parameters.At(1); }
0056     float z() const { return parameters.At(2); }
0057     float posR() const { return getHypot(x(), y()); }
0058     float posRsq() const { return x() * x() + y() * y(); }
0059     float posPhi() const { return getPhi(x(), y()); }
0060     float posEta() const { return getEta(posR(), z()); }
0061 
0062     // track state position errors
0063     float exx() const { return std::sqrt(errors.At(0, 0)); }
0064     float eyy() const { return std::sqrt(errors.At(1, 1)); }
0065     float ezz() const { return std::sqrt(errors.At(2, 2)); }
0066     float exy() const { return std::sqrt(errors.At(0, 1)); }
0067     float exz() const { return std::sqrt(errors.At(0, 2)); }
0068     float eyz() const { return std::sqrt(errors.At(1, 2)); }
0069 
0070     float eposR() const { return std::sqrt(getRadErr2(x(), y(), errors.At(0, 0), errors.At(1, 1), errors.At(0, 1))); }
0071     float eposPhi() const { return std::sqrt(getPhiErr2(x(), y(), errors.At(0, 0), errors.At(1, 1), errors.At(0, 1))); }
0072     float eposEta() const {
0073       return std::sqrt(getEtaErr2(x(),
0074                                   y(),
0075                                   z(),
0076                                   errors.At(0, 0),
0077                                   errors.At(1, 1),
0078                                   errors.At(2, 2),
0079                                   errors.At(0, 1),
0080                                   errors.At(0, 2),
0081                                   errors.At(1, 2)));
0082     }
0083 
0084     // track state momentum
0085     float invpT() const { return parameters.At(3); }
0086     float momPhi() const { return parameters.At(4); }
0087     float theta() const { return parameters.At(5); }
0088     float pT() const { return std::abs(1.f / parameters.At(3)); }
0089     float px() const { return pT() * std::cos(parameters.At(4)); }
0090     float py() const { return pT() * std::sin(parameters.At(4)); }
0091     float pz() const { return pT() / std::tan(parameters.At(5)); }
0092     float momEta() const { return getEta(theta()); }
0093     float p() const { return pT() / std::sin(parameters.At(5)); }
0094 
0095     float einvpT() const { return std::sqrt(errors.At(3, 3)); }
0096     float emomPhi() const { return std::sqrt(errors.At(4, 4)); }
0097     float etheta() const { return std::sqrt(errors.At(5, 5)); }
0098     float epT() const { return std::sqrt(errors.At(3, 3)) / (parameters.At(3) * parameters.At(3)); }
0099     float emomEta() const { return std::sqrt(errors.At(5, 5)) / std::sin(parameters.At(5)); }
0100     float epxpx() const { return std::sqrt(getPxPxErr2(invpT(), momPhi(), errors.At(3, 3), errors.At(4, 4))); }
0101     float epypy() const { return std::sqrt(getPyPyErr2(invpT(), momPhi(), errors.At(3, 3), errors.At(4, 4))); }
0102     float epzpz() const { return std::sqrt(getPyPyErr2(invpT(), theta(), errors.At(3, 3), errors.At(5, 5))); }
0103 
0104     void convertFromCartesianToCCS();
0105     void convertFromCCSToCartesian();
0106     SMatrix66 jacobianCCSToCartesian(float invpt, float phi, float theta) const;
0107     SMatrix66 jacobianCartesianToCCS(float px, float py, float pz) const;
0108 
0109     void convertFromGlbCurvilinearToCCS();
0110     void convertFromCCSToGlbCurvilinear();
0111     //last row/column are zeros
0112     SMatrix66 jacobianCCSToCurvilinear(float invpt, float cosP, float sinP, float cosT, float sinT, short charge) const;
0113     SMatrix66 jacobianCurvilinearToCCS(float px, float py, float pz, short charge) const;
0114   };
0115 
0116   //==============================================================================
0117   // TrackBase
0118   //==============================================================================
0119 
0120   class TrackBase {
0121   public:
0122     TrackBase() {}
0123 
0124     TrackBase(const TrackState& state, float chi2, int label) : state_(state), chi2_(chi2), label_(label) {}
0125 
0126     TrackBase(int charge, const SVector3& position, const SVector3& momentum, const SMatrixSym66& errors, float chi2)
0127         : state_(charge, position, momentum, errors), chi2_(chi2) {}
0128 
0129     const TrackState& state() const { return state_; }
0130     void setState(const TrackState& newState) { state_ = newState; }
0131 
0132     const SVector6& parameters() const { return state_.parameters; }
0133     const SMatrixSym66& errors() const { return state_.errors; }
0134 
0135     const float* posArray() const { return state_.parameters.Array(); }
0136     const float* errArray() const { return state_.errors.Array(); }
0137 
0138     // Non-const versions needed for CopyOut of Matriplex.
0139     SVector6& parameters_nc() { return state_.parameters; }
0140     SMatrixSym66& errors_nc() { return state_.errors; }
0141     TrackState& state_nc() { return state_; }
0142 
0143     SVector3 position() const { return SVector3(state_.parameters[0], state_.parameters[1], state_.parameters[2]); }
0144     SVector3 momentum() const { return SVector3(state_.parameters[3], state_.parameters[4], state_.parameters[5]); }
0145 
0146     float x() const { return state_.parameters[0]; }
0147     float y() const { return state_.parameters[1]; }
0148     float z() const { return state_.parameters[2]; }
0149     float posR() const { return getHypot(state_.parameters[0], state_.parameters[1]); }
0150     float posRsq() const { return state_.posRsq(); }
0151     float posPhi() const { return getPhi(state_.parameters[0], state_.parameters[1]); }
0152     float posEta() const { return getEta(state_.parameters[0], state_.parameters[1], state_.parameters[2]); }
0153 
0154     float px() const { return state_.px(); }
0155     float py() const { return state_.py(); }
0156     float pz() const { return state_.pz(); }
0157     float pT() const { return state_.pT(); }
0158     float invpT() const { return state_.invpT(); }
0159     float p() const { return state_.p(); }
0160     float momPhi() const { return state_.momPhi(); }
0161     float momEta() const { return state_.momEta(); }
0162     float theta() const { return state_.theta(); }
0163 
0164     // track state momentum errors
0165     float epT() const { return state_.epT(); }
0166     float emomPhi() const { return state_.emomPhi(); }
0167     float emomEta() const { return state_.emomEta(); }
0168 
0169     // ------------------------------------------------------------------------
0170 
0171     int charge() const { return state_.charge; }
0172     float chi2() const { return chi2_; }
0173     float score() const { return score_; }
0174     int label() const { return label_; }
0175 
0176     void setCharge(int chg) { state_.charge = chg; }
0177     void setChi2(float chi2) { chi2_ = chi2; }
0178     void setScore(float s) { score_ = s; }
0179     void setLabel(int lbl) { label_ = lbl; }
0180 
0181     bool hasSillyValues(bool dump, bool fix, const char* pref = "");
0182     bool hasNanNSillyValues() const;
0183 
0184     // ------------------------------------------------------------------------
0185 
0186     float d0BeamSpot(const float x_bs, const float y_bs, bool linearize = false) const;
0187 
0188     // used for swimming cmssw rec tracks to mkFit position
0189     float swimPhiToR(const float x, const float y) const;
0190 
0191     bool canReachRadius(float R) const;
0192     float maxReachRadius() const;
0193     float zAtR(float R, float* r_reached = nullptr) const;
0194     float rAtZ(float Z) const;
0195 
0196     // ------------------------------------------------------------------------
0197 
0198     struct Status {
0199       static constexpr int kNSeedHitBits = 4;
0200       static constexpr int kMaxSeedHits = (1 << kNSeedHitBits) - 1;
0201 
0202       // Set to true for short, low-pt CMS tracks. They do not generate mc seeds and
0203       // do not enter the efficiency denominator.
0204       bool not_findable : 1;
0205 
0206       // Set to true when number of holes would exceed an external limit, Config::maxHolesPerCand.
0207       // XXXXMT Not used yet, -2 last hit idx is still used! Need to add it to MkFi**r classes.
0208       // Problem is that I have to carry bits in/out of the MkFinder, too.
0209       bool stopped : 1;
0210 
0211       // Production type (most useful for sim tracks): 0, 1, 2, 3 for unset, signal, in-time PU, oot PU
0212       unsigned int prod_type : 2;
0213 
0214       unsigned int align_was_seed_type : 2;
0215 
0216       // Whether or not the track matched to another track and had the lower cand score
0217       bool duplicate : 1;
0218 
0219       // Tracking iteration/algorithm
0220       unsigned int algorithm : 6;
0221 
0222       // Temporary store number of overlaps for Track here
0223       int n_overlaps : 8;
0224 
0225       // Number of seed hits at import time
0226       unsigned int n_seed_hits : kNSeedHitBits;
0227 
0228       // mkFit tracking region TrackerInfo::EtaRegion, determined by seed partition function
0229       unsigned int eta_region : 3;
0230 
0231       // The remaining bits.
0232       unsigned int _free_bits_ : 4;
0233 
0234       Status()
0235           : not_findable(false),
0236             stopped(false),
0237             prod_type(0),
0238             align_was_seed_type(0),
0239             duplicate(false),
0240             algorithm(0),
0241             n_overlaps(0),
0242             n_seed_hits(0),
0243             eta_region(0),
0244             _free_bits_(0) {}
0245     };
0246     static_assert(sizeof(Status) == sizeof(int));
0247 
0248     Status getStatus() const { return status_; }
0249     void setStatus(Status s) { status_ = s; }
0250 
0251     bool isFindable() const { return !status_.not_findable; }
0252     bool isNotFindable() const { return status_.not_findable; }
0253     void setNotFindable() { status_.not_findable = true; }
0254 
0255     void setDuplicateValue(bool d) { status_.duplicate = d; }
0256     bool getDuplicateValue() const { return status_.duplicate; }
0257     enum class ProdType { NotSet = 0, Signal = 1, InTimePU = 2, OutOfTimePU = 3 };
0258     ProdType prodType() const { return ProdType(status_.prod_type); }
0259     void setProdType(ProdType ptyp) { status_.prod_type = static_cast<unsigned int>(ptyp); }
0260 
0261     int getNSeedHits() const { return status_.n_seed_hits; }
0262     void setNSeedHits(int n) { status_.n_seed_hits = n; }
0263     int getEtaRegion() const { return status_.eta_region; }
0264     void setEtaRegion(int r) { status_.eta_region = r; }
0265 
0266     // Those are defined in Track, TrackCand has separate member. To be consolidated but
0267     // it's a binary format change.
0268     // int  nOverlapHits()  const  { return status_.n_overlaps; }
0269     // void setNOverlapHits(int n) { status_.n_overlaps = n; }
0270 
0271     /// track algorithm; copy from TrackBase.h to keep in standalone builds
0272     enum class TrackAlgorithm {
0273       undefAlgorithm = 0,
0274       ctf = 1,
0275       duplicateMerge = 2,
0276       cosmics = 3,
0277       initialStep = 4,
0278       lowPtTripletStep = 5,
0279       pixelPairStep = 6,
0280       detachedTripletStep = 7,
0281       mixedTripletStep = 8,
0282       pixelLessStep = 9,
0283       tobTecStep = 10,
0284       jetCoreRegionalStep = 11,
0285       conversionStep = 12,
0286       muonSeededStepInOut = 13,
0287       muonSeededStepOutIn = 14,
0288       outInEcalSeededConv = 15,
0289       inOutEcalSeededConv = 16,
0290       nuclInter = 17,
0291       standAloneMuon = 18,
0292       globalMuon = 19,
0293       cosmicStandAloneMuon = 20,
0294       cosmicGlobalMuon = 21,
0295       // Phase1
0296       highPtTripletStep = 22,
0297       lowPtQuadStep = 23,
0298       detachedQuadStep = 24,
0299       reservedForUpgrades1 = 25,
0300       reservedForUpgrades2 = 26,
0301       bTagGhostTracks = 27,
0302       beamhalo = 28,
0303       gsf = 29,
0304       // HLT algo name
0305       hltPixel = 30,
0306       // steps used by PF
0307       hltIter0 = 31,
0308       hltIter1 = 32,
0309       hltIter2 = 33,
0310       hltIter3 = 34,
0311       hltIter4 = 35,
0312       // steps used by all other objects @HLT
0313       hltIterX = 36,
0314       // steps used by HI muon regional iterative tracking
0315       hiRegitMuInitialStep = 37,
0316       hiRegitMuLowPtTripletStep = 38,
0317       hiRegitMuPixelPairStep = 39,
0318       hiRegitMuDetachedTripletStep = 40,
0319       hiRegitMuMixedTripletStep = 41,
0320       hiRegitMuPixelLessStep = 42,
0321       hiRegitMuTobTecStep = 43,
0322       hiRegitMuMuonSeededStepInOut = 44,
0323       hiRegitMuMuonSeededStepOutIn = 45,
0324       algoSize = 46
0325     };
0326 
0327     int algoint() const { return status_.algorithm; }
0328     TrackAlgorithm algorithm() const { return TrackAlgorithm(status_.algorithm); }
0329     void setAlgorithm(TrackAlgorithm algo) { status_.algorithm = static_cast<unsigned int>(algo); }
0330     void setAlgoint(int algo) { status_.algorithm = algo; }
0331     // To be used later
0332     // bool isStopped() const { return status_.stopped; }
0333     // void setStopped()      { status_.stopped = true; }
0334 
0335     static const char* algoint_to_cstr(int algo);
0336 
0337     // ------------------------------------------------------------------------
0338 
0339   protected:
0340     TrackState state_;
0341     float chi2_ = 0.;
0342     float score_ = 0.;
0343     short int lastHitIdx_ = -1;
0344     short int nFoundHits_ = 0;
0345     Status status_;
0346     int label_ = -1;
0347   };
0348 
0349   //==============================================================================
0350   // TrackCand
0351   //==============================================================================
0352 
0353   // TrackCand defined in TrackStructures.h along with CombCandidate.
0354 
0355   // class TrackCand : public TrackBase { ... };
0356 
0357   //==============================================================================
0358   // Track
0359   //==============================================================================
0360 
0361   class Track : public TrackBase {
0362   public:
0363     Track() {}
0364 
0365     explicit Track(const TrackBase& base) : TrackBase(base) {
0366       // Reset hit counters -- caller has to initialize hits.
0367       lastHitIdx_ = -1;
0368       nFoundHits_ = 0;
0369     }
0370 
0371     Track(const TrackState& state, float chi2, int label, int nHits, const HitOnTrack* hits)
0372         : TrackBase(state, chi2, label) {
0373       reserveHits(nHits);
0374       for (int h = 0; h < nHits; ++h) {
0375         addHitIdx(hits[h].index, hits[h].layer, 0.0f);
0376       }
0377     }
0378 
0379     Track(int charge, const SVector3& position, const SVector3& momentum, const SMatrixSym66& errors, float chi2)
0380         : TrackBase(charge, position, momentum, errors, chi2) {}
0381 
0382     // This function is very inefficient, use only for debug and validation!
0383     HitVec hitsVector(const std::vector<HitVec>& globalHitVec) const {
0384       HitVec hitsVec;
0385       for (int ihit = 0; ihit < Config::nMaxTrkHits; ++ihit) {
0386         const HitOnTrack& hot = hitsOnTrk_[ihit];
0387         if (hot.index >= 0) {
0388           hitsVec.push_back(globalHitVec[hot.layer][hot.index]);
0389         }
0390       }
0391       return hitsVec;
0392     }
0393 
0394     void mcHitIDsVec(const std::vector<HitVec>& globalHitVec,
0395                      const MCHitInfoVec& globalMCHitInfo,
0396                      std::vector<int>& mcHitIDs) const {
0397       for (int ihit = 0; ihit <= lastHitIdx_; ++ihit) {
0398         const HitOnTrack& hot = hitsOnTrk_[ihit];
0399         if ((hot.index >= 0) && (static_cast<size_t>(hot.index) < globalHitVec[hot.layer].size())) {
0400           mcHitIDs.push_back(globalHitVec[hot.layer][hot.index].mcTrackID(globalMCHitInfo));
0401         } else {
0402           mcHitIDs.push_back(hot.index);
0403         }
0404       }
0405     }
0406 
0407     int mcHitIDofFirstHit(const std::vector<HitVec>& globalHitVec, const MCHitInfoVec& globalMCHitInfo) const {
0408       const HitOnTrack& hot = hitsOnTrk_[0];
0409       if ((hot.index >= 0) && (static_cast<size_t>(hot.index) < globalHitVec[hot.layer].size())) {
0410         return globalHitVec[hot.layer][hot.index].mcTrackID(globalMCHitInfo);
0411       } else {
0412         return hot.index;
0413       }
0414     }
0415 
0416     // The following 2 (well, 3) funcs to be fixed once we move lastHitIdx_ and nFoundHits_
0417     // out of TrackBase. If we do it.
0418     void reserveHits(int nHits) { hitsOnTrk_.reserve(nHits); }
0419 
0420     void resetHits() {
0421       lastHitIdx_ = -1;
0422       nFoundHits_ = 0;
0423       hitsOnTrk_.clear();
0424     }
0425 
0426     // For MkFinder::copy_out and TrackCand::ExportTrack
0427     void resizeHits(int nHits, int nFoundHits) {
0428       hitsOnTrk_.resize(nHits);
0429       lastHitIdx_ = nHits - 1;
0430       nFoundHits_ = nFoundHits;
0431     }
0432     // Used by TrackCand::ExportTrack
0433     void setHitIdxAtPos(int pos, const HitOnTrack& hot) { hitsOnTrk_[pos] = hot; }
0434 
0435     void resizeHitsForInput();
0436 
0437     void addHitIdx(int hitIdx, int hitLyr, float chi2) {
0438       hitsOnTrk_.push_back({hitIdx, hitLyr});
0439       ++lastHitIdx_;
0440       if (hitIdx >= 0 || hitIdx == -9) {
0441         ++nFoundHits_;
0442         chi2_ += chi2;
0443       }
0444     }
0445 
0446     void addHitIdx(const HitOnTrack& hot, float chi2) { addHitIdx(hot.index, hot.layer, chi2); }
0447 
0448     HitOnTrack getHitOnTrack(int posHitIdx) const { return hitsOnTrk_[posHitIdx]; }
0449 
0450     int getHitIdx(int posHitIdx) const { return hitsOnTrk_[posHitIdx].index; }
0451     int getHitLyr(int posHitIdx) const { return hitsOnTrk_[posHitIdx].layer; }
0452 
0453     HitOnTrack getLastHitOnTrack() const { return hitsOnTrk_[lastHitIdx_]; }
0454     int getLastHitIdx() const { return hitsOnTrk_[lastHitIdx_].index; }
0455     int getLastHitLyr() const { return hitsOnTrk_[lastHitIdx_].layer; }
0456 
0457     int getLastFoundHitPos() const {
0458       int hi = lastHitIdx_;
0459       while (hi >= 0 && hitsOnTrk_[hi].index < 0)
0460         --hi;
0461       return hi;
0462     }
0463 
0464     HitOnTrack getLastFoundHitOnTrack() const {
0465       int p = getLastFoundHitPos();
0466       return p >= 0 ? hitsOnTrk_[p] : HitOnTrack(-1, -1);
0467     }
0468     int getLastFoundHitIdx() const {
0469       int p = getLastFoundHitPos();
0470       return p >= 0 ? hitsOnTrk_[p].index : -1;
0471     }
0472     int getLastFoundHitLyr() const {
0473       int p = getLastFoundHitPos();
0474       return p >= 0 ? hitsOnTrk_[p].layer : -1;
0475     }
0476 
0477     int getLastFoundMCHitID(const std::vector<HitVec>& globalHitVec) const {
0478       HitOnTrack hot = getLastFoundHitOnTrack();
0479       return globalHitVec[hot.layer][hot.index].mcHitID();
0480     }
0481 
0482     int getMCHitIDFromLayer(const std::vector<HitVec>& globalHitVec, int layer) const {
0483       int mcHitID = -1;
0484       for (int ihit = 0; ihit <= lastHitIdx_; ++ihit) {
0485         if (hitsOnTrk_[ihit].layer == layer) {
0486           mcHitID = globalHitVec[hitsOnTrk_[ihit].layer][hitsOnTrk_[ihit].index].mcHitID();
0487           break;
0488         }
0489       }
0490       return mcHitID;
0491     }
0492 
0493     const HitOnTrack* getHitsOnTrackArray() const { return hitsOnTrk_.data(); }
0494     const HitOnTrack* beginHitsOnTrack() const { return hitsOnTrk_.data(); }
0495     const HitOnTrack* endHitsOnTrack() const { return hitsOnTrk_.data() + (lastHitIdx_ + 1); }
0496 
0497     HitOnTrack* beginHitsOnTrack_nc() { return hitsOnTrk_.data(); }
0498 
0499     void setHitIdx(int posHitIdx, int newIdx) { hitsOnTrk_[posHitIdx].index = newIdx; }
0500 
0501     void setHitIdxLyr(int posHitIdx, int newIdx, int newLyr) { hitsOnTrk_[posHitIdx] = {newIdx, newLyr}; }
0502 
0503     void countAndSetNFoundHits() {
0504       nFoundHits_ = 0;
0505       for (int i = 0; i <= lastHitIdx_; i++) {
0506         if (hitsOnTrk_[i].index >= 0 || hitsOnTrk_[i].index == -9)
0507           nFoundHits_++;
0508       }
0509     }
0510 
0511     int nFoundHits() const { return nFoundHits_; }
0512     int nTotalHits() const { return lastHitIdx_ + 1; }
0513 
0514     int nOverlapHits() const { return status_.n_overlaps; }
0515     void setNOverlapHits(int n) { status_.n_overlaps = n; }
0516 
0517     int nInsideMinusOneHits() const {
0518       int n = 0;
0519       bool insideValid = false;
0520       for (int i = lastHitIdx_; i >= 0; --i) {
0521         if (hitsOnTrk_[i].index >= 0)
0522           insideValid = true;
0523         if (insideValid && hitsOnTrk_[i].index == -1)
0524           ++n;
0525       }
0526       return n;
0527     }
0528 
0529     int nTailMinusOneHits() const {
0530       int n = 0;
0531       for (int i = lastHitIdx_; i >= 0; --i) {
0532         if (hitsOnTrk_[i].index >= 0)
0533           return n;
0534         if (hitsOnTrk_[i].index == -1)
0535           ++n;
0536       }
0537       return n;
0538     }
0539 
0540     int nUniqueLayers() const {
0541       // make local copy in vector: sort it in place
0542       std::vector<HitOnTrack> tmp_hitsOnTrk(hitsOnTrk_.begin(), hitsOnTrk_.end());
0543       std::sort(tmp_hitsOnTrk.begin(), tmp_hitsOnTrk.end(), [](const auto& h1, const auto& h2) {
0544         return h1.layer < h2.layer;
0545       });
0546 
0547       // local counters
0548       auto lyr_cnt = 0;
0549       auto prev_lyr = -1;
0550 
0551       // loop over copy of hitsOnTrk
0552       for (auto ihit = 0; ihit <= lastHitIdx_; ++ihit) {
0553         const auto& hot = tmp_hitsOnTrk[ihit];
0554         const auto lyr = hot.layer;
0555         const auto idx = hot.index;
0556         if (lyr >= 0 && (idx >= 0 || idx == -9) && lyr != prev_lyr) {
0557           ++lyr_cnt;
0558           prev_lyr = lyr;
0559         }
0560       }
0561       return lyr_cnt;
0562     }
0563 
0564     // this method sorts the data member hitOnTrk_ and is ONLY to be used by sim track seeding
0565     void sortHitsByLayer();
0566 
0567     // used by fittest only (NOT mplex)
0568     std::vector<int> foundLayers() const {
0569       std::vector<int> layers;
0570       for (int ihit = 0; ihit <= lastHitIdx_; ++ihit) {
0571         if (hitsOnTrk_[ihit].index >= 0 || hitsOnTrk_[ihit].index == -9) {
0572           layers.push_back(hitsOnTrk_[ihit].layer);
0573         }
0574       }
0575       return layers;
0576     }
0577 
0578   private:
0579     std::vector<HitOnTrack> hitsOnTrk_;
0580   };
0581 
0582   typedef std::vector<Track> TrackVec;
0583   typedef std::vector<TrackVec> TrackVecVec;
0584 
0585   inline bool sortByHitsChi2(const Track& cand1, const Track& cand2) {
0586     if (cand1.nFoundHits() == cand2.nFoundHits())
0587       return cand1.chi2() < cand2.chi2();
0588     return cand1.nFoundHits() > cand2.nFoundHits();
0589   }
0590 
0591   inline bool sortByScoreCand(const Track& cand1, const Track& cand2) { return cand1.score() > cand2.score(); }
0592 
0593   inline bool sortByScoreStruct(const IdxChi2List& cand1, const IdxChi2List& cand2) {
0594     return cand1.score > cand2.score;
0595   }
0596 
0597   inline float getScoreWorstPossible() {
0598     return -std::numeric_limits<float>::max();  // used for handling of best short track during finding
0599   }
0600 
0601   inline float getScoreCand(const track_score_func& score_func,
0602                             const Track& cand1,
0603                             bool penalizeTailMissHits = false,
0604                             bool inFindCandidates = false) {
0605     int nfoundhits = cand1.nFoundHits();
0606     int noverlaphits = cand1.nOverlapHits();
0607     int nmisshits = cand1.nInsideMinusOneHits();
0608     float ntailmisshits = penalizeTailMissHits ? cand1.nTailMinusOneHits() : 0;
0609     float pt = cand1.pT();
0610     float chi2 = cand1.chi2();
0611     // Do not allow for chi2<0 in score calculation
0612     if (chi2 < 0)
0613       chi2 = 0.f;
0614     return score_func(nfoundhits, ntailmisshits, noverlaphits, nmisshits, chi2, pt, inFindCandidates);
0615   }
0616 
0617   inline float getScoreStruct(const track_score_func& score_func, const IdxChi2List& cand1) {
0618     int nfoundhits = cand1.nhits;
0619     int ntailholes = cand1.ntailholes;
0620     int noverlaphits = cand1.noverlaps;
0621     int nmisshits = cand1.nholes;
0622     float pt = cand1.pt;
0623     float chi2 = cand1.chi2;
0624     // Do not allow for chi2<0 in score calculation
0625     if (chi2 < 0)
0626       chi2 = 0.f;
0627     return score_func(nfoundhits, ntailholes, noverlaphits, nmisshits, chi2, pt, true /*inFindCandidates*/);
0628   }
0629 
0630   template <typename Vector>
0631   inline void squashPhiGeneral(Vector& v) {
0632     const int i = v.kSize - 2;  // phi index
0633     v[i] = squashPhiGeneral(v[i]);
0634   }
0635 
0636   //https://github.com/cms-sw/cmssw/blob/09c3fce6626f70fd04223e7dacebf0b485f73f54/SimTracker/TrackAssociatorProducers/plugins/getChi2.cc#L23
0637   template <typename Vector, typename Matrix>
0638   float computeHelixChi2(const Vector& simV, const Vector& recoV, const Matrix& recoM, const bool diagOnly = false) {
0639     Vector diffV = recoV - simV;
0640     if (diffV.kSize > 2)
0641       squashPhiGeneral(diffV);
0642 
0643     Matrix recoM_tmp = recoM;
0644     if (diagOnly)
0645       diagonalOnly(recoM_tmp);
0646     int invFail(0);
0647     const Matrix recoMI = recoM_tmp.InverseFast(invFail);
0648 
0649     return ROOT::Math::Dot(diffV * recoMI, diffV) / (diffV.kSize - 1);
0650   }
0651 
0652   void print(const TrackState& s);
0653   void print(std::string pfx, int itrack, const Track& trk, bool print_hits = false);
0654   void print(std::string pfx, const TrackState& s);
0655 
0656 }  // end namespace mkfit
0657 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!