CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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.at(0).position()));
73  // recover cases where the primary vertex is split
74  if (useVertexFit_ && (iVertex > 0) && (iVertex <= 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
121  (!useTime or dtmin / timeReso < maxDtSigForPrimaryAssignment_))
122  iVertex = vtxIdMinSignif;
123  } else {
124  // consider only distances to first vertex for association of pileup vertices (originally used in PUPPI)
125  if ((vtxIdMinSignif >= 0) && (std::abs(track->eta()) > fEtaMinUseDz_))
126  iVertex =
127  ((firstVertexDz < maxDzForPrimaryAssignment_ and firstVertexDz / dzE < maxDzSigForPrimaryAssignment_ and
129  (!useTime or std::abs(time - vertices.at(0).t()) / timeReso < maxDtSigForPrimaryAssignment_))
130  ? 0
131  : vtxIdMinSignif;
132  }
133  }
134 
135  if (iVertex >= 0)
136  return {iVertex, PrimaryVertexAssignment::PrimaryDz};
137 
138  // if track not assigned yet, it could be a b-decay secondary , use jet axis dist criterion
139  // first find the closest jet within maxJetDeltaR_
140  int jetIdx = -1;
141  double minDeltaR2 = maxJetDeltaR_ * maxJetDeltaR_;
142  for (edm::View<reco::Candidate>::const_iterator ij = jets.begin(); ij != jets.end(); ++ij) {
143  if (ij->pt() < minJetPt_)
144  continue; // skip jets below the jet Pt threshold
145 
146  double deltaR2 = reco::deltaR2(*ij, *track);
147  if (deltaR2 < minDeltaR2 and track->dzError() < maxDzErrorForPrimaryAssignment_) {
148  minDeltaR2 = deltaR2;
149  jetIdx = std::distance(jets.begin(), ij);
150  }
151  }
152  // if jet found
153  if (jetIdx != -1) {
154  reco::TransientTrack transientTrack = builder.build(*track);
155  GlobalVector direction(jets.at(jetIdx).px(), jets.at(jetIdx).py(), jets.at(jetIdx).pz());
156  // find the vertex with the smallest distanceToJetAxis that is still within maxDistaneToJetAxis_
157  int vtxIdx = -1;
158  double minDistanceToJetAxis = maxDistanceToJetAxis_;
159  for (IV iv = vertices.begin(); iv != vertices.end(); ++iv) {
160  // only check for vertices that are close enough in Z and for tracks that have not too high dXY
161  if (std::abs(track->dz(iv->position())) > maxDzForJetAxisAssigment_ ||
162  std::abs(track->dxy(iv->position())) > maxDxyForJetAxisAssigment_)
163  continue;
164 
165  double distanceToJetAxis = IPTools::jetTrackDistance(transientTrack, direction, *iv).second.value();
166  if (distanceToJetAxis < minDistanceToJetAxis) {
167  minDistanceToJetAxis = distanceToJetAxis;
168  vtxIdx = std::distance(vertices.begin(), iv);
169  }
170  }
171  if (vtxIdx >= 0) {
172  iVertex = vtxIdx;
173  }
174  }
175  if (iVertex >= 0)
176  return {iVertex, PrimaryVertexAssignment::BTrack};
177 
178  // if the track is not compatible with other PVs but is compatible with the BeamSpot, we may simply have not reco'ed the PV!
179  // we still point it to the closest in Z, but flag it as possible orphan-primary
180  if (!vertices.empty() && std::abs(track->dxy(vertices[0].position())) < maxDxyForNotReconstructedPrimary_ &&
181  std::abs(track->dxy(vertices[0].position()) / track->dxyError()) < maxDxySigForNotReconstructedPrimary_)
182  return {vtxIdMinSignif, PrimaryVertexAssignment::NotReconstructedPrimary};
183 
184  //FIXME: here we could better handle V0s and NucInt
185 
186  // all other tracks could be non-B secondaries and we just attach them with closest Z
187  if (vtxIdMinSignif >= 0)
188  return {vtxIdMinSignif, PrimaryVertexAssignment::OtherDz};
189  //If for some reason even the dz failed (when?) we consider the track not assigned
191 }
float dt
Definition: AMPTWrapper.h:136
int32_t *__restrict__ iv
const double w
Definition: UKUtility.cc:23
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::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
double dxyError() const
error on dxy
Definition: TrackBase.h:769
reco::TransientTrack build(const reco::Track *p) const
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
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
const_iterator begin() const
T sqrt(T t)
Definition: SSEVec.h:19
vector< PseudoJet > jets
double pt() const
track transverse momentum
Definition: TrackBase.h:637
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
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< LinkConnSpec >::const_iterator IT
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:622
double dzError() const
error on dz
Definition: TrackBase.h:778
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:86
const_reference at(size_type pos) const
const_iterator end() const
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:608
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector&lt;TrackRef&gt;
Definition: Vertex.h:38