File indexing completed on 2023-01-23 02:42:44
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;
0035 int hitIdx;
0036 unsigned int module;
0037 int nhits;
0038 int ntailholes;
0039 int noverlaps;
0040 int nholes;
0041 float pt;
0042 float chi2;
0043 float chi2_hit;
0044 float score;
0045 };
0046
0047
0048
0049
0050
0051 struct TrackState
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
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
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
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
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
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
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
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
0207
0208 bool not_findable : 1;
0209
0210
0211
0212
0213 bool stopped : 1;
0214
0215
0216 unsigned int prod_type : 2;
0217
0218 unsigned int align_was_seed_type : 2;
0219
0220
0221 bool duplicate : 1;
0222
0223
0224 unsigned int algorithm : 6;
0225
0226
0227 int n_overlaps : 8;
0228
0229
0230 unsigned int n_seed_hits : kNSeedHitBits;
0231
0232
0233 unsigned int eta_region : 3;
0234
0235
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
0271
0272
0273
0274
0275
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
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
0309 hltPixel = 30,
0310
0311 hltIter0 = 31,
0312 hltIter1 = 32,
0313 hltIter2 = 33,
0314 hltIter3 = 34,
0315 hltIter4 = 35,
0316
0317 hltIterX = 36,
0318
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
0336
0337
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
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366 class Track : public TrackBase {
0367 public:
0368 Track() {}
0369
0370 explicit Track(const TrackBase& base) : TrackBase(base) {
0371
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
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
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
0423
0424 void reserveHits(int nHits) { hitsOnTrk_.reserve(nHits); }
0425
0426 void resetHits() {
0427 lastHitIdx_ = -1;
0428 nFoundHits_ = 0;
0429 hitsOnTrk_.clear();
0430 }
0431
0432
0433 void resizeHits(int nHits, int nFoundHits) {
0434 hitsOnTrk_.resize(nHits);
0435 lastHitIdx_ = nHits - 1;
0436 nFoundHits_ = nFoundHits;
0437 }
0438
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
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
0554 auto lyr_cnt = 0;
0555 auto prev_lyr = -1;
0556
0557
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
0571 void sortHitsByLayer();
0572
0573
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;
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
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
0631 if (chi2 < 0)
0632 chi2 = 0.f;
0633 return score_func(nfoundhits, ntailholes, noverlaphits, nmisshits, chi2, pt, true );
0634 }
0635
0636 template <typename Vector>
0637 inline void squashPhiGeneral(Vector& v) {
0638 const int i = v.kSize - 2;
0639 v[i] = squashPhiGeneral(v[i]);
0640 }
0641
0642
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 }
0663 #endif