TrackFindingTMTT/interface/L1fittedTrack.h

0001 #ifndef L1Trigger_TrackFindingTMTT_L1fittedTrack_h
0002 #define L1Trigger_TrackFindingTMTT_L1fittedTrack_h
0003
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/L1trackBase.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
0008 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/Utility.h"
0010 #include "L1Trigger/TrackFindingTMTT/interface/TP.h"
0011 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0012 #include "L1Trigger/TrackFindingTMTT/interface/Sector.h"
0013 #include "L1Trigger/TrackFindingTMTT/interface/HTrphi.h"
0014 #include "L1Trigger/TrackFindingTMTT/interface/DigitalTrack.h"
0015 #include "L1Trigger/TrackFindingTMTT/interface/KFTrackletTrack.h"
0016
0017 #include <vector>
0018 #include <set>
0019 #include <utility>
0020 #include <string>
0021 #include <memory>
0022
0023 //=== This represents a fitted L1 track candidate found in 3 dimensions.
0024 //=== It gives access to the fitted helix parameters & chi2 etc.
0025 //=== It also calculates & gives access to associated truth particle (Tracking Particle) if any.
0026 //=== It also gives access to the 3D hough-transform track candidate (L1track3D) on which the fit was run.
0027
0028 namespace tmtt {
0029
0030   class L1fittedTrack : public L1trackBase {
0031   public:
0032     // Store a new fitted track, specifying the input Hough transform track, the stubs used for the fit,
0033     // bit-encoded hit layer pattern (numbered by increasing distance from origin),
0034     // the fitted helix parameters & chi2,
0035     // and the number of helix parameters being fitted (=5 if d0 is fitted, or =4 if d0 is not fitted).
0036     // And if track fit declared this to be a valid track (enough stubs left on track after fit etc.).
0037     L1fittedTrack(const Settings* settings,
0038                   const L1track3D* l1track3D,
0039                   const std::vector<Stub*>& stubs,
0040                   unsigned int hitPattern,
0041                   float qOverPt,
0042                   float d0,
0043                   float phi0,
0044                   float z0,
0045                   float tanLambda,
0046                   float chi2rphi,
0047                   float chi2rz,
0048                   unsigned int nHelixParam,
0049                   bool accepted = true)
0050         : L1trackBase(),
0051           settings_(settings),
0052           l1track3D_(l1track3D),
0053           stubs_(stubs),
0054           stubsConst_(stubs_.begin(), stubs_.end()),
0055           hitPattern_(hitPattern),
0056           qOverPt_(qOverPt),
0057           d0_(d0),
0058           phi0_(phi0),
0059           z0_(z0),
0060           tanLambda_(tanLambda),
0061           chi2rphi_(chi2rphi),
0062           chi2rz_(chi2rz),
0063           done_bcon_(false),
0064           qOverPt_bcon_(qOverPt),
0065           d0_bcon_(d0),
0066           phi0_bcon_(phi0),
0067           chi2rphi_bcon_(chi2rphi),
0068           nHelixParam_(nHelixParam),
0069           nSkippedLayers_(0),
0070           numUpdateCalls_(0),
0071           numIterations_(0),
0072           accepted_(accepted) {
0073       if (l1track3D != nullptr) {
0074         iPhiSec_ = l1track3D->iPhiSec();
0075         iEtaReg_ = l1track3D->iEtaReg();
0076         optoLinkID_ = l1track3D->optoLinkID();
0077       } else {  // Rejected track
0078         iPhiSec_ = 0;
0079         iEtaReg_ = 0;
0080         optoLinkID_ = 0;
0081       }
0082       if (settings != nullptr) {
0083         // Count tracker layers these stubs are in
0084         nLayers_ = Utility::countLayers(settings, stubs_);
0085         // Find associated truth particle & calculate info about match.
0086         matchedTP_ = Utility::matchingTP(settings, stubs_, nMatchedLayers_, matchedStubs_);
0087       } else {  // Rejected track
0088         nLayers_ = 0;
0089         matchedTP_ = nullptr;
0090       }
0091       // Set d0 = 0 for 4 param fit, in case fitter didn't do it.
0092       if (nHelixParam == 4) {
0093         d0_ = 0.;
0094         d0_bcon_ = 0.;
0095       }
0096       if (settings != nullptr && not settings->hybrid()) {
0097         //Sector class used to check if fitted track trajectory is in expected sector.
0098         secTmp_ = std::make_shared<Sector>(settings, iPhiSec_, iEtaReg_);
0099         // HT class used to identify HT cell that corresponds to fitted helix parameters.
0100         htRphiTmp_ = std::make_shared<HTrphi>(
0101             settings, iPhiSec_, iEtaReg_, secTmp_->etaMin(), secTmp_->etaMax(), secTmp_->phiCentre());
0102         this->setConsistentHTcell();
0103       } else {
0104         consistentCell_ = false;
0105       }
0106     }
0107
0108     // Creates track rejected by fitter.
0109     L1fittedTrack() : L1fittedTrack(nullptr, nullptr, noStubs_, 0, 0., 0., 0., 0., 0., 0., 0., 0, false) {}
0110
0111     ~L1fittedTrack() override = default;
0112
0113     //--- Optionally std::set track helix params & chi2 if beam-spot constraint is used (for 5-parameter fit).
0114     void setBeamConstr(float qOverPt_bcon, float phi0_bcon, float chi2rphi_bcon, bool accepted) {
0115       done_bcon_ = true;
0116       qOverPt_bcon_ = qOverPt_bcon;
0117       d0_bcon_ = 0.0, phi0_bcon_ = phi0_bcon;
0118       chi2rphi_bcon_ = chi2rphi_bcon;
0119       accepted_ = accepted;
0120     }
0121
0122     //--- Set/get additional info about fitted track that is specific to individual track fit algorithms (KF, LR, chi2)
0123     //--- and is used for debugging/histogramming purposes.
0124
0125     void setInfoKF(unsigned int nSkippedLayers, unsigned int numUpdateCalls) {
0126       nSkippedLayers_ = nSkippedLayers;
0127       numUpdateCalls_ = numUpdateCalls;
0128     }
0129     void setInfoLR(int numIterations, std::string lostMatchingState, std::unordered_map<std::string, int> stateCalls) {
0130       numIterations_ = numIterations;
0131       lostMatchingState_ = lostMatchingState;
0132       stateCalls_ = stateCalls;
0133     }
0134     void setInfoCHI2() {}
0135
0136     void infoKF(unsigned int& nSkippedLayers, unsigned int& numUpdateCalls) const {
0137       nSkippedLayers = nSkippedLayers_;
0138       numUpdateCalls = numUpdateCalls_;
0139     }
0140     void infoLR(int& numIterations,
0141                 std::string& lostMatchingState,
0142                 std::unordered_map<std::string, int>& stateCalls) const {
0143       numIterations = numIterations_;
0144       lostMatchingState = lostMatchingState_;
0145       stateCalls = stateCalls_;
0146     }
0147     void infoCHI2() const {}
0148
0149     //--- Convert fitted track to KFTrackletTrack format, for use with HYBRID.
0150
0151     KFTrackletTrack returnKFTrackletTrack() {
0152       KFTrackletTrack trk_(l1track3D(),
0153                            stubsConst(),
0154                            hitPattern(),
0155                            qOverPt(),
0156                            d0(),
0157                            phi0(),
0158                            z0(),
0159                            tanLambda(),
0160                            chi2rphi(),
0161                            chi2rz(),
0162                            nHelixParam(),
0163                            iPhiSec(),
0164                            iEtaReg(),
0165                            accepted(),
0166                            done_bcon(),
0167                            qOverPt_bcon(),
0168                            d0_bcon(),
0169                            phi0_bcon(),
0170                            chi2rphi_bcon());
0171       return trk_;
0172     }
0173
0174     //--- Get the 3D Hough transform track candididate corresponding to the fitted track,
0175     //--- Provide direct access to some of the info it contains.
0176
0177     // Get track candidate from HT (before fit).
0178     const L1track3D* l1track3D() const { return l1track3D_; }
0179
0180     // Get stubs on fitted track (can differ from those on HT track if track fit kicked out stubs with bad residuals)
0181     const std::vector<const Stub*>& stubsConst() const override { return stubsConst_; }
0182     const std::vector<Stub*>& stubs() const override { return stubs_; }
0183     // Get number of stubs on fitted track.
0184     unsigned int numStubs() const override { return stubs_.size(); }
0185     // Get number of tracker layers these stubs are in.
0186     unsigned int numLayers() const override { return nLayers_; }
0187     // Get number of stubs deleted from track candidate by fitter (because they had large residuals)
0188     unsigned int numKilledStubs() const { return l1track3D_->numStubs() - this->numStubs(); }
0189     // Get bit-encoded hit pattern (where layer number assigned by increasing distance from origin, according to layers track expected to cross).
0190     unsigned int hitPattern() const { return hitPattern_; }
0191
0192     // Get Hough transform cell locations in units of bin number, corresponding to the fitted helix parameters of the track.
0193     // Always uses the beam-spot constrained helix params if they are available.
0194     // (If fitted track is outside HT array, it it put in the closest bin inside it).
0195     std::pair<unsigned int, unsigned int> cellLocationFit() const { return htRphiTmp_->cell(this); }
0196     // Also get HT cell determined by Hough transform.
0197     std::pair<unsigned int, unsigned int> cellLocationHT() const override { return l1track3D_->cellLocationHT(); }
0198
0199     //--- Get information about its association (if any) to a truth Tracking Particle.
0200     //--- Can differ from that of corresponding HT track, if track fit kicked out stubs with bad residuals.
0201
0202     // Get best matching tracking particle (=nullptr if none).
0203     const TP* matchedTP() const override { return matchedTP_; }
0204     // Get the matched stubs with this Tracking Particle
0205     const std::vector<const Stub*>& matchedStubs() const override { return matchedStubs_; }
0206     // Get number of matched stubs with this Tracking Particle
0207     unsigned int numMatchedStubs() const override { return matchedStubs_.size(); }
0208     // Get number of tracker layers with matched stubs with this Tracking Particle
0209     unsigned int numMatchedLayers() const override { return nMatchedLayers_; }
0210     // Get purity of stubs on track (i.e. fraction matching best Tracking Particle)
0211     float purity() const { return numMatchedStubs() / float(numStubs()); }
0212     // Get number of stubs matched to correct TP that were deleted from track candidate by fitter.
0213     unsigned int numKilledMatchedStubs() const {
0214       unsigned int nStubCount = l1track3D_->numMatchedStubs();
0215       if (nStubCount > 0) {  // Original HT track candidate did match a truth particle
0216         const TP* tp = l1track3D_->matchedTP();
0217         for (const Stub* s : stubs_) {
0218           std::set<const TP*> assTPs = s->assocTPs();
0219           if (assTPs.find(tp) != assTPs.end())
0220             nStubCount--;  // We found a stub matched to original truth particle that survived fit.
0221         }
0222       }
0223       return nStubCount;
0224     }
0225
0226     //--- Get the fitted track helix parameters.
0227
0228     float qOverPt() const override { return qOverPt_; }
0229     float charge() const { return (qOverPt_ > 0 ? 1 : -1); }
0230     float invPt() const { return std::abs(qOverPt_); }
0231     // Protect pt against 1/pt = 0.
0232     float pt() const {
0233       constexpr float small = 1.0e-6;
0234       return 1. / (small + this->invPt());
0235     }
0236     float d0() const { return d0_; }
0237     float phi0() const override { return phi0_; }
0238     float z0() const { return z0_; }
0239     float tanLambda() const { return tanLambda_; }
0240     float theta() const { return atan2(1., tanLambda_); }  // Use atan2 to ensure 0 < theta < pi.
0241     float eta() const { return -log(tan(0.5 * this->theta())); }
0242
0243     //--- Get the fitted helix parameters with beam-spot constraint.
0244     //--- If constraint not applied (e.g. 4 param fit) then these are identical to unconstrained values.
0245
0246     bool done_bcon() const { return done_bcon_; }  // Was beam-spot constraint aplied?
0247     float qOverPt_bcon() const { return qOverPt_bcon_; }
0248     float charge_bcon() const { return (qOverPt_bcon_ > 0 ? 1 : -1); }
0249     float invPt_bcon() const { return std::abs(qOverPt_bcon_); }
0250     // Protect pt against 1/pt = 0.
0251     float pt_bcon() const {
0252       constexpr float small = 1.0e-6;
0253       return 1. / (small + this->invPt_bcon());
0254     }
0255     float phi0_bcon() const { return phi0_bcon_; }
0256     float d0_bcon() const { return d0_bcon_; }
0257
0258     // Phi and z coordinates at which track crosses "chosenR" values used by r-phi HT and rapidity sectors respectively.
0259     // (Optionally with beam-spot constraint applied).
0260     float phiAtChosenR(bool beamConstraint) const {
0261       if (beamConstraint) {
0262         return reco::deltaPhi(phi0_bcon_ - ((settings_->invPtToDphi() * settings_->chosenRofPhi()) * qOverPt_bcon_) -
0263                                   d0_bcon_ / (settings_->chosenRofPhi()),
0264                               0.);
0265       } else {
0266         return reco::deltaPhi(phi0_ - ((settings_->invPtToDphi() * settings_->chosenRofPhi()) * qOverPt_) -
0267                                   d0_ / (settings_->chosenRofPhi()),
0268                               0.);
0269       }
0270     }
0271     float zAtChosenR() const {
0272       return (z0_ + (settings_->chosenRofZ()) * tanLambda_);
0273     }  // neglects transverse impact parameter & track curvature.
0274
0275     // Get the number of helix parameters being fitted (=5 if d0 is fitted or =4 if d0 is not fitted).
0276     float nHelixParam() const { return nHelixParam_; }
0277
0278     // Get the fit degrees of freedom, chi2 & chi2/DOF (also in r-phi & r-z planes).
0279     unsigned int numDOF() const { return 2 * this->numStubs() - nHelixParam_; }
0280     unsigned int numDOFrphi() const { return this->numStubs() - (nHelixParam_ - 2); }
0281     unsigned int numDOFrz() const { return this->numStubs() - 2; }
0282     float chi2rphi() const { return chi2rphi_; }
0283     float chi2rz() const { return chi2rz_; }
0284     float chi2() const { return chi2rphi_ + chi2rz_; }
0285     float chi2dof() const { return (this->chi2()) / this->numDOF(); }
0286
0287     //--- Ditto, but if beam-spot constraint is applied.
0288     //--- If constraint not applied (e.g. 4 param fit) then these are identical to unconstrained values.
0289     unsigned int numDOF_bcon() const { return (this->numDOF() - 1); }
0290     unsigned int numDOFrphi_bcon() const { return (this->numDOFrphi() - 1); }
0291     float chi2rphi_bcon() const { return chi2rphi_bcon_; }
0292     float chi2_bcon() const { return chi2rphi_bcon_ + chi2rz_; }
0293     float chi2dof_bcon() const { return (this->chi2_bcon()) / this->numDOF_bcon(); }
0294
0295     //--- Get phi sector and eta region used by track finding code that this track is in.
0296     unsigned int iPhiSec() const override { return iPhiSec_; }
0297     unsigned int iEtaReg() const override { return iEtaReg_; }
0298
0299     //--- Opto-link ID used to send this track from HT to Track Fitter
0300     unsigned int optoLinkID() const override { return optoLinkID_; }
0301
0302     //--- Get whether the track has been rejected or accepted by the fit
0303
0304     bool accepted() const { return accepted_; }
0305
0306     //--- Functions to help eliminate duplicate tracks.
0307
0308     // Is the fitted track trajectory should lie within the same HT cell in which the track was originally found?
0309     bool consistentHTcell() const { return consistentCell_; }
0310
0311     // Determine if the fitted track trajectory should lie within the same HT cell in which the track was originally found?
0312     void setConsistentHTcell() {
0313       // Use helix params with beam-spot constaint if done in case of 5 param fit.
0314
0315       std::pair<unsigned int, unsigned int> htCell = this->cellLocationHT();
0316       bool consistent = (htCell == this->cellLocationFit());
0317
0318       if (l1track3D_->mergedHTcell()) {
0319         // If this is a merged cell, check other elements of merged cell.
0320         std::pair<unsigned int, unsigned int> htCell10(htCell.first + 1, htCell.second);
0321         std::pair<unsigned int, unsigned int> htCell01(htCell.first, htCell.second + 1);
0322         std::pair<unsigned int, unsigned int> htCell11(htCell.first + 1, htCell.second + 1);
0323         if (htCell10 == this->cellLocationFit())
0324           consistent = true;
0325         if (htCell01 == this->cellLocationFit())
0326           consistent = true;
0327         if (htCell11 == this->cellLocationFit())
0328           consistent = true;
0329       }
0330
0331       consistentCell_ = consistent;
0332     }
0333
0334     // Is the fitted track trajectory within the same (eta,phi) sector of the HT used to find it?
0335     bool consistentSector() const {
0336       if (settings_->hybrid()) {
0337         float phiCentre = 2. * M_PI * iPhiSec() / settings_->numPhiSectors();
0338         float sectorHalfWidth = M_PI / settings_->numPhiSectors();
0339         bool insidePhi = (std::abs(reco::deltaPhi(this->phiAtChosenR(done_bcon_), phiCentre)) < sectorHalfWidth);
0340         return insidePhi;
0341       } else {
0342         bool insidePhi = (std::abs(reco::deltaPhi(this->phiAtChosenR(done_bcon_), secTmp_->phiCentre())) <
0343                           secTmp_->sectorHalfWidth());
0344         bool insideEta =
0345             (this->zAtChosenR() > secTmp_->zAtChosenR_Min() && this->zAtChosenR() < secTmp_->zAtChosenR_Max());
0346         return (insidePhi && insideEta);
0347       }
0348     }
0349
0350     // Digitize track and degrade helix parameter resolution according to effect of digitisation.
0351     void digitizeTrack(const std::string& fitterName);
0352
0353     // Access to detailed info about digitized track. (Gets nullptr if trk not digitized)
0354     const DigitalTrack* digitaltrack() const { return digitalTrack_.get(); }
0355
0356   private:
0357     //--- Configuration parameters
0358     const Settings* settings_;
0359
0360     //--- The 3D hough-transform track candidate which was fitted.
0361     const L1track3D* l1track3D_;
0362
0363     //--- The stubs on the fitted track (can differ from those on HT track if fit kicked off stubs with bad residuals)
0364     std::vector<Stub*> stubs_;
0365     std::vector<const Stub*> stubsConst_;
0366     unsigned int nLayers_;
0367
0368     //--- Bit-encoded hit pattern (where layer number assigned by increasing distance from origin, according to layers track expected to cross).
0369     unsigned int hitPattern_;
0370
0371     //--- The fitted helix parameters and fit chi-squared.
0372     float qOverPt_;
0373     float d0_;
0374     float phi0_;
0375     float z0_;
0376     float tanLambda_;
0377     float chi2rphi_;
0378     float chi2rz_;
0379
0380     //--- Ditto with beam-spot constraint applied in case of 5-parameter fit, plus boolean to indicate
0381     bool done_bcon_;
0382     float qOverPt_bcon_;
0383     float d0_bcon_;
0384     float phi0_bcon_;
0385     float chi2rphi_bcon_;
0386
0387     //--- The number of helix parameters being fitted (=5 if d0 is fitted or =4 if d0 is not fitted).
0388     unsigned int nHelixParam_;
0389
0390     //--- Phi sector and eta region used track finding code that this track was in.
0391     unsigned int iPhiSec_;
0392     unsigned int iEtaReg_;
0393     //--- Opto-link ID from HT to Track Fitter.
0394     unsigned int optoLinkID_;
0395
0396     //--- Information about its association (if any) to a truth Tracking Particle.
0397     const TP* matchedTP_;
0398     std::vector<const Stub*> matchedStubs_;
0399     unsigned int nMatchedLayers_;
0400
0401     //--- Sector class used to check if fitted track trajectory is in same sector as HT used to find it.
0402     std::shared_ptr<Sector> secTmp_;  // shared so as to allow copy of L1fittedTrack.
0403     //--- r-phi HT class used to determine HT cell location that corresponds to fitted track helix parameters.
0404     std::shared_ptr<HTrphi> htRphiTmp_;
0405
0406     //--- Info specific to KF fitter.
0407     unsigned int nSkippedLayers_;
0408     unsigned int numUpdateCalls_;
0409     //--- Info specific to LR fitter.
0410     int numIterations_;
0411     std::string lostMatchingState_;
0412     std::unordered_map<std::string, int> stateCalls_;
0413
0414     std::shared_ptr<DigitalTrack> digitalTrack_;  // Class used to digitize track if required.
0415
0416     static const std::vector<Stub*> noStubs_;  // Empty vector used to initialize rejected tracks.
0417
0418     bool consistentCell_;
0419
0420     //--- Has the track fit declared this to be a valid track?
0421     bool accepted_;
0422   };
0423
0424 }  // namespace tmtt
0425
0426 #endif