TrackFindingTMTT/interface/KFbase.h

0001 #ifndef L1Trigger_TrackFindingTMTT_KFbase_h
0002 #define L1Trigger_TrackFindingTMTT_KFbase_h
0003
0004 #include "L1Trigger/TrackFindingTMTT/interface/TrackFitGeneric.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
0008 #include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/KalmanState.h"
0010 #include <TMatrixD.h>
0011 #include <TVectorD.h>
0012
0013 #include <map>
0014 #include <vector>
0015 #include <fstream>
0016 #include <memory>
0017 #include <TString.h>
0018
0019 class TH1F;
0020 class TH2F;
0021
0022 ///=== This is the base class for the Kalman Combinatorial Filter track fit algorithm.
0023 ///=== All variable names & equations come from Fruhwirth KF paper
0024 ///=== http://dx.doi.org/10.1016/0168-9002%2887%2990887-4
0025
0026 namespace tmtt {
0027
0028   class TP;
0029   class KalmanState;
0030
0031   class KFbase : public TrackFitGeneric {
0032   public:
0033     enum PAR_IDS { INV2R, PHI0, T, Z0, D0 };
0034     enum PAR_IDS_VAR { QOVERPT };
0035     enum MEAS_IDS { PHI, Z };
0036
0037   public:
0038     // Initialize configuration
0039     KFbase(const Settings *settings, const uint nHelixPar, const std::string &fitterName = "", const uint nMeas = 2);
0040
0041     ~KFbase() override { this->resetStates(); }
0042
0043     KFbase(const KFbase &kf) = delete;  // Disable copy & move of this class.
0044     KFbase(KFbase &&kf) = delete;
0045     KFbase &operator=(const KFbase &kf) = delete;
0046     KFbase &operator=(KFbase &&kf) = delete;
0047
0048     // Do KF fit
0049     L1fittedTrack fit(const L1track3D &l1track3D) override;
0050
0051     static const unsigned int nKFlayer_ = 7;
0052     static const unsigned int nEta_ = 16;
0053     static const unsigned int nGPlayer_ = 7;
0054     static constexpr unsigned int invalidKFlayer_ = nKFlayer_;
0055
0056     // index across is GP encoded layer ID (where barrel layers=1,2,7,5,4,3 & endcap wheels=3,4,5,6,7 & 0 never occurs)
0057     // index down is eta reg
0058     // element.first is kalman layer when stub is from barrel, 7 is invalid
0059     // element.second is kalman layer when stub is from endcap, 7 is invalid
0060
0061     static constexpr std::pair<unsigned, unsigned> layerMap_[nEta_ / 2][nGPlayer_ + 1] = {
0062         {{7, 7}, {0, 7}, {1, 7}, {5, 7}, {4, 7}, {3, 7}, {7, 7}, {2, 7}},  // B1 B2 B3 B4 B5 B6
0063         {{7, 7}, {0, 7}, {1, 7}, {5, 7}, {4, 7}, {3, 7}, {7, 7}, {2, 7}},  // B1 B2 B3 B4 B5 B6
0064         {{7, 7}, {0, 7}, {1, 7}, {5, 7}, {4, 7}, {3, 7}, {7, 7}, {2, 7}},  // B1 B2 B3 B4 B5 B6
0065         {{7, 7}, {0, 7}, {1, 7}, {5, 7}, {4, 7}, {3, 7}, {7, 7}, {2, 7}},  // B1 B2 B3 B4 B5 B6
0066         {{7, 7}, {0, 7}, {1, 7}, {5, 4}, {4, 5}, {3, 6}, {7, 7}, {2, 7}},  // B1 B2 B3 B4(/D3) B5(/D2) B6(/D1)
0067         {{7, 7}, {0, 7}, {1, 7}, {7, 3}, {7, 4}, {2, 5}, {7, 6}, {2, 6}},  // B1 B2 B3(/D5)+B4(/D3) D1 D2 X D4
0068         {{7, 7}, {0, 7}, {1, 7}, {7, 2}, {7, 3}, {7, 4}, {7, 5}, {7, 6}},  // B1 B2 D1 D2 D3 D4 D5
0069         {{7, 7}, {0, 7}, {7, 7}, {7, 1}, {7, 2}, {7, 3}, {7, 4}, {7, 5}},  // B1 D1 D2 D3 D4 D5
0070     };
0071
0072   protected:
0073     // Do KF fit (internal)
0074     const KalmanState *doKF(const L1track3D &l1track3D, const std::vector<Stub *> &stubs, const TP *tpa);
0075
0076     // Add one stub to a helix state
0077     virtual const KalmanState *kalmanUpdate(
0078         unsigned nSkipped, unsigned layer, Stub *stub, const KalmanState *state, const TP *tp);
0079
0080     // Create a KalmanState, containing a helix state & next stub it is to be updated with.
0081     const KalmanState *mkState(const L1track3D &candidate,
0082                                unsigned nSkipped,
0083                                unsigned layer,
0084                                const KalmanState *last_state,
0085                                const TVectorD &x,
0086                                const TMatrixD &pxx,
0087                                const TMatrixD &K,
0088                                const TMatrixD &dcov,
0089                                Stub *stub,
0090                                double chi2rphi,
0091                                double chi2rz);
0092
0093     //--- Input data
0094
0095     // Seed track helix params & covariance matrix
0096     virtual TVectorD seedX(const L1track3D &l1track3D) const = 0;
0097     virtual TMatrixD seedC(const L1track3D &l1track3D) const = 0;
0098
0099     // Stub coordinate measurements & resolution
0100     virtual TVectorD vectorM(const Stub *stub) const = 0;
0101     virtual TMatrixD matrixV(const Stub *stub, const KalmanState *state) const = 0;
0102
0103     //--- KF maths matrix multiplications
0104
0105     // Derivate of helix intercept point w.r.t. helix params.
0106     virtual TMatrixD matrixH(const Stub *stub) const = 0;
0107     // Kalman helix ref point extrapolation matrix
0108     virtual TMatrixD matrixF(const Stub *stub = nullptr, const KalmanState *state = nullptr) const = 0;
0109     // Product of H*C*H(transpose) (where C = helix covariance matrix)
0110     TMatrixD matrixHCHt(const TMatrixD &h, const TMatrixD &c) const;
0111     // Get inverted Kalman R matrix: inverse(V + HCHt)
0112     TMatrixD matrixRinv(const TMatrixD &matH, const TMatrixD &matCref, const TMatrixD &matV) const;
0113     // Kalman gain matrix
0114     TMatrixD getKalmanGainMatrix(const TMatrixD &h, const TMatrixD &pxcov, const TMatrixD &covRinv) const;
0115
0116     // Residuals of stub with respect to helix.
0117     virtual TVectorD residual(const Stub *stub, const TVectorD &x, double candQoverPt) const;
0118
0119     // Update helix state & its covariance matrix with new stub
0120     void adjustState(const TMatrixD &K,
0121                      const TMatrixD &pxcov,
0122                      const TVectorD &x,
0123                      const TMatrixD &h,
0124                      const TVectorD &delta,
0125                      TVectorD &new_x,
0126                      TMatrixD &new_xcov) const;
0127     // Update track fit chi2 with new stub
0128     virtual void adjustChi2(const KalmanState *state,
0129                             const TMatrixD &covRinv,
0130                             const TVectorD &delta,
0131                             double &chi2rphi,
0132                             double &chi2rz) const;
0133
0134     //--- Utilities
0135
0136     // Reset internal data ready for next track.
0137     void resetStates();
0138
0139     // Convert to physical helix params instead of local ones used by KF
0140     virtual TVectorD trackParams(const KalmanState *state) const = 0;
0141     // Ditto after applying beam-spot constraint.
0142     virtual TVectorD trackParams_BeamConstr(const KalmanState *state, double &chi2rphi_bcon) const = 0;
0143
0144     // Get phi of centre of sector containing track.
0145     double sectorPhi() const {
0146       float phiCentreSec0 = -M_PI / float(settings_->numPhiNonants()) + M_PI / float(settings_->numPhiSectors());
0147       return 2. * M_PI * float(iPhiSec_) / float(settings_->numPhiSectors()) + phiCentreSec0;
0148     }
0149
0150     // Get KF layer (which is integer representing order layers cross)
0151     virtual unsigned int kalmanLayer(
0152         unsigned int iEtaReg, unsigned int layerIDreduced, bool barrel, float r, float z) const;
0153     // Check if it is unclear whether a particle is expect to cross this layer.
0154     virtual bool kalmanAmbiguousLayer(unsigned int iEtaReg, unsigned int kfLayer);
0155     // KF algo mods to cope with dead tracker layers.
0156     std::set<unsigned> kalmanDeadLayers(bool &remove2PSCut) const;
0157
0158     // Function to calculate approximation for tilted barrel modules (aka B) copied from Stub class.
0159     float approxB(float z, float r) const;
0160
0161     // Is this HLS code?
0162     virtual bool isHLS() { return false; };
0163
0164     // Helix state pases cuts.
0165     virtual bool isGoodState(const KalmanState &state) const = 0;
0166
0167     //--- Debug printout
0168     void printTP(const TP *tp) const;
0169     void printStubLayers(const std::vector<Stub *> &stubs, unsigned int iEtaReg) const;
0170     void printStub(const Stub *stub) const;
0171     void printStubs(const std::vector<Stub *> &stubs) const;
0172
0173   protected:
0174     unsigned nHelixPar_;
0175     unsigned nMeas_;
0176     unsigned numEtaRegions_;
0177
0178     unsigned int iPhiSec_;
0179     unsigned int iEtaReg_;
0180
0181     unsigned int numUpdateCalls_;
0182
0183     // All helix states KF produces for current track.
0184     std::vector<std::unique_ptr<const KalmanState>> listAllStates_;
0185
0186     const TP *tpa_;
0187   };
0188
0189 }  // namespace tmtt
0190
0191 #endif