TrackFindingTMTT/src/ChiSquaredFitBase.cc

0001 ///=== This is the base class for the linearised chi-squared track fit algorithms.
0002
0003 ///=== Written by: Sioni Summers and Alexander D. Morton
0004
0005 #include "L1Trigger/TrackFindingTMTT/interface/ChiSquaredFitBase.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
0008 #include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/Utility.h"
0010
0011 #include <algorithm>
0012 #include <functional>
0013 #include <set>
0014
0015 namespace tmtt {
0016
0017   ChiSquaredFitBase::ChiSquaredFitBase(const Settings* settings, const uint nPar)
0018       : TrackFitGeneric(settings), chiSq_(0.0) {
0019     // Bad stub killing settings
0020     numFittingIterations_ = settings_->numTrackFitIterations();
0021     killTrackFitWorstHit_ = settings_->killTrackFitWorstHit();
0022     generalResidualCut_ = settings_->generalResidualCut();  // The cut used to remove bad stubs (if nStubs > minLayers)
0023     killingResidualCut_ = settings_->killingResidualCut();  // The cut used to kill off tracks entirely
0024
0025     //--- These two parameters are used to check if after the fit, there are still enough stubs on the track
0026     minStubLayers_ = settings_->minStubLayers();
0027     nPar_ = nPar;
0028   }
0029
0030   void ChiSquaredFitBase::calculateChiSq(const TVectorD& resids) {
0031     chiSq_ = 0.0;
0032     uint j = 0;
0033     for (uint i = 0; i < stubs_.size(); i++) {
0034       chiSq_ += resids[j] * resids[j] + resids[j + 1] * resids[j + 1];
0035       j = j + 2;
0036     }
0037   }
0038
0039   void ChiSquaredFitBase::calculateDeltaChiSq(const TVectorD& delX, const TVectorD& covX) {
0040     for (int i = 0; i < covX.GetNrows(); i++) {
0041       chiSq_ -= (delX[i]) * covX[i];
0042     }
0043   }
0044
0045   L1fittedTrack ChiSquaredFitBase::fit(const L1track3D& l1track3D) {
0046     qOverPt_seed_ = l1track3D.qOverPt();
0047     stubs_ = l1track3D.stubs();
0048
0049     // Get cut on number of layers including variation due to dead sectors, pt dependence etc.
0050     minStubLayersRed_ = Utility::numLayerCut(Utility::AlgoStep::FIT,
0051                                              settings_,
0052                                              l1track3D.iPhiSec(),
0053                                              l1track3D.iEtaReg(),
0054                                              std::abs(l1track3D.qOverPt()),
0055                                              l1track3D.eta());
0056
0057     TVectorD x = seed(l1track3D);
0058
0059     for (int i = 0; i < numFittingIterations_; i++) {
0060       TMatrixD d = D(x);
0061       TMatrixD dTrans(TMatrixD::kTransposed, d);
0062       TMatrixD dtVinv = dTrans * Vinv();
0063       TMatrixD dtVinvTrans(TMatrixD::kTransposed, dtVinv);
0064       //TMatrixD M = dtVinv * d; // Must insert extra factor Vinv, due to unconventional Vinv() definition.
0065       TMatrixD M = dtVinv * dtVinvTrans;
0066       TMatrixD Minv(TMatrixD::kInverted, M);
0067       TVectorD resids = residuals(x);
0068       TVectorD deltaX = Minv * dtVinv * resids;
0069       x = x - deltaX;
0070       TVectorD covX = dTrans * Vinv() * resids;
0071       calculateChiSq(resids);
0072       calculateDeltaChiSq(deltaX, covX);
0073
0074       if (i < numFittingIterations_ - 1) {  // Don't kill stub if will not refit.
0075
0076         resids = residuals(x);  // update resids & largestresid_
0077
0078         bool killWorstStub = false;
0079         if (killTrackFitWorstHit_) {
0080           if (largestresid_ > killingResidualCut_) {
0081             killWorstStub = true;
0082           } else if (largestresid_ > generalResidualCut_) {
0083             std::vector<Stub*> stubsTmp = stubs_;
0084             stubsTmp.erase(stubsTmp.begin() + ilargestresid_);
0085             if (Utility::countLayers(settings_, stubsTmp) >= minStubLayersRed_)
0086               killWorstStub = true;
0087           } else {
0088             // Get better apparent tracking performance by always killing worst stub until only 4 layers left.
0089             if (Utility::countLayers(settings_, stubs_) > minStubLayersRed_)
0090               killWorstStub = true;
0091           }
0092         }
0093
0094         if (killWorstStub) {
0095           stubs_.erase(stubs_.begin() + ilargestresid_);
0096
0097           // Reject tracks with too many killed stubs & stop iterating.
0098           unsigned int nLayers = Utility::countLayers(settings_, stubs_);  // Count tracker layers with stubs
0099           bool valid = nLayers >= minStubLayersRed_;
0100
0101           if (not valid) {
0102             L1fittedTrack rejectedTrk;
0103             return rejectedTrk;
0104           }
0105         } else {
0106           break;
0107         }
0108       }
0109     }
0110
0111     // Reject tracks with too many killed stubs
0112     unsigned int nLayers = Utility::countLayers(settings_, stubs_);  // Count tracker layers with stubs
0113     bool valid = nLayers >= minStubLayersRed_;
0114
0115     if (valid) {
0116       const unsigned int hitPattern = 0;  // FIX: Needs setting
0117       const float chi2rz = 0;             // FIX: Needs setting
0118       return L1fittedTrack(settings_,
0119                            &l1track3D,
0120                            stubs_,
0121                            hitPattern,
0122                            x[INVR] / (settings_->invPtToInvR()),
0123                            0,
0124                            x[PHI0],
0125                            x[Z0],
0126                            x[T],
0127                            chiSq_,
0128                            chi2rz,
0129                            nPar_);
0130     } else {
0131       L1fittedTrack rejectedTrk;
0132       return rejectedTrk;
0133     }
0134   }
0135
0136 }  // namespace tmtt