ChiSquaredFitBase.cc

Go to the documentation of this file.

 
 
#include "L1Trigger/TrackFindingTMTT/interface/ChiSquaredFitBase.h"
#include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
#include "L1Trigger/TrackFindingTMTT/interface/L1fittedTrack.h"
#include "L1Trigger/TrackFindingTMTT/interface/L1track3D.h"
#include "L1Trigger/TrackFindingTMTT/interface/Utility.h"
 
#include <algorithm>
#include <functional>
#include <set>
 
namespace tmtt {
 
  ChiSquaredFitBase::ChiSquaredFitBase(const Settings* settings, const uint nPar)
      : TrackFitGeneric(settings), chiSq_(0.0) {
    // Bad stub killing settings
    numFittingIterations_ = settings_->numTrackFitIterations();
    killTrackFitWorstHit_ = settings_->killTrackFitWorstHit();
    generalResidualCut_ = settings_->generalResidualCut();  // The cut used to remove bad stubs (if nStubs > minLayers)
    killingResidualCut_ = settings_->killingResidualCut();  // The cut used to kill off tracks entirely
 
    //--- These two parameters are used to check if after the fit, there are still enough stubs on the track
    minStubLayers_ = settings_->minStubLayers();
    nPar_ = nPar;
  }
 
  void ChiSquaredFitBase::calculateChiSq(const TVectorD& resids) {
    chiSq_ = 0.0;
    uint j = 0;
    for (uint i = 0; i < stubs_.size(); i++) {
      chiSq_ += resids[j] * resids[j] + resids[j + 1] * resids[j + 1];
      j = j + 2;
    }
  }
 
  void ChiSquaredFitBase::calculateDeltaChiSq(const TVectorD& delX, const TVectorD& covX) {
    for (int i = 0; i < covX.GetNrows(); i++) {
      chiSq_ -= (delX[i]) * covX[i];
    }
  }
 
  L1fittedTrack ChiSquaredFitBase::fit(const L1track3D& l1track3D) {
    qOverPt_seed_ = l1track3D.qOverPt();
    stubs_ = l1track3D.stubs();
 
    // Get cut on number of layers including variation due to dead sectors, pt dependence etc.
    minStubLayersRed_ = Utility::numLayerCut(Utility::AlgoStep::FIT,
                                             settings_,
                                             l1track3D.iPhiSec(),
                                             l1track3D.iEtaReg(),
                                             std::abs(l1track3D.qOverPt()),
                                             l1track3D.eta());
 
    TVectorD x = seed(l1track3D);
 
    for (int i = 0; i < numFittingIterations_; i++) {
      TMatrixD d = D(x);
      TMatrixD dTrans(TMatrixD::kTransposed, d);
      TMatrixD dtVinv = dTrans * Vinv();
      TMatrixD dtVinvTrans(TMatrixD::kTransposed, dtVinv);
      //TMatrixD M = dtVinv * d; // Must insert extra factor Vinv, due to unconventional Vinv() definition.
      TMatrixD M = dtVinv * dtVinvTrans;
      TMatrixD Minv(TMatrixD::kInverted, M);
      TVectorD resids = residuals(x);
      TVectorD deltaX = Minv * dtVinv * resids;
      x = x - deltaX;
      TVectorD covX = dTrans * Vinv() * resids;
      calculateChiSq(resids);
      calculateDeltaChiSq(deltaX, covX);
 
      if (i < numFittingIterations_ - 1) {  // Don't kill stub if will not refit.
 
        resids = residuals(x);  // update resids & largestresid_
 
        bool killWorstStub = false;
        if (killTrackFitWorstHit_) {
          if (largestresid_ > killingResidualCut_) {
            killWorstStub = true;
          } else if (largestresid_ > generalResidualCut_) {
            std::vector<Stub*> stubsTmp = stubs_;
            stubsTmp.erase(stubsTmp.begin() + ilargestresid_);
            if (Utility::countLayers(settings_, stubsTmp) >= minStubLayersRed_)
              killWorstStub = true;
          } else {
            // Get better apparent tracking performance by always killing worst stub until only 4 layers left.
            if (Utility::countLayers(settings_, stubs_) > minStubLayersRed_)
              killWorstStub = true;
          }
        }
 
        if (killWorstStub) {
          stubs_.erase(stubs_.begin() + ilargestresid_);
 
          // Reject tracks with too many killed stubs & stop iterating.
          unsigned int nLayers = Utility::countLayers(settings_, stubs_);  // Count tracker layers with stubs
          bool valid = nLayers >= minStubLayersRed_;
 
          if (not valid) {
            L1fittedTrack rejectedTrk;
            return rejectedTrk;
          }
        } else {
          break;
        }
      }
    }
 
    // Reject tracks with too many killed stubs
    unsigned int nLayers = Utility::countLayers(settings_, stubs_);  // Count tracker layers with stubs
    bool valid = nLayers >= minStubLayersRed_;
 
    if (valid) {
      const unsigned int hitPattern = 0;  // FIX: Needs setting
      const float chi2rz = 0;             // FIX: Needs setting
      return L1fittedTrack(settings_,
                           &l1track3D,
                           stubs_,
                           hitPattern,
                           x[INVR] / (settings_->invPtToInvR()),
                           0,
                           x[PHI0],
                           x[Z0],
                           x[T],
                           chiSq_,
                           chi2rz,
                           nPar_);
    } else {
      L1fittedTrack rejectedTrk;
      return rejectedTrk;
    }
  }
 
}  // namespace tmtt