TrackFindingTMTT/src/KFParamsComb.cc

0001 ///=== This is the Kalman Combinatorial Filter for 4 helix parameters track fit algorithm.
0002
0003 #include "L1Trigger/TrackFindingTMTT/interface/KFParamsComb.h"
0004 #include "L1Trigger/TrackFindingTMTT/interface/KalmanState.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007
0008 #include <array>
0009 #include <sstream>
0010
0011 using namespace std;
0012
0013 namespace tmtt {
0014
0015   /* Initialize */
0016
0017   KFParamsComb::KFParamsComb(const Settings* settings, const uint nHelixPar, const std::string& fitterName)
0018       : KFbase(settings, nHelixPar, fitterName),
0019         // Initialize cuts applied to helix states vs KF layer number of last added stub.
0020         kfLayerVsPtToler_(settings->kfLayerVsPtToler()),
0021         kfLayerVsD0Cut5_(settings->kfLayerVsD0Cut5()),
0022         kfLayerVsZ0Cut5_(settings->kfLayerVsZ0Cut5()),
0023         kfLayerVsZ0Cut4_(settings->kfLayerVsZ0Cut4()),
0024         kfLayerVsChiSq5_(settings->kfLayerVsChiSq5()),
0025         kfLayerVsChiSq4_(settings->kfLayerVsChiSq4()) {}
0026
0027   /* Helix state seed  */
0028
0029   TVectorD KFParamsComb::seedX(const L1track3D& l1track3D) const {
0030     TVectorD vecX(nHelixPar_);
0031     vecX[INV2R] = settings_->invPtToInvR() * l1track3D.qOverPt() / 2;
0032     vecX[PHI0] = reco::deltaPhi(l1track3D.phi0() - sectorPhi(), 0.);
0033     vecX[Z0] = l1track3D.z0();
0034     vecX[T] = l1track3D.tanLambda();
0035     if (nHelixPar_ == 5) {  // fit without d0 constraint
0036       vecX[D0] = l1track3D.d0();
0037     }
0038
0039     return vecX;
0040   }
0041
0042   /* Helix state seed covariance matrix */
0043
0044   TMatrixD KFParamsComb::seedC(const L1track3D& l1track3D) const {
0045     TMatrixD matC(nHelixPar_, nHelixPar_);
0046
0047     double invPtToInv2R = settings_->invPtToInvR() / 2;
0048
0049     // Assumed track seed (from HT) uncertainty in transverse impact parameter.
0050
0051     // Constants optimised by hand for TMTT algo.
0052     const float inv2Rsigma = 0.0314 * invPtToInv2R;
0053     constexpr float phi0sigma = 0.0102;
0054     constexpr float z0sigma = 5.0;
0055     constexpr float tanLsigma = 0.5;
0056     constexpr float d0Sigma = 1.0;
0057     // (z0, tanL, d0) uncertainties could be smaller for Hybrid, if seeded in PS? -- To check!
0058     // if (L1track3D.seedPS() > 0) z0sigma /= 4; ???
0059     matC[INV2R][INV2R] = pow(inv2Rsigma, 2);
0060     matC[PHI0][PHI0] = pow(phi0sigma, 2);
0061     matC[Z0][Z0] = pow(z0sigma, 2);
0062     matC[T][T] = pow(tanLsigma, 2);
0063     if (nHelixPar_ == 5) {  // fit without d0 constraint
0064       matC[D0][D0] = pow(d0Sigma, 2);
0065     }
0066     return matC;
0067   }
0068
0069   /* Stub position measurements in (phi,z) */
0070
0071   TVectorD KFParamsComb::vectorM(const Stub* stub) const {
0072     TVectorD meas(2);
0073     meas[PHI] = reco::deltaPhi(stub->phi(), sectorPhi());
0074     meas[Z] = stub->z();
0075     return meas;
0076   }
0077
0078   // Stub position resolution in (phi,z)
0079
0080   TMatrixD KFParamsComb::matrixV(const Stub* stub, const KalmanState* state) const {
0081     // Take Pt from input track candidate as more stable.
0082     double inv2R = (settings_->invPtToInvR()) * 0.5 * state->candidate().qOverPt();
0083     double inv2R2 = inv2R * inv2R;
0084
0085     double tanl = state->vectorX()(T);  // factor of 0.9 improves rejection
0086     double tanl2 = tanl * tanl;
0087
0088     double vphi(0);
0089     double vz(0);
0090     double vcorr(0);
0091
0092     double a = stub->sigmaPerp() * stub->sigmaPerp();
0093     double b = stub->sigmaPar() * stub->sigmaPar();
0094     double r2 = stub->r() * stub->r();
0095     double invr2 = 1. / r2;
0096
0097     // Scattering term scaling as 1/Pt.
0098     double sigmaScat = settings_->kalmanMultiScattTerm() / (state->candidate().pt());
0099     double sigmaScat2 = sigmaScat * sigmaScat;
0100
0101     if (stub->barrel()) {
0102       vphi = (a * invr2) + sigmaScat2;
0103
0104       if (stub->tiltedBarrel()) {
0105         // Convert uncertainty in (r,phi) to (z,phi).
0106         float scaleTilted = 1.;
0107         if (settings_->kalmanHOtilted()) {
0108           if (settings_->useApproxB()) {  // Simple firmware approximation
0109             scaleTilted = approxB(stub->z(), stub->r());
0110           } else {  // Exact C++ implementation.
0111             float tilt = stub->tiltAngle();
0112             scaleTilted = sin(tilt) + cos(tilt) * tanl;
0113           }
0114         }
0115         float scaleTilted2 = scaleTilted * scaleTilted;
0116         // This neglects the non-radial strip effect, assumed negligeable for PS.
0117         vz = b * scaleTilted2;
0118       } else {
0119         vz = b;
0120       }
0121
0122       if (settings_->kalmanHOfw()) {
0123         // Use approximation corresponding to current firmware.
0124         vz = b;
0125       }
0126
0127     } else {
0128       vphi = a * invr2 + sigmaScat2;
0129       vz = (b * tanl2);
0130
0131       if (not stub->psModule()) {  // Neglect these terms in PS
0132         double beta = 0.;
0133         // Add correlation term related to conversion of stub residuals from (r,phi) to (z,phi).
0134         if (settings_->kalmanHOprojZcorr() == 2)
0135           beta += -inv2R;
0136         // Add alpha correction for non-radial 2S endcap strips..
0137         if (settings_->kalmanHOalpha() == 2)
0138           beta += -stub->alpha();  // alpha is 0 except in endcap 2S disks
0139
0140         double beta2 = beta * beta;
0141         vphi += b * beta2;
0142         vcorr = b * (beta * tanl);
0143
0144         if (settings_->kalmanHOfw()) {
0145           // Use approximation corresponding to current firmware.
0146           vphi = (a * invr2) + (b * inv2R2) + sigmaScat2;
0147           vcorr = 0.;
0148           vz = (b * tanl2);
0149         }
0150       }
0151     }
0152
0153     TMatrixD matV(2, 2);
0154     matV(PHI, PHI) = vphi;
0155     matV(Z, Z) = vz;
0156     matV(PHI, Z) = vcorr;
0157     matV(Z, PHI) = vcorr;
0158
0159     return matV;
0160   }
0161
0162   /* The Kalman measurement matrix = derivative of helix intercept w.r.t. helix params */
0163   /* Here I always measure phi(r), and z(r) */
0164
0165   TMatrixD KFParamsComb::matrixH(const Stub* stub) const {
0166     TMatrixD matH(2, nHelixPar_);
0167     double r = stub->r();
0168     matH(PHI, INV2R) = -r;
0169     matH(PHI, PHI0) = 1;
0170     if (nHelixPar_ == 5) {  // fit without d0 constraint
0171       matH(PHI, D0) = -1. / r;
0172     }
0173     matH(Z, Z0) = 1;
0174     matH(Z, T) = r;
0175     return matH;
0176   }
0177
0178   /* Kalman helix ref point extrapolation matrix */
0179
0180   TMatrixD KFParamsComb::matrixF(const Stub* stub, const KalmanState* state) const {
0181     const TMatrixD unitMatrix(TMatrixD::kUnit, TMatrixD(nHelixPar_, nHelixPar_));
0182     return unitMatrix;
0183   }
0184
0185   /* Get physical helix params */
0186
0187   TVectorD KFParamsComb::trackParams(const KalmanState* state) const {
0188     TVectorD vecX = state->vectorX();
0189     TVectorD vecY(nHelixPar_);
0190     vecY[QOVERPT] = 2. * vecX[INV2R] / settings_->invPtToInvR();
0191     vecY[PHI0] = reco::deltaPhi(vecX[PHI0] + sectorPhi(), 0.);
0192     vecY[Z0] = vecX[Z0];
0193     vecY[T] = vecX[T];
0194     if (nHelixPar_ == 5) {  // fit without d0 constraint
0195       vecY[D0] = vecX[D0];
0196     }
0197     return vecY;
0198   }
0199
0200   /* If using 5 param helix fit, get track params with beam-spot constraint & track fit chi2 from applying it. */
0201   /* (N.B. chi2rz unchanged by constraint) */
0202
0203   TVectorD KFParamsComb::trackParams_BeamConstr(const KalmanState* state, double& chi2rphi) const {
0204     if (nHelixPar_ == 5) {  // fit without d0 constraint
0205       TVectorD vecX = state->vectorX();
0206       TMatrixD matC = state->matrixC();
0207       TVectorD vecY(nHelixPar_);
0208       double delChi2rphi = (vecX[D0] * vecX[D0]) / matC[D0][D0];
0209       chi2rphi = state->chi2rphi() + delChi2rphi;
0210       // Apply beam-spot constraint to helix params in transverse plane only, as most sensitive to it.
0211       vecX[INV2R] -= vecX[D0] * (matC[INV2R][D0] / matC[D0][D0]);
0212       vecX[PHI0] -= vecX[D0] * (matC[PHI0][D0] / matC[D0][D0]);
0213       vecX[D0] = 0.0;
0214       vecY[QOVERPT] = 2. * vecX[INV2R] / settings_->invPtToInvR();
0215       vecY[PHI0] = reco::deltaPhi(vecX[PHI0] + sectorPhi(), 0.);
0216       vecY[Z0] = vecX[Z0];
0217       vecY[T] = vecX[T];
0218       vecY[D0] = vecX[D0];
0219       return vecY;
0220     } else {
0221       return (this->trackParams(state));
0222     }
0223   }
0224
0225   /* Check if helix state passes cuts */
0226
0227   bool KFParamsComb::isGoodState(const KalmanState& state) const {
0228     // Set cut values that are different for 4 & 5 param helix fits.
0229     vector<double> kfLayerVsZ0Cut = (nHelixPar_ == 5) ? kfLayerVsZ0Cut5_ : kfLayerVsZ0Cut4_;
0230     vector<double> kfLayerVsChiSqCut = (nHelixPar_ == 5) ? kfLayerVsChiSq5_ : kfLayerVsChiSq4_;
0231
0232     unsigned nStubLayers = state.nStubLayers();
0233     bool goodState(true);
0234
0235     TVectorD vecY = trackParams(&state);
0236     double qOverPt = vecY[QOVERPT];
0237     double pt = std::abs(1 / qOverPt);
0238     double z0 = std::abs(vecY[Z0]);
0239
0240     // state parameter selections
0241
0242     if (z0 > kfLayerVsZ0Cut[nStubLayers])
0243       goodState = false;
0244     if (pt < settings_->houghMinPt() - kfLayerVsPtToler_[nStubLayers])
0245       goodState = false;
0246     if (nHelixPar_ == 5) {  // fit without d0 constraint
0247       double d0 = std::abs(state.vectorX()[D0]);
0248       if (d0 > kfLayerVsD0Cut5_[nStubLayers])
0249         goodState = false;
0250     }
0251
0252     // chi2 selection
0253
0254     double chi2scaled = state.chi2scaled();  // chi2(r-phi) scaled down to improve electron performance.
0255
0256     if (chi2scaled > kfLayerVsChiSqCut[nStubLayers])
0257       goodState = false;  // No separate pT selection needed
0258
0259     const bool countUpdateCalls = false;  // Print statement to count calls to Updator.
0260
0261     if (countUpdateCalls || (settings_->kalmanDebugLevel() >= 2 && tpa_ != nullptr) ||
0262         (settings_->kalmanDebugLevel() >= 2 && settings_->hybrid())) {
0263       std::stringstream text;
0264       text << std::fixed << std::setprecision(4);
0265       if (not goodState)
0266         text << "State veto:";
0267       if (goodState)
0268         text << "State kept:";
0269       text << " nlay=" << nStubLayers << " nskip=" << state.nSkippedLayers() << " chi2_scaled=" << chi2scaled;
0270       if (tpa_ != nullptr)
0271         text << " pt(mc)=" << tpa_->pt();
0272       text << " pt=" << pt << " q/pt=" << qOverPt << " tanL=" << vecY[T] << " z0=" << vecY[Z0]
0273            << " phi0=" << vecY[PHI0];
0274       if (nHelixPar_ == 5)  // fit without d0 constraint
0275         text << " d0=" << vecY[D0];
0276       text << " fake" << (tpa_ == nullptr);
0277       if (tpa_ != nullptr)
0278         text << " pt(mc)=" << tpa_->pt();
0279       PrintL1trk() << text.str();
0280     }
0281
0282     return goodState;
0283   }
0284
0285 }  // namespace tmtt