MkFitCore/interface/Track.h

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