TrackAssociatorProducers/plugins/TrackAssociatorByChi2Impl.cc

0001 #include "TrackAssociatorByChi2Impl.h"
0002 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0003
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "DataFormats/GeometrySurface/interface/Line.h"
0006 #include "SimTracker/TrackAssociation/interface/trackAssociationChi2.h"
0007 #include "TrackingTools/PatternTools/interface/trackingParametersAtClosestApproachToBeamSpot.h"
0008
0009 using namespace edm;
0010 using namespace reco;
0011 using namespace std;
0012 using namespace track_associator;
0013
0014 RecoToSimCollection TrackAssociatorByChi2Impl::associateRecoToSim(
0015     const RefToBaseVector<Track>& tC, const RefVector<TrackingParticleCollection>& tPCH) const {
0016   const BeamSpot& bs = *beamSpot_;
0017
0018   RecoToSimCollection outputCollection(productGetter_);
0019
0020   //dereference the Refs only once and precompute params
0021   std::vector<TrackingParticle const*> tPC;
0022   std::vector<std::pair<bool, TrackBase::ParameterVector>> tpParams;
0023   tPC.reserve(tPCH.size());
0024   tpParams.reserve(tPCH.size());
0025   for (auto const& ref : tPCH) {
0026     auto const& tp = *ref;
0027     tPC.push_back(&tp);
0028
0029     int charge = tp.charge();
0030     if (charge == 0)
0031       tpParams.emplace_back(false, TrackBase::ParameterVector());
0032     else {
0033       using BVec = Basic3DVector<double>;
0034       tpParams.emplace_back(
0035           trackingParametersAtClosestApproachToBeamSpot(BVec(tp.vertex()), BVec(tp.momentum()), charge, *mF_, bs));
0036     }
0037   }
0038
0039   int tindex = 0;
0040   for (RefToBaseVector<Track>::const_iterator rt = tC.begin(); rt != tC.end(); rt++, tindex++) {
0041     LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
0042                                 << "\n"
0043                                 << "rec::Track #" << tindex << " with pt=" << (*rt)->pt() << "\n"
0044                                 << "==========================================="
0045                                 << "\n";
0046
0047     TrackBase::ParameterVector rParameters = (*rt)->parameters();
0048
0049     TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
0050     if (onlyDiagonal_) {
0051       for (unsigned int i = 0; i < 5; i++) {
0052         for (unsigned int j = 0; j < 5; j++) {
0053           if (i != j)
0054             recoTrackCovMatrix(i, j) = 0;
0055         }
0056       }
0057     }
0058
0059     recoTrackCovMatrix.Invert();
0060
0061     int tpindex = 0;
0062     for (auto tp = tPC.begin(); tp != tPC.end(); tp++, ++tpindex) {
0063       //skip tps with a very small pt
0064       //if (sqrt((*tp)->momentum().perp2())<0.5) continue;
0065       if (!tpParams[tpindex].first)
0066         continue;
0067
0068       double chi2 = trackAssociationChi2(rParameters, recoTrackCovMatrix, tpParams[tpindex].second);
0069
0070       if (chi2 < chi2cut_) {
0071         //-chi2 because the Association Map is ordered using std::greater
0072         outputCollection.insert(tC[tindex], std::make_pair(tPCH[tpindex], -chi2));
0073       }
0074     }
0075   }
0076   outputCollection.post_insert();
0077   return outputCollection;
0078 }
0079
0080 SimToRecoCollection TrackAssociatorByChi2Impl::associateSimToReco(
0081     const RefToBaseVector<Track>& tC, const RefVector<TrackingParticleCollection>& tPCH) const {
0082   const BeamSpot& bs = *beamSpot_;
0083
0084   SimToRecoCollection outputCollection(productGetter_);
0085
0086   //compute track parameters only once
0087   std::vector<TrackBase::ParameterVector> tPars;
0088   tPars.reserve(tC.size());
0089   std::vector<TrackBase::CovarianceMatrix> tCovs;
0090   tCovs.reserve(tC.size());
0091   for (auto const& ref : tC) {
0092     auto const& aTk = *ref;
0093     tPars.emplace_back(aTk.parameters());
0094
0095     TrackBase::CovarianceMatrix recoTrackCovMatrix = aTk.covariance();
0096     if (onlyDiagonal_) {
0097       for (unsigned int i = 0; i < 5; i++) {
0098         for (unsigned int j = 0; j < 5; j++) {
0099           if (i != j)
0100             recoTrackCovMatrix(i, j) = 0;
0101         }
0102       }
0103     }
0104     recoTrackCovMatrix.Invert();
0105     tCovs.emplace_back(recoTrackCovMatrix);
0106   }
0107
0108   int tpindex = 0;
0109   for (auto tp = tPCH.begin(); tp != tPCH.end(); tp++, ++tpindex) {
0110     //skip tps with a very small pt
0111     //if (sqrt(tp->momentum().perp2())<0.5) continue;
0112     auto const& aTP = **tp;
0113     int charge = aTP.charge();
0114     if (charge == 0)
0115       continue;
0116
0117     LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
0118                                 << "\n"
0119                                 << "TrackingParticle #" << tpindex << " with pt=" << sqrt(aTP.momentum().perp2())
0120                                 << "\n"
0121                                 << "==========================================="
0122                                 << "\n";
0123
0124     using BVec = Basic3DVector<double>;
0125     auto const tpBoolParams =
0126         trackingParametersAtClosestApproachToBeamSpot(BVec(aTP.vertex()), BVec(aTP.momentum()), charge, *mF_, bs);
0127     if (!tpBoolParams.first)
0128       continue;
0129
0130     for (unsigned int tindex = 0; tindex < tC.size(); tindex++) {
0131       TrackBase::ParameterVector const& rParameters = tPars[tindex];
0132       TrackBase::CovarianceMatrix const& recoTrackCovMatrix = tCovs[tindex];
0133
0134       double chi2 = trackAssociationChi2(rParameters, recoTrackCovMatrix, tpBoolParams.second);
0135
0136       if (chi2 < chi2cut_) {
0137         //-chi2 because the Association Map is ordered using std::greater
0138         outputCollection.insert(*tp, std::make_pair(tC[tindex], -chi2));
0139       }
0140     }
0141   }
0142   outputCollection.post_insert();
0143   return outputCollection;
0144 }