Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:30:41

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   protected:
0052     // Do KF fit (internal)
0053     const KalmanState *doKF(const L1track3D &l1track3D, const std::vector<Stub *> &stubs, const TP *tpa);
0054 
0055     // Add one stub to a helix state
0056     virtual const KalmanState *kalmanUpdate(
0057         unsigned nSkipped, unsigned layer, Stub *stub, const KalmanState *state, const TP *tp);
0058 
0059     // Create a KalmanState, containing a helix state & next stub it is to be updated with.
0060     const KalmanState *mkState(const L1track3D &candidate,
0061                                unsigned nSkipped,
0062                                unsigned layer,
0063                                const KalmanState *last_state,
0064                                const TVectorD &x,
0065                                const TMatrixD &pxx,
0066                                const TMatrixD &K,
0067                                const TMatrixD &dcov,
0068                                Stub *stub,
0069                                double chi2rphi,
0070                                double chi2rz);
0071 
0072     //--- Input data
0073 
0074     // Seed track helix params & covariance matrix
0075     virtual TVectorD seedX(const L1track3D &l1track3D) const = 0;
0076     virtual TMatrixD seedC(const L1track3D &l1track3D) const = 0;
0077 
0078     // Stub coordinate measurements & resolution
0079     virtual TVectorD vectorM(const Stub *stub) const = 0;
0080     virtual TMatrixD matrixV(const Stub *stub, const KalmanState *state) const = 0;
0081 
0082     //--- KF maths matrix multiplications
0083 
0084     // Derivate of helix intercept point w.r.t. helix params.
0085     virtual TMatrixD matrixH(const Stub *stub) const = 0;
0086     // Kalman helix ref point extrapolation matrix
0087     virtual TMatrixD matrixF(const Stub *stub = nullptr, const KalmanState *state = nullptr) const = 0;
0088     // Product of H*C*H(transpose) (where C = helix covariance matrix)
0089     TMatrixD matrixHCHt(const TMatrixD &h, const TMatrixD &c) const;
0090     // Get inverted Kalman R matrix: inverse(V + HCHt)
0091     TMatrixD matrixRinv(const TMatrixD &matH, const TMatrixD &matCref, const TMatrixD &matV) const;
0092     // Kalman gain matrix
0093     TMatrixD getKalmanGainMatrix(const TMatrixD &h, const TMatrixD &pxcov, const TMatrixD &covRinv) const;
0094 
0095     // Residuals of stub with respect to helix.
0096     virtual TVectorD residual(const Stub *stub, const TVectorD &x, double candQoverPt) const;
0097 
0098     // Update helix state & its covariance matrix with new stub
0099     void adjustState(const TMatrixD &K,
0100                      const TMatrixD &pxcov,
0101                      const TVectorD &x,
0102                      const TMatrixD &h,
0103                      const TVectorD &delta,
0104                      TVectorD &new_x,
0105                      TMatrixD &new_xcov) const;
0106     // Update track fit chi2 with new stub
0107     virtual void adjustChi2(const KalmanState *state,
0108                             const TMatrixD &covRinv,
0109                             const TVectorD &delta,
0110                             double &chi2rphi,
0111                             double &chi2rz) const;
0112 
0113     //--- Utilities
0114 
0115     // Reset internal data ready for next track.
0116     void resetStates();
0117 
0118     // Convert to physical helix params instead of local ones used by KF
0119     virtual TVectorD trackParams(const KalmanState *state) const = 0;
0120     // Ditto after applying beam-spot constraint.
0121     virtual TVectorD trackParams_BeamConstr(const KalmanState *state, double &chi2rphi_bcon) const = 0;
0122 
0123     // Get phi of centre of sector containing track.
0124     double sectorPhi() const {
0125       float phiCentreSec0 = -M_PI / float(settings_->numPhiNonants()) + M_PI / float(settings_->numPhiSectors());
0126       return 2. * M_PI * float(iPhiSec_) / float(settings_->numPhiSectors()) + phiCentreSec0;
0127     }
0128 
0129     // Get KF layer (which is integer representing order layers cross)
0130     virtual unsigned int kalmanLayer(
0131         unsigned int iEtaReg, unsigned int layerIDreduced, bool barrel, float r, float z) const;
0132     // Check if it is unclear whether a particle is expect to cross this layer.
0133     virtual bool kalmanAmbiguousLayer(unsigned int iEtaReg, unsigned int kfLayer);
0134     // KF algo mods to cope with dead tracker layers.
0135     std::set<unsigned> kalmanDeadLayers(bool &remove2PSCut) const;
0136 
0137     // Function to calculate approximation for tilted barrel modules (aka B) copied from Stub class.
0138     float approxB(float z, float r) const;
0139 
0140     // Is this HLS code?
0141     virtual bool isHLS() { return false; };
0142 
0143     // Helix state pases cuts.
0144     virtual bool isGoodState(const KalmanState &state) const = 0;
0145 
0146     //--- Debug printout
0147     void printTP(const TP *tp) const;
0148     void printStubLayers(const std::vector<Stub *> &stubs, unsigned int iEtaReg) const;
0149     void printStub(const Stub *stub) const;
0150     void printStubs(const std::vector<Stub *> &stubs) const;
0151 
0152   protected:
0153     unsigned nHelixPar_;
0154     unsigned nMeas_;
0155     unsigned numEtaRegions_;
0156 
0157     unsigned int iPhiSec_;
0158     unsigned int iEtaReg_;
0159 
0160     unsigned int numUpdateCalls_;
0161 
0162     // All helix states KF produces for current track.
0163     std::vector<std::unique_ptr<const KalmanState>> listAllStates_;
0164 
0165     const TP *tpa_;
0166   };
0167 
0168 }  // namespace tmtt
0169 
0170 #endif