Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-01 02:38:49

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