TrackAssociatorByChi2Impl.cc

Go to the documentation of this file.

#include "TrackAssociatorByChi2Impl.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
 
#include "DataFormats/Math/interface/deltaPhi.h"
#include "DataFormats/GeometrySurface/interface/Line.h"
#include "SimTracker/TrackAssociation/interface/trackAssociationChi2.h"
#include "TrackingTools/PatternTools/interface/trackingParametersAtClosestApproachToBeamSpot.h"
 
using namespace edm;
using namespace reco;
using namespace std;
using namespace track_associator;
 
RecoToSimCollection TrackAssociatorByChi2Impl::associateRecoToSim(
    const RefToBaseVector<Track>& tC, const RefVector<TrackingParticleCollection>& tPCH) const {
  const BeamSpot& bs = *beamSpot_;
 
  RecoToSimCollection outputCollection(productGetter_);
 
  //dereference the Refs only once and precompute params
  std::vector<TrackingParticle const*> tPC;
  std::vector<std::pair<bool, TrackBase::ParameterVector>> tpParams;
  tPC.reserve(tPCH.size());
  tpParams.reserve(tPCH.size());
  for (auto const& ref : tPCH) {
    auto const& tp = *ref;
    tPC.push_back(&tp);
 
    int charge = tp.charge();
    if (charge == 0)
      tpParams.emplace_back(false, TrackBase::ParameterVector());
    else {
      using BVec = Basic3DVector<double>;
      tpParams.emplace_back(
          trackingParametersAtClosestApproachToBeamSpot(BVec(tp.vertex()), BVec(tp.momentum()), charge, *mF_, bs));
    }
  }
 
  int tindex = 0;
  for (RefToBaseVector<Track>::const_iterator rt = tC.begin(); rt != tC.end(); rt++, tindex++) {
    LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
                                << "\n"
                                << "rec::Track #" << tindex << " with pt=" << (*rt)->pt() << "\n"
                                << "==========================================="
                                << "\n";
 
    TrackBase::ParameterVector rParameters = (*rt)->parameters();
 
    TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
    if (onlyDiagonal_) {
      for (unsigned int i = 0; i < 5; i++) {
        for (unsigned int j = 0; j < 5; j++) {
          if (i != j)
            recoTrackCovMatrix(i, j) = 0;
        }
      }
    }
 
    recoTrackCovMatrix.Invert();
 
    int tpindex = 0;
    for (auto tp = tPC.begin(); tp != tPC.end(); tp++, ++tpindex) {
      //skip tps with a very small pt
      //if (sqrt((*tp)->momentum().perp2())<0.5) continue;
      if (!tpParams[tpindex].first)
        continue;
 
      double chi2 = trackAssociationChi2(rParameters, recoTrackCovMatrix, tpParams[tpindex].second);
 
      if (chi2 < chi2cut_) {
        //-chi2 because the Association Map is ordered using std::greater
        outputCollection.insert(tC[tindex], std::make_pair(tPCH[tpindex], -chi2));
      }
    }
  }
  outputCollection.post_insert();
  return outputCollection;
}
 
SimToRecoCollection TrackAssociatorByChi2Impl::associateSimToReco(
    const RefToBaseVector<Track>& tC, const RefVector<TrackingParticleCollection>& tPCH) const {
  const BeamSpot& bs = *beamSpot_;
 
  SimToRecoCollection outputCollection(productGetter_);
 
  //compute track parameters only once
  std::vector<TrackBase::ParameterVector> tPars;
  tPars.reserve(tC.size());
  std::vector<TrackBase::CovarianceMatrix> tCovs;
  tCovs.reserve(tC.size());
  for (auto const& ref : tC) {
    auto const& aTk = *ref;
    tPars.emplace_back(aTk.parameters());
 
    TrackBase::CovarianceMatrix recoTrackCovMatrix = aTk.covariance();
    if (onlyDiagonal_) {
      for (unsigned int i = 0; i < 5; i++) {
        for (unsigned int j = 0; j < 5; j++) {
          if (i != j)
            recoTrackCovMatrix(i, j) = 0;
        }
      }
    }
    recoTrackCovMatrix.Invert();
    tCovs.emplace_back(recoTrackCovMatrix);
  }
 
  int tpindex = 0;
  for (auto tp = tPCH.begin(); tp != tPCH.end(); tp++, ++tpindex) {
    //skip tps with a very small pt
    //if (sqrt(tp->momentum().perp2())<0.5) continue;
    auto const& aTP = **tp;
    int charge = aTP.charge();
    if (charge == 0)
      continue;
 
    LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
                                << "\n"
                                << "TrackingParticle #" << tpindex << " with pt=" << sqrt(aTP.momentum().perp2())
                                << "\n"
                                << "==========================================="
                                << "\n";
 
    using BVec = Basic3DVector<double>;
    auto const tpBoolParams =
        trackingParametersAtClosestApproachToBeamSpot(BVec(aTP.vertex()), BVec(aTP.momentum()), charge, *mF_, bs);
    if (!tpBoolParams.first)
      continue;
 
    for (unsigned int tindex = 0; tindex < tC.size(); tindex++) {
      TrackBase::ParameterVector const& rParameters = tPars[tindex];
      TrackBase::CovarianceMatrix const& recoTrackCovMatrix = tCovs[tindex];
 
      double chi2 = trackAssociationChi2(rParameters, recoTrackCovMatrix, tpBoolParams.second);
 
      if (chi2 < chi2cut_) {
        //-chi2 because the Association Map is ordered using std::greater
        outputCollection.insert(*tp, std::make_pair(tC[tindex], -chi2));
      }
    }
  }
  outputCollection.post_insert();
  return outputCollection;
}