TrackAssociatorProducers/plugins/TrackGenAssociatorByChi2Impl.cc

0001 #include "TrackGenAssociatorByChi2Impl.h"
0002 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0003
0004 #include "DataFormats/GeometrySurface/interface/Line.h"
0005
0006 #include "SimTracker/TrackAssociation/interface/trackAssociationChi2.h"
0007
0008 using namespace edm;
0009 using namespace reco;
0010 using namespace std;
0011
0012 double TrackGenAssociatorByChi2Impl::getChi2(const TrackBase::ParameterVector& rParameters,
0013                                              const TrackBase::CovarianceMatrix& recoTrackCovMatrix,
0014                                              const Basic3DVector<double>& momAtVtx,
0015                                              const Basic3DVector<double>& vert,
0016                                              int charge,
0017                                              const reco::BeamSpot& bs) const {
0018   return track_associator::trackAssociationChi2(rParameters, recoTrackCovMatrix, momAtVtx, vert, charge, *theMF, bs);
0019 }
0020
0021 RecoToGenCollection TrackGenAssociatorByChi2Impl::associateRecoToGen(
0022     const edm::RefToBaseVector<reco::Track>& tC, const edm::RefVector<GenParticleCollection>& tPCH) const {
0023   reco::BeamSpot const& bs = *theBeamSpot;
0024
0025   RecoToGenCollection outputCollection;
0026
0027   //dereference the edm::Refs only once
0028   std::vector<GenParticle const*> tPC;
0029   tPC.reserve(tPCH.size());
0030   for (auto const& ref : tPCH) {
0031     tPC.push_back(&(*ref));
0032   }
0033
0034   int tindex = 0;
0035   for (RefToBaseVector<reco::Track>::const_iterator rt = tC.begin(); rt != tC.end(); rt++, tindex++) {
0036     LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
0037                                 << "\n"
0038                                 << "rec::Track #" << tindex << " with pt=" << (*rt)->pt() << "\n"
0039                                 << "==========================================="
0040                                 << "\n";
0041
0042     TrackBase::ParameterVector rParameters = (*rt)->parameters();
0043
0044     TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
0045     if (onlyDiagonal) {
0046       for (unsigned int i = 0; i < 5; i++) {
0047         for (unsigned int j = 0; j < 5; j++) {
0048           if (i != j)
0049             recoTrackCovMatrix(i, j) = 0;
0050         }
0051       }
0052     }
0053
0054     recoTrackCovMatrix.Invert();
0055
0056     int tpindex = 0;
0057     for (auto tp = tPC.begin(); tp != tPC.end(); tp++, ++tpindex) {
0058       //skip tps with a very small pt
0059       //if (sqrt((*tp)->momentum().perp2())<0.5) continue;
0060       int charge = (*tp)->charge();
0061       if (charge == 0)
0062         continue;
0063       Basic3DVector<double> momAtVtx((*tp)->momentum().x(), (*tp)->momentum().y(), (*tp)->momentum().z());
0064       Basic3DVector<double> vert = (Basic3DVector<double>)(*tp)->vertex();
0065
0066       double chi2 = getChi2(rParameters, recoTrackCovMatrix, momAtVtx, vert, charge, bs);
0067
0068       if (chi2 < chi2cut) {
0069         //NOTE: tPC and tPCH use the same index for the same object
0070         outputCollection.insert(
0071             tC[tindex],
0072             std::make_pair(tPCH[tpindex],
0073                            -chi2));  //-chi2 because the Association Map is ordered using std::greater
0074       }
0075     }
0076   }
0077   outputCollection.post_insert();
0078   return outputCollection;
0079 }
0080
0081 GenToRecoCollection TrackGenAssociatorByChi2Impl::associateGenToReco(
0082     const edm::RefToBaseVector<reco::Track>& tC, const edm::RefVector<GenParticleCollection>& tPCH) const {
0083   reco::BeamSpot const& bs = *theBeamSpot;
0084
0085   GenToRecoCollection outputCollection;
0086
0087   int tpindex = 0;
0088   for (auto tp = tPCH.begin(); tp != tPCH.end(); tp++, ++tpindex) {
0089     //skip tps with a very small pt
0090     //if (sqrt((*tp)->momentum().perp2())<0.5) continue;
0091     int charge = (*tp)->charge();
0092     if (charge == 0)
0093       continue;
0094
0095     LogDebug("TrackAssociator") << "=========LOOKING FOR ASSOCIATION==========="
0096                                 << "\n"
0097                                 << "TrackingParticle #" << tpindex << " with pt=" << sqrt((*tp)->momentum().perp2())
0098                                 << "\n"
0099                                 << "==========================================="
0100                                 << "\n";
0101
0102     Basic3DVector<double> momAtVtx((*tp)->momentum().x(), (*tp)->momentum().y(), (*tp)->momentum().z());
0103     Basic3DVector<double> vert((*tp)->vertex().x(), (*tp)->vertex().y(), (*tp)->vertex().z());
0104
0105     int tindex = 0;
0106     for (RefToBaseVector<reco::Track>::const_iterator rt = tC.begin(); rt != tC.end(); rt++, tindex++) {
0107       TrackBase::ParameterVector rParameters = (*rt)->parameters();
0108       TrackBase::CovarianceMatrix recoTrackCovMatrix = (*rt)->covariance();
0109       if (onlyDiagonal) {
0110         for (unsigned int i = 0; i < 5; i++) {
0111           for (unsigned int j = 0; j < 5; j++) {
0112             if (i != j)
0113               recoTrackCovMatrix(i, j) = 0;
0114           }
0115         }
0116       }
0117       recoTrackCovMatrix.Invert();
0118
0119       double chi2 = getChi2(rParameters, recoTrackCovMatrix, momAtVtx, vert, charge, bs);
0120
0121       if (chi2 < chi2cut) {
0122         outputCollection.insert(
0123             *tp,
0124             std::make_pair(tC[tindex],
0125                            -chi2));  //-chi2 because the Association Map is ordered using std::greater
0126       }
0127     }
0128   }
0129   outputCollection.post_insert();
0130   return outputCollection;
0131 }