CMS 3D CMS Logo

PrimaryVertexAssignment.cc
Go to the documentation of this file.
8 
9 std::pair<int, PrimaryVertexAssignment::Quality> PrimaryVertexAssignment::chargedHadronVertex(
11  const reco::TrackRef& trackRef,
12  const reco::Track* track,
13  float time,
14  float timeReso, // <0 if timing not available for this object
16  const TransientTrackBuilder& builder) const {
17  typedef reco::VertexCollection::const_iterator IV;
19 
20  int iVertex = -1;
21  size_t index = 0;
22  float bestweight = 0;
23  for (auto const& vtx : vertices) {
24  float w = vtx.trackWeight(trackRef);
25  if (w > bestweight) {
26  bestweight = w;
27  iVertex = index;
28  }
29  index++;
30  }
31  return chargedHadronVertex(vertices, iVertex, track, time, timeReso, jets, builder);
32 }
33 
34 std::pair<int, PrimaryVertexAssignment::Quality> PrimaryVertexAssignment::chargedHadronVertex(
36  int iVertex,
37  const reco::Track* track,
38  float time,
39  float timeReso, // <0 if timing not available for this object
41  const TransientTrackBuilder& builder) const {
42  typedef reco::VertexCollection::const_iterator IV;
44 
45  bool useTime = useTiming_;
46  if (edm::isNotFinite(time) || timeReso < 1e-6) {
47  useTime = false;
48  time = 0.;
49  timeReso = -1.;
50  }
51 
52  if (preferHighRanked_) {
53  for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
54  int ivtx = iv - vertices.begin();
55  if (useVertexFit_ && (iVertex == ivtx))
57 
58  double dz = std::abs(track->dz(iv->position()));
59  double dt = std::abs(time - iv->t());
60 
61  bool useTimeVtx = useTime && iv->tError() > 0.;
62 
63  if ((dz < maxDzForPrimaryAssignment_ or dz / track->dzError() < maxDzSigForPrimaryAssignment_) and
64  (!useTimeVtx or dt / timeReso < maxDtSigForPrimaryAssignment_)) {
66  }
67  }
68  }
69 
70  double firstVertexDz = std::numeric_limits<double>::max();
71  if (!vertices.empty())
72  firstVertexDz = std::abs(track->dz(vertices[0].position()));
73  // recover cases where the primary vertex is split
74  if (useVertexFit_ && (iVertex > 0) && (iVertex <= int(fNumOfPUVtxsForCharged_)) &&
75  firstVertexDz < fDzCutForChargedFromPUVtxs_)
77 
78  if (useVertexFit_ && (iVertex >= 0))
79  return {iVertex, PrimaryVertexAssignment::UsedInFit};
80 
81  double distmin = std::numeric_limits<double>::max();
82  double dzmin = std::numeric_limits<double>::max();
83  double dtmin = std::numeric_limits<double>::max();
84  int vtxIdMinSignif = -1;
85  for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
86  double dz = std::abs(track->dz(iv->position()));
87  double dt = std::abs(time - iv->t());
88 
89  double dzsig = dz / track->dzError();
90  double dist = dzsig * dzsig;
91 
92  bool useTimeVtx = useTime && iv->tError() > 0.;
93  if (useTime && !useTimeVtx) {
94  dt = timeReso;
95  }
96 
97  if (useTime) {
98  double dtsig = dt / timeReso;
99 
100  dist += dtsig * dtsig;
101  }
102  if (dist < distmin) {
103  distmin = dist;
104  dzmin = dz;
105  dtmin = dt;
106  vtxIdMinSignif = iv - vertices.begin();
107  }
108  }
109 
110  // protect high pT particles from association to pileup vertices and assign them to the first vertex
111  if ((fPtMaxCharged_ > 0) && (vtxIdMinSignif >= 0) && (track->pt() > fPtMaxCharged_)) {
112  iVertex = 0;
113  } else {
114  // first use "closest in Z" with tight cuts (targetting primary particles)
115  const float add_cov = vtxIdMinSignif >= 0 ? vertices[vtxIdMinSignif].covariance(2, 2) : 0.f;
116  const float dzE = sqrt(track->dzError() * track->dzError() + add_cov);
117  if (!fOnlyUseFirstDz_) {
118  if (vtxIdMinSignif >= 0 and
119  (dzmin < maxDzForPrimaryAssignment_ and dzmin / dzE < maxDzSigForPrimaryAssignment_ and
120  track->dzError() < maxDzErrorForPrimaryAssignment_) and
121  (!useTime or dtmin / timeReso < maxDtSigForPrimaryAssignment_))
122  iVertex = vtxIdMinSignif;
123  } else {
124  // consider only distances to first vertex for association (originally used in PUPPI)
125  if ((vtxIdMinSignif >= 0) && (std::abs(track->eta()) > fEtaMinUseDz_) &&
126  ((firstVertexDz < maxDzForPrimaryAssignment_ && firstVertexDz / dzE < maxDzSigForPrimaryAssignment_ &&
127  track->dzError() < maxDzErrorForPrimaryAssignment_) &&
128  (!useTime || std::abs(time - vertices[0].t()) / timeReso < maxDtSigForPrimaryAssignment_)))
129  iVertex = 0;
130  }
131  }
132 
133  if (iVertex >= 0)
134  return {iVertex, PrimaryVertexAssignment::PrimaryDz};
135 
136  // if track not assigned yet, it could be a b-decay secondary , use jet axis dist criterion
137  // first find the closest jet within maxJetDeltaR_
138  int jetIdx = -1;
139  double minDeltaR2 = maxJetDeltaR_ * maxJetDeltaR_;
140  for (edm::View<reco::Candidate>::const_iterator ij = jets.begin(); ij != jets.end(); ++ij) {
141  if (ij->pt() < minJetPt_)
142  continue; // skip jets below the jet Pt threshold
143 
144  double deltaR2 = reco::deltaR2(*ij, *track);
145  if (deltaR2 < minDeltaR2 and track->dzError() < maxDzErrorForPrimaryAssignment_) {
146  minDeltaR2 = deltaR2;
147  jetIdx = std::distance(jets.begin(), ij);
148  }
149  }
150  // if jet found
151  if (jetIdx != -1) {
152  reco::TransientTrack transientTrack = builder.build(*track);
153  GlobalVector direction(jets.at(jetIdx).px(), jets.at(jetIdx).py(), jets.at(jetIdx).pz());
154  // find the vertex with the smallest distanceToJetAxis that is still within maxDistaneToJetAxis_
155  int vtxIdx = -1;
156  double minDistanceToJetAxis = maxDistanceToJetAxis_;
157  for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
158  // only check for vertices that are close enough in Z and for tracks that have not too high dXY
159  if (std::abs(track->dz(iv->position())) > maxDzForJetAxisAssigment_ ||
160  std::abs(track->dxy(iv->position())) > maxDxyForJetAxisAssigment_)
161  continue;
162 
163  double distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *iv).second.value();
164  if (distanceToJetAxis < minDistanceToJetAxis) {
165  minDistanceToJetAxis = distanceToJetAxis;
166  vtxIdx = std::distance(vertices.begin(), iv);
167  }
168  }
169  if (vtxIdx >= 0) {
170  iVertex = vtxIdx;
171  }
172  }
173  if (iVertex >= 0)
174  return {iVertex, PrimaryVertexAssignment::BTrack};
175 
176  // if the track is not compatible with other PVs but is compatible with the BeamSpot, we may simply have not reco'ed the PV!
177  // we still point it to the closest in Z, but flag it as possible orphan-primary
178  if (!vertices.empty() && std::abs(track->dxy(vertices[0].position())) < maxDxyForNotReconstructedPrimary_ &&
179  std::abs(track->dxy(vertices[0].position()) / track->dxyError()) < maxDxySigForNotReconstructedPrimary_)
180  return {vtxIdMinSignif, PrimaryVertexAssignment::NotReconstructedPrimary};
181 
182  //FIXME: here we could better handle V0s and NucInt
183 
184  // all other tracks could be non-B secondaries and we just attach them with closest Z
185  if (vtxIdMinSignif >= 0)
186  return {vtxIdMinSignif, PrimaryVertexAssignment::OtherDz};
187  //If for some reason even the dz failed (when?) we consider the track not assigned
189 }
float dt
Definition: AMPTWrapper.h:136
int32_t *__restrict__ iv
T w() const
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
std::pair< double, Measurement1D > jetTrackDistance(const reco::TransientTrack &track, const GlobalVector &direction, const reco::Vertex &vertex)
Definition: IPTools.cc:206
reco::TransientTrack build(const reco::Track *p) const
T sqrt(T t)
Definition: SSEVec.h:19
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< LinkConnSpec >::const_iterator IT
std::pair< int, PrimaryVertexAssignment::Quality > chargedHadronVertex(const reco::VertexCollection &vertices, const reco::TrackRef &trackRef, const reco::Track *track, float trackTime, float trackTimeResolution, const edm::View< reco::Candidate > &jets, const TransientTrackBuilder &builder) const
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:88
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
Definition: Vertex.h:38