VertexAssociation/src/VertexAssociatorByPositionAndTracks.cc

0001 #include <CLHEP/Units/SystemOfUnits.h>
0002 #include "SimTracker/VertexAssociation/interface/VertexAssociatorByPositionAndTracks.h"
0003 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 VertexAssociatorByPositionAndTracks::VertexAssociatorByPositionAndTracks(
0008     const edm::EDProductGetter *productGetter,
0009     double absZ,
0010     double sigmaZ,
0011     double maxRecoZ,
0012     double absT,
0013     double sigmaT,
0014     double maxRecoT,
0015     double sharedTrackFraction,
0016     const reco::RecoToSimCollection *trackRecoToSimAssociation,
0017     const reco::SimToRecoCollection *trackSimToRecoAssociation)
0018     : productGetter_(productGetter),
0019       absZ_(absZ),
0020       sigmaZ_(sigmaZ),
0021       maxRecoZ_(maxRecoZ),
0022       absT_(absT),
0023       sigmaT_(sigmaT),
0024       maxRecoT_(maxRecoT),
0025       sharedTrackFraction_(sharedTrackFraction),
0026       trackRecoToSimAssociation_(trackRecoToSimAssociation),
0027       trackSimToRecoAssociation_(trackSimToRecoAssociation) {}
0028
0029 VertexAssociatorByPositionAndTracks::VertexAssociatorByPositionAndTracks(
0030     const edm::EDProductGetter *productGetter,
0031     double absZ,
0032     double sigmaZ,
0033     double maxRecoZ,
0034     double sharedTrackFraction,
0035     const reco::RecoToSimCollection *trackRecoToSimAssociation,
0036     const reco::SimToRecoCollection *trackSimToRecoAssociation)
0037     : productGetter_(productGetter),
0038       absZ_(absZ),
0039       sigmaZ_(sigmaZ),
0040       maxRecoZ_(maxRecoZ),
0041       absT_(std::numeric_limits<double>::max()),
0042       sigmaT_(std::numeric_limits<double>::max()),
0043       maxRecoT_(std::numeric_limits<double>::max()),
0044       sharedTrackFraction_(sharedTrackFraction),
0045       trackRecoToSimAssociation_(trackRecoToSimAssociation),
0046       trackSimToRecoAssociation_(trackSimToRecoAssociation) {}
0047
0048 VertexAssociatorByPositionAndTracks::~VertexAssociatorByPositionAndTracks() {}
0049
0050 reco::VertexRecoToSimCollection VertexAssociatorByPositionAndTracks::associateRecoToSim(
0051     const edm::Handle<edm::View<reco::Vertex>> &vCH, const edm::Handle<TrackingVertexCollection> &tVCH) const {
0052   reco::VertexRecoToSimCollection ret(productGetter_);
0053
0054   const edm::View<reco::Vertex> &recoVertices = *vCH;
0055   const TrackingVertexCollection &simVertices = *tVCH;
0056
0057   LogDebug("VertexAssociation") << "VertexAssociatorByPositionAndTracks::"
0058                                    "associateRecoToSim(): associating "
0059                                 << recoVertices.size() << " reco::Vertices to" << simVertices.size()
0060                                 << " TrackingVertices";
0061
0062   // filter sim PVs
0063   std::vector<size_t> simPVindices;
0064   simPVindices.reserve(recoVertices.size());
0065   {
0066     int current_event = -1;
0067     for (size_t iSim = 0; iSim != simVertices.size(); ++iSim) {
0068       const TrackingVertex &simVertex = simVertices[iSim];
0069
0070       // Associate only to primary vertices of the in-time pileup
0071       // events (BX=0, first vertex in each of the events)
0072       if (simVertex.eventId().bunchCrossing() != 0)
0073         continue;
0074       if (simVertex.eventId().event() != current_event) {
0075         current_event = simVertex.eventId().event();
0076         simPVindices.push_back(iSim);
0077       }
0078     }
0079   }
0080
0081   for (size_t iReco = 0; iReco != recoVertices.size(); ++iReco) {
0082     const reco::Vertex &recoVertex = recoVertices[iReco];
0083
0084     // skip fake vertices
0085     if (std::abs(recoVertex.z()) > maxRecoZ_ || recoVertex.isFake() || !recoVertex.isValid() || recoVertex.ndof() < 0.)
0086       continue;
0087
0088     LogTrace("VertexAssociation") << " reco::Vertex at Z " << recoVertex.z();
0089
0090     for (const size_t iSim : simPVindices) {
0091       const TrackingVertex &simVertex = simVertices[iSim];
0092       LogTrace("VertexAssociation") << "  Considering TrackingVertex at Z " << simVertex.position().z();
0093
0094       //  recoVertex.t() == 0.  is a special value
0095       // need to change this to std::numeric_limits<double>::max() or something
0096       // more clear
0097       const bool useTiming = (absT_ != std::numeric_limits<double>::max() && recoVertex.t() != 0.);
0098       if (useTiming) {
0099         LogTrace("VertexAssociation") << " and T " << recoVertex.t() * CLHEP::second << std::endl;
0100       }
0101
0102       const double tdiff = std::abs(recoVertex.t() - simVertex.position().t() * CLHEP::second);
0103       const double zdiff = std::abs(recoVertex.z() - simVertex.position().z());
0104       if (zdiff < absZ_ && zdiff / recoVertex.zError() < sigmaZ_ &&
0105           (!useTiming || (tdiff < absT_ && tdiff / recoVertex.tError() < sigmaT_))) {
0106         auto sharedTracks = calculateVertexSharedTracks(recoVertex, simVertex, *trackRecoToSimAssociation_);
0107         auto fraction = double(sharedTracks) / recoVertex.tracksSize();
0108         if (sharedTrackFraction_ < 0 || fraction > sharedTrackFraction_) {
0109           LogTrace("VertexAssociation") << "   Matched with significance " << zdiff / recoVertex.zError() << " "
0110                                         << tdiff / recoVertex.tError() << " shared tracks " << sharedTracks
0111                                         << " reco Tracks " << recoVertex.tracksSize() << " TrackingParticles "
0112                                         << simVertex.nDaughterTracks();
0113
0114           ret.insert(reco::VertexBaseRef(vCH, iReco), std::make_pair(TrackingVertexRef(tVCH, iSim), sharedTracks));
0115         }
0116       }
0117     }
0118   }
0119
0120   ret.post_insert();
0121
0122   LogDebug("VertexAssociation") << "VertexAssociatorByPositionAndTracks::associateRecoToSim(): finished";
0123
0124   return ret;
0125 }
0126
0127 reco::VertexSimToRecoCollection VertexAssociatorByPositionAndTracks::associateSimToReco(
0128     const edm::Handle<edm::View<reco::Vertex>> &vCH, const edm::Handle<TrackingVertexCollection> &tVCH) const {
0129   reco::VertexSimToRecoCollection ret(productGetter_);
0130
0131   const edm::View<reco::Vertex> &recoVertices = *vCH;
0132   const TrackingVertexCollection &simVertices = *tVCH;
0133
0134   LogDebug("VertexAssociation") << "VertexAssociatorByPositionAndTracks::"
0135                                    "associateSimToReco(): associating "
0136                                 << simVertices.size() << " TrackingVertices to " << recoVertices.size()
0137                                 << " reco::Vertices";
0138
0139   int current_event = -1;
0140   for (size_t iSim = 0; iSim != simVertices.size(); ++iSim) {
0141     const TrackingVertex &simVertex = simVertices[iSim];
0142
0143     // Associate only primary vertices of the in-time pileup
0144     // events (BX=0, first vertex in each of the events)
0145     if (simVertex.eventId().bunchCrossing() != 0)
0146       continue;
0147     if (simVertex.eventId().event() != current_event) {
0148       current_event = simVertex.eventId().event();
0149     } else {
0150       continue;
0151     }
0152
0153     LogTrace("VertexAssociation") << " TrackingVertex at Z " << simVertex.position().z();
0154
0155     for (size_t iReco = 0; iReco != recoVertices.size(); ++iReco) {
0156       const reco::Vertex &recoVertex = recoVertices[iReco];
0157
0158       // skip fake vertices
0159       if (std::abs(recoVertex.z()) > maxRecoZ_ || recoVertex.isFake() || !recoVertex.isValid() ||
0160           recoVertex.ndof() < 0.)
0161         continue;
0162
0163       LogTrace("VertexAssociation") << "  Considering reco::Vertex at Z " << recoVertex.z();
0164       const bool useTiming = (absT_ != std::numeric_limits<double>::max() && recoVertex.t() != 0.);
0165       if (useTiming) {
0166         LogTrace("VertexAssociation") << " and T " << recoVertex.t() * CLHEP::second << std::endl;
0167       }
0168
0169       const double tdiff = std::abs(recoVertex.t() - simVertex.position().t() * CLHEP::second);
0170       const double zdiff = std::abs(recoVertex.z() - simVertex.position().z());
0171       if (zdiff < absZ_ && zdiff / recoVertex.zError() < sigmaZ_ &&
0172           (!useTiming || (tdiff < absT_ && tdiff / recoVertex.tError() < sigmaT_))) {
0173         auto sharedTracks = calculateVertexSharedTracks(simVertex, recoVertex, *trackSimToRecoAssociation_);
0174         auto fraction = double(sharedTracks) / recoVertex.tracksSize();
0175         if (sharedTrackFraction_ < 0 || fraction > sharedTrackFraction_) {
0176           LogTrace("VertexAssociation") << "   Matched with significance " << zdiff / recoVertex.zError() << " "
0177                                         << tdiff / recoVertex.tError() << " shared tracks " << sharedTracks
0178                                         << " reco Tracks " << recoVertex.tracksSize() << " TrackingParticles "
0179                                         << simVertex.nDaughterTracks();
0180
0181           ret.insert(TrackingVertexRef(tVCH, iSim), std::make_pair(reco::VertexBaseRef(vCH, iReco), sharedTracks));
0182         }
0183       }
0184     }
0185   }
0186
0187   ret.post_insert();
0188
0189   LogDebug("VertexAssociation") << "VertexAssociatorByPositionAndTracks::associateSimToReco(): finished";
0190
0191   return ret;
0192 }