RecoAlgos/src/PrimaryVertexAssignment.cc

0001 #include "CommonTools/RecoAlgos/interface/PrimaryVertexAssignment.h"
0002 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0003 #include "DataFormats/VertexReco/interface/Vertex.h"
0004 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006 #include "TrackingTools/IPTools/interface/IPTools.h"
0007 #include "FWCore/Utilities/interface/isFinite.h"
0008
0009 std::pair<int, PrimaryVertexAssignment::Quality> PrimaryVertexAssignment::chargedHadronVertex(
0010     const reco::VertexCollection& vertices,
0011     const reco::TrackRef& trackRef,
0012     const reco::Track* track,
0013     float time,
0014     float timeReso,  // <0 if timing not available for this object
0015     const edm::View<reco::Candidate>& jets,
0016     const TransientTrackBuilder& builder) const {
0017   typedef reco::VertexCollection::const_iterator IV;
0018   typedef reco::Vertex::trackRef_iterator IT;
0019
0020   int iVertex = -1;
0021   size_t index = 0;
0022   float bestweight = 0;
0023   for (auto const& vtx : vertices) {
0024     float w = vtx.trackWeight(trackRef);
0025     if (w > bestweight) {
0026       bestweight = w;
0027       iVertex = index;
0028     }
0029     index++;
0030   }
0031   return chargedHadronVertex(vertices, iVertex, track, time, timeReso, jets, builder);
0032 }
0033
0034 std::pair<int, PrimaryVertexAssignment::Quality> PrimaryVertexAssignment::chargedHadronVertex(
0035     const reco::VertexCollection& vertices,
0036     int iVertex,
0037     const reco::Track* track,
0038     float time,
0039     float timeReso,  // <0 if timing not available for this object
0040     const edm::View<reco::Candidate>& jets,
0041     const TransientTrackBuilder& builder) const {
0042   typedef reco::VertexCollection::const_iterator IV;
0043   typedef reco::Vertex::trackRef_iterator IT;
0044
0045   bool useTime = useTiming_;
0046   if (edm::isNotFinite(time) || timeReso < 1e-6) {
0047     useTime = false;
0048     time = 0.;
0049     timeReso = -1.;
0050   }
0051
0052   if (preferHighRanked_) {
0053     for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
0054       int ivtx = iv - vertices.begin();
0055       if (useVertexFit_ && (iVertex == ivtx))
0056         return {ivtx, PrimaryVertexAssignment::UsedInFit};
0057
0058       double dz = std::abs(track->dz(iv->position()));
0059       double dt = std::abs(time - iv->t());
0060
0061       bool useTimeVtx = useTime && iv->tError() > 0.;
0062
0063       if ((dz < maxDzForPrimaryAssignment_ or dz / track->dzError() < maxDzSigForPrimaryAssignment_) and
0064           (!useTimeVtx or dt / timeReso < maxDtSigForPrimaryAssignment_)) {
0065         return {ivtx, PrimaryVertexAssignment::PrimaryDz};
0066       }
0067     }
0068   }
0069
0070   double firstVertexDz = std::numeric_limits<double>::max();
0071   if (!vertices.empty())
0072     firstVertexDz = std::abs(track->dz(vertices[0].position()));
0073   // recover cases where the primary vertex is split
0074   if (useVertexFit_ && (iVertex > 0) && (iVertex <= int(fNumOfPUVtxsForCharged_)) &&
0075       firstVertexDz < fDzCutForChargedFromPUVtxs_)
0076     return {0, PrimaryVertexAssignment::PrimaryDz};
0077
0078   if (useVertexFit_ && (iVertex >= 0))
0079     return {iVertex, PrimaryVertexAssignment::UsedInFit};
0080
0081   double distmin = std::numeric_limits<double>::max();
0082   double dzmin = std::numeric_limits<double>::max();
0083   double dtmin = std::numeric_limits<double>::max();
0084   int vtxIdMinSignif = -1;
0085   for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
0086     double dz = std::abs(track->dz(iv->position()));
0087     double dt = std::abs(time - iv->t());
0088
0089     double dzsig = dz / track->dzError();
0090     double dist = dzsig * dzsig;
0091
0092     bool useTimeVtx = useTime && iv->tError() > 0.;
0093     if (useTime && !useTimeVtx) {
0094       dt = timeReso;
0095     }
0096
0097     if (useTime) {
0098       double dtsig = dt / timeReso;
0099
0100       dist += dtsig * dtsig;
0101     }
0102     if (dist < distmin) {
0103       distmin = dist;
0104       dzmin = dz;
0105       dtmin = dt;
0106       vtxIdMinSignif = iv - vertices.begin();
0107     }
0108   }
0109
0110   // protect high pT particles from association to pileup vertices and assign them to the first vertex
0111   if ((fPtMaxCharged_ > 0) && (vtxIdMinSignif >= 0) && (track->pt() > fPtMaxCharged_)) {
0112     iVertex = 0;
0113   } else {
0114     // first use "closest in Z" with tight cuts (targetting primary particles)
0115     const float add_cov = vtxIdMinSignif >= 0 ? vertices[vtxIdMinSignif].covariance(2, 2) : 0.f;
0116     const float dzE = sqrt(track->dzError() * track->dzError() + add_cov);
0117     if (!fOnlyUseFirstDz_) {
0118       if (vtxIdMinSignif >= 0 and
0119           (dzmin < maxDzForPrimaryAssignment_ and dzmin / dzE < maxDzSigForPrimaryAssignment_ and
0120            track->dzError() < maxDzErrorForPrimaryAssignment_) and
0121           (!useTime or dtmin / timeReso < maxDtSigForPrimaryAssignment_))
0122         iVertex = vtxIdMinSignif;
0123     } else {
0124       // consider only distances to first vertex for association (originally used in PUPPI)
0125       if ((vtxIdMinSignif >= 0) && (std::abs(track->eta()) > fEtaMinUseDz_) &&
0126           ((firstVertexDz < maxDzForPrimaryAssignment_ && firstVertexDz / dzE < maxDzSigForPrimaryAssignment_ &&
0127             track->dzError() < maxDzErrorForPrimaryAssignment_) &&
0128            (!useTime || std::abs(time - vertices[0].t()) / timeReso < maxDtSigForPrimaryAssignment_)))
0129         iVertex = 0;
0130     }
0131   }
0132
0133   if (iVertex >= 0)
0134     return {iVertex, PrimaryVertexAssignment::PrimaryDz};
0135
0136   // if track not assigned yet, it could be a b-decay secondary , use jet axis dist criterion
0137   // first find the closest jet within maxJetDeltaR_
0138   int jetIdx = -1;
0139   double minDeltaR2 = maxJetDeltaR_ * maxJetDeltaR_;
0140   for (edm::View<reco::Candidate>::const_iterator ij = jets.begin(); ij != jets.end(); ++ij) {
0141     if (ij->pt() < minJetPt_)
0142       continue;  // skip jets below the jet Pt threshold
0143
0144     double deltaR2 = reco::deltaR2(*ij, *track);
0145     if (deltaR2 < minDeltaR2 and track->dzError() < maxDzErrorForPrimaryAssignment_) {
0146       minDeltaR2 = deltaR2;
0147       jetIdx = std::distance(jets.begin(), ij);
0148     }
0149   }
0150   // if jet found
0151   if (jetIdx != -1) {
0152     reco::TransientTrack transientTrack = builder.build(*track);
0153     GlobalVector direction(jets.at(jetIdx).px(), jets.at(jetIdx).py(), jets.at(jetIdx).pz());
0154     // find the vertex with the smallest distanceToJetAxis that is still within maxDistaneToJetAxis_
0155     int vtxIdx = -1;
0156     double minDistanceToJetAxis = maxDistanceToJetAxis_;
0157     for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
0158       // only check for vertices that are close enough in Z and for tracks that have not too high dXY
0159       if (std::abs(track->dz(iv->position())) > maxDzForJetAxisAssigment_ ||
0160           std::abs(track->dxy(iv->position())) > maxDxyForJetAxisAssigment_)
0161         continue;
0162
0163       double distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *iv).second.value();
0164       if (distanceToJetAxis < minDistanceToJetAxis) {
0165         minDistanceToJetAxis = distanceToJetAxis;
0166         vtxIdx = std::distance(vertices.begin(), iv);
0167       }
0168     }
0169     if (vtxIdx >= 0) {
0170       iVertex = vtxIdx;
0171     }
0172   }
0173   if (iVertex >= 0)
0174     return {iVertex, PrimaryVertexAssignment::BTrack};
0175
0176   // if the track is not compatible with other PVs but is compatible with the BeamSpot, we may simply have not reco'ed the PV!
0177   //  we still point it to the closest in Z, but flag it as possible orphan-primary
0178   if (!vertices.empty() && std::abs(track->dxy(vertices[0].position())) < maxDxyForNotReconstructedPrimary_ &&
0179       std::abs(track->dxy(vertices[0].position()) / track->dxyError()) < maxDxySigForNotReconstructedPrimary_)
0180     return {vtxIdMinSignif, PrimaryVertexAssignment::NotReconstructedPrimary};
0181
0182   //FIXME: here we could better handle V0s and NucInt
0183
0184   // all other tracks could be non-B secondaries and we just attach them with closest Z
0185   if (vtxIdMinSignif >= 0)
0186     return {vtxIdMinSignif, PrimaryVertexAssignment::OtherDz};
0187   //If for some reason even the dz failed (when?) we consider the track not assigned
0188   return {-1, PrimaryVertexAssignment::Unassigned};
0189 }