CMS 3D CMS Logo

RecoTauVertexAssociator.cc
Go to the documentation of this file.
2 
3 #include <functional>
4 #include <boost/foreach.hpp>
5 
14 
15 namespace reco { namespace tau {
16 
17 namespace {
18  inline const reco::Track* getTrack(const Candidate& cand)
19  {
20  const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
21  if (pfCandPtr != nullptr) {
22  if ( pfCandPtr->trackRef().isNonnull() ) return pfCandPtr->trackRef().get();
23  else if ( pfCandPtr->gsfTrackRef().isNonnull() ) return pfCandPtr->gsfTrackRef().get();
24  else return nullptr;
25  }
26  const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
27  if (packedCand != nullptr && packedCand->hasTrackDetails())
28  return &packedCand->pseudoTrack();
29 
30  return nullptr;
31  }
32 
33  inline const reco::TrackBaseRef getTrackRef(const Candidate& cand)
34  {
35  // TauReco@MiniAOD: This version does not work on top of MiniAOD, however,
36  // it is only used for non-default track-vertex associations
37  const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
38  if (pfCandPtr != nullptr) {
39  if ( pfCandPtr->trackRef().isNonnull() ) return reco::TrackBaseRef(pfCandPtr->trackRef());
40  else if ( pfCandPtr->gsfTrackRef().isNonnull() ) return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
41  else return reco::TrackBaseRef();
42  }
43 
44  return reco::TrackBaseRef();
45  }
46 }
47 
48 // Get the highest pt track in a jet.
49 // Get the KF track if it exists. Otherwise, see if it has a GSF track.
51 {
52  std::vector<CandidatePtr> chargedPFCands = pfChargedCands(jet, true);
53  if ( verbosity_ >= 1 ) {
54  std::cout << "<RecoTauVertexAssociator::getLeadTrack>:" << std::endl;
55  std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
56  std::cout << " num. chargedPFCands = " << chargedPFCands.size() << std::endl;
57  std::cout << " vxTrkFiltering = " << vxTrkFiltering_ << std::endl;
58  }
59 
60  if ( chargedPFCands.empty() ) {
61  return reco::CandidatePtr(nullptr, 0);
62  }
63 
64  std::vector<CandidatePtr> selectedPFCands;
65  if ( vxTrkFiltering_ ) {
66  selectedPFCands = qcuts_->filterCandRefs(chargedPFCands);
67  } else {
68  selectedPFCands = chargedPFCands;
69  }
70  if ( verbosity_ >= 1 ) {
71  std::cout << " num. selectedPFCands = " << selectedPFCands.size() << std::endl;
72  }
73 
75  if ( !selectedPFCands.empty() ) {
76  double leadTrackPt = 0.;
77  if ( leadingTrkOrPFCandOption_ == kFirstTrack){ leadCand=selectedPFCands[0];}
78  else {
79  for ( std::vector<CandidatePtr>::const_iterator pfCand = selectedPFCands.begin();
80  pfCand != selectedPFCands.end(); ++pfCand ) {
81  const reco::Track* track = getTrack(**pfCand);
82  double actualTrackPt = 0. , actualTrackPtError = 0.;
83  if ( track != nullptr ) {
84  actualTrackPt = track->pt();
85  actualTrackPtError = track->ptError();
86  }
87  double trackPt = 0.;
89  trackPt = actualTrackPt - 2.*actualTrackPtError;
90  } else if ( leadingTrkOrPFCandOption_ == kLeadPFCand ) {
91  trackPt = (*pfCand)->pt();
93  trackPt = std::min(actualTrackPt, (double)(*pfCand)->pt());
94  } else assert(0);
95  if ( trackPt > leadTrackPt ) {
96  leadCand = (*pfCand);
97  leadTrackPt = trackPt;
98  }
99  }
100  }
101  }
102  if ( leadCand.isNull() ) {
103  if ( recoverLeadingTrk_ ) {
104  leadCand = chargedPFCands[0];
105  } else {
106  return reco::CandidatePtr(nullptr, 0);
107  }
108  }
109  if ( verbosity_ >= 1 ) {
110  std::cout << "leadCand: Pt = " << leadCand->pt() << ", eta = " << leadCand->eta() << ", phi = " << leadCand->phi() << std::endl;
111  }
112  return leadCand;
113 }
114 
116  auto leadCand = getLeadCand(jet);
117  if(leadCand.isNull()) return nullptr;
118  const reco::Track* track = getTrack(*leadCand);
119  return track;
120 }
121 
123  auto leadCand = getLeadCand(jet);
124  if(leadCand.isNull()) return reco::TrackBaseRef();
125  return getTrackRef(*leadCand);
126 }
127 
128 namespace {
129  // Define functors which extract the relevant information from a collection of
130  // vertices.
131  double dzToTrack(const reco::VertexRef& vtx, const reco::Track* trk)
132  {
133  if ( !trk || !vtx ) {
135  }
136  return std::abs(trk->dz(vtx->position()));
137  }
138 
139  double trackWeightInVertex(const reco::VertexRef& vtx, const reco::TrackBaseRef& trk)
140  {
141  if ( !trk || !vtx ) {
142  return 0.0;
143  }
144  return vtx->trackWeight(trk);
145  }
146 }
147 
150  qcuts_(nullptr),
152  lastEvent_(0)
153 {
154  //std::cout << "<RecoTauVertexAssociator::RecoTauVertexAssociator>:" << std::endl;
155 
156  vertexTag_ = edm::InputTag("offlinePrimaryVertices", "");
157  algorithm_ = "highestPtInEvent";
158  // Sanity check, will remove once HLT module configs are updated.
159  if ( !pset.exists("primaryVertexSrc") || !pset.exists("pvFindingAlgo") ) {
160  edm::LogWarning("NoVertexFindingMethodSpecified")
161  << "The PSet passed to the RecoTauVertexAssociator was incorrectly configured."
162  << " The vertex will be taken as the highest Pt vertex from the 'offlinePrimaryVertices' collection."
163  << std::endl;
164  } else {
165  vertexTag_ = pset.getParameter<edm::InputTag>("primaryVertexSrc");
166  algorithm_ = pset.getParameter<std::string>("pvFindingAlgo");
167  }
168 
169  if ( pset.exists("vxAssocQualityCuts") ) {
170  //std::cout << " reading 'vxAssocQualityCuts'" << std::endl;
171  qcuts_ = new RecoTauQualityCuts(pset.getParameterSet("vxAssocQualityCuts"));
172  } else {
173  //std::cout << " reading 'signalQualityCuts'" << std::endl;
174  qcuts_ = new RecoTauQualityCuts(pset.getParameterSet("signalQualityCuts"));
175  }
176  assert(qcuts_);
177 
178  vxTrkFiltering_ = false;
179  if ( !pset.exists("vertexTrackFiltering") && pset.exists("vxAssocQualityCuts") ) {
180  edm::LogWarning("NoVertexTrackFilteringSpecified")
181  << "The PSet passed to the RecoTauVertexAssociator was incorrectly configured."
182  << " Please define vertexTrackFiltering in config file."
183  << " No filtering of tracks to vertices will be applied."
184  << std::endl;
185  } else {
186  vxTrkFiltering_ = pset.exists("vertexTrackFiltering") ?
187  pset.getParameter<bool>("vertexTrackFiltering") : false;
188  }
189  if ( pset.exists("vertexSelection") ) {
190  std::string vertexSelection = pset.getParameter<std::string>("vertexSelection");
191  if ( !vertexSelection.empty() ) {
193  }
194  }
195 
196  if ( algorithm_ == "highestPtInEvent" ) {
198  } else if ( algorithm_ == "closestInDeltaZ" ) {
200  } else if ( algorithm_ == "highestWeightForLeadTrack" ) {
202  } else if ( algorithm_ == "combined" ) {
203  algo_ = kCombined;
204  } else {
205  throw cms::Exception("BadVertexAssociatorConfig")
206  << "Invalid Configuration parameter 'algorithm' " << algorithm_ << "."
207  << " Valid options are: 'highestPtInEvent', 'closestInDeltaZ', 'highestWeightForLeadTrack' and 'combined'.\n";
208  }
209 
211 
212  recoverLeadingTrk_ = pset.exists("recoverLeadingTrk") ?
213  pset.getParameter<bool>("recoverLeadingTrk") : false;
214 
215  std::string leadingTrkOrPFCandOption_string = pset.exists("leadingTrkOrPFCandOption") ? pset.getParameter<std::string>("leadingTrkOrPFCandOption") : "firstTrack" ;
216  if ( leadingTrkOrPFCandOption_string == "leadTrack" ) leadingTrkOrPFCandOption_ = kLeadTrack;
217  else if ( leadingTrkOrPFCandOption_string == "leadPFCand" ) leadingTrkOrPFCandOption_ = kLeadPFCand;
218  else if ( leadingTrkOrPFCandOption_string == "minLeadTrackOrPFCand" ) leadingTrkOrPFCandOption_ = kMinLeadTrackOrPFCand;
219  else if ( leadingTrkOrPFCandOption_string == "firstTrack" ) leadingTrkOrPFCandOption_ = kFirstTrack;
220  else throw cms::Exception("BadVertexAssociatorConfig")
221  << "Invalid Configuration parameter 'leadingTrkOrPFCandOption' " << leadingTrkOrPFCandOption_string << "."
222  << " Valid options are: 'leadTrack', 'leadPFCand', 'firstTrack'.\n";
223 
224  verbosity_ = ( pset.exists("verbosity") ) ?
225  pset.getParameter<int>("verbosity") : 0;
226 }
227 
229 {
230  delete vertexSelector_;
231  delete qcuts_;
233 }
234 
236 {
238  evt.getByToken(vxToken_, vertices);
239  selectedVertices_.clear();
240  selectedVertices_.reserve(vertices->size());
241  for ( size_t idxVertex = 0; idxVertex < vertices->size(); ++idxVertex ) {
242  reco::VertexRef vertex(vertices, idxVertex);
243  if ( vertexSelector_ && !(*vertexSelector_)(*vertex) ) continue;
244  selectedVertices_.push_back(vertex);
245  }
246  if ( !selectedVertices_.empty() ) {
248  }
249  edm::EventNumber_t currentEvent = evt.id().event();
250  if ( currentEvent != lastEvent_ || !jetToVertexAssociation_ ) {
251  if ( !jetToVertexAssociation_ ) jetToVertexAssociation_ = new std::map<const reco::Jet*, reco::VertexRef>;
252  else jetToVertexAssociation_->clear();
253  lastEvent_ = currentEvent;
254  }
255 }
256 
259 {
260  if ( !useJet ) {
261  if ( tau.leadChargedHadrCand().isNonnull() ) {
263  if (track != nullptr)
264  return associatedVertex(track);
265  }
266  }
267  // MB: use vertex associated to a given jet if explicitely requested or in case of missing leading track
268  reco::JetBaseRef jetRef = tau.jetRef();
269  // FIXME workaround for HLT which does not use updated data format
270  if ( jetRef.isNull() ) jetRef = tau.pfTauTagInfoRef()->pfjetRef();
271  return associatedVertex(*jetRef);
272 }
273 
276 {
277  reco::VertexRef trkVertex = ( !selectedVertices_.empty() ) ? selectedVertices_[0] : reco::VertexRef();
278 
279  // algos kHighestWeigtForLeadTrack and kCombined not supported
280  if ( algo_ == kHighestPtInEvent ) {
281  if ( !selectedVertices_.empty() ) trkVertex = selectedVertices_[0];
282  } else if ( algo_ == kClosestDeltaZ ) {
283  if ( track ) {
284  double closestDistance = 1.e+6;
285  // Find the vertex that has the lowest dZ to the track
286  int idxVertex = 0;
287  for ( std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
288  selectedVertex != selectedVertices_.end(); ++selectedVertex ) {
289  double dZ = dzToTrack(*selectedVertex, track);
290  if ( verbosity_ ) {
291  std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x() << ", y = " << (*selectedVertex)->position().y() << ", z = " << (*selectedVertex)->position().z()
292  << " --> dZ = " << dZ << std::endl;
293  }
294  if ( dZ < closestDistance ) {
295  trkVertex = (*selectedVertex);
296  closestDistance = dZ;
297  }
298  ++idxVertex;
299  }
300  }
301  }
302 
303  if ( verbosity_ >= 1 ) {
304  std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y() << ", z = " << trkVertex->position().z() << std::endl;
305  }
306 
307  return trkVertex;
308 }
309 
310 
313 {
314 
315  reco::VertexRef trkVertex = ( !selectedVertices_.empty() ) ? selectedVertices_[0] : reco::VertexRef();
316 
318  if ( track.isNonnull() ) {
319  double largestWeight = -1.;
320  // Find the vertex that has the highest association probability to the track
321  int idxVertex = 0;
322  for ( std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
323  selectedVertex != selectedVertices_.end(); ++selectedVertex ) {
324  double weight = trackWeightInVertex(*selectedVertex, track);
325  if ( verbosity_ ) {
326  std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x() << ", y = " << (*selectedVertex)->position().y() << ", z = " << (*selectedVertex)->position().z()
327  << " --> weight = " << weight << std::endl;
328  }
329  if ( weight > largestWeight ) {
330  trkVertex = (*selectedVertex);
331  largestWeight = weight;
332  }
333  ++idxVertex;
334  }
335  // the weight was never larger than zero
336  if ( algo_ == kCombined && largestWeight < 1.e-7 ) {
337  if ( verbosity_ ) {
338  std::cout << "No vertex had positive weight! Trying dZ instead... " << std::endl;
339  }
340  double closestDistance = 1.e+6;
341  // Find the vertex that has the lowest dZ to the leading track
342  int idxVertex = 0;
343  for ( std::vector<reco::VertexRef>::const_iterator selectedVertex = selectedVertices_.begin();
344  selectedVertex != selectedVertices_.end(); ++selectedVertex ) {
345  double dZ = dzToTrack(*selectedVertex, track.get());
346  if ( verbosity_ ) {
347  std::cout << "vertex #" << idxVertex << ": x = " << (*selectedVertex)->position().x() << ", y = " << (*selectedVertex)->position().y() << ", z = " << (*selectedVertex)->position().z()
348  << " --> dZ = " << dZ << std::endl;
349  }
350  if ( dZ < closestDistance ) {
351  trkVertex = (*selectedVertex);
352  closestDistance = dZ;
353  }
354  ++idxVertex;
355  }
356  }
357  }
358  }
359 
360  if ( verbosity_ >= 1 ) {
361  std::cout << "--> returning vertex: x = " << trkVertex->position().x() << ", y = " << trkVertex->position().y() << ", z = " << trkVertex->position().z() << std::endl;
362  }
363 
364  return trkVertex;
365 }
366 
369 {
370 if ( verbosity_ >= 1 ) {
371  std::cout << "<RecoTauVertexAssociator::associatedVertex>:" << std::endl;
372  std::cout << " jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi() << std::endl;
373  std::cout << " num. Vertices = " << selectedVertices_.size() << std::endl;
374  std::cout << " size(jetToVertexAssociation) = " << jetToVertexAssociation_->size() << std::endl;
375  std::cout << " vertexTag = " << vertexTag_ << std::endl;
376  std::cout << " algorithm = " << algorithm_ << std::endl;
377  std::cout << " recoverLeadingTrk = " << recoverLeadingTrk_ << std::endl;
378  }
379 
380  reco::VertexRef jetVertex = ( !selectedVertices_.empty() ) ? selectedVertices_[0] : reco::VertexRef();
381  const Jet* jetPtr = &jet;
382 
383  // check if jet-vertex association has been determined for this jet before
384  std::map<const reco::Jet*, reco::VertexRef>::iterator vertexPtr = jetToVertexAssociation_->find(jetPtr);
385  if ( vertexPtr != jetToVertexAssociation_->end() ) {
386  jetVertex = vertexPtr->second;
387  } else {
388  // no jet-vertex association exists for this jet yet, compute it!
389  if ( algo_ == kHighestPtInEvent ) {
390  if ( !selectedVertices_.empty() ) jetVertex = selectedVertices_[0];
391  } else if ( algo_ == kClosestDeltaZ ) {
392  // find "leading" (highest Pt) track in jet
393  const reco::Track* leadTrack = getLeadTrack(jet);
394  if ( verbosity_ ) {
395  if ( leadTrack != nullptr ) std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta() << ", phi = " << leadTrack->phi() << std::endl;
396  else std::cout << "leadTrack: N/A" << std::endl;
397  }
398  if ( leadTrack != nullptr ) {
399  jetVertex = associatedVertex(leadTrack);
400  }
401  } else if (algo_ == kHighestWeigtForLeadTrack ||
402  algo_ == kCombined ) {
403  // find "leading" (highest Pt) track in jet
405  if ( verbosity_ ) {
406  if ( leadTrack.isNonnull() ) std::cout << "leadTrack: Pt = " << leadTrack->pt() << ", eta = " << leadTrack->eta() << ", phi = " << leadTrack->phi() << std::endl;
407  else std::cout << "leadTrack: N/A" << std::endl;
408  }
409  if ( leadTrack.isNonnull()) {
410  jetVertex = associatedVertex(leadTrack);
411  }
412  }
413 
414  jetToVertexAssociation_->insert(std::pair<const Jet*, reco::VertexRef>(jetPtr, jetVertex));
415  }
416 
417  if ( verbosity_ >= 1 ) {
418  std::cout << "--> returning vertex: x = " << jetVertex->position().x() << ", y = " << jetVertex->position().y() << ", z = " << jetVertex->position().z() << std::endl;
419  }
420 
421  return jetVertex;
422 }
423 
424 }}
value_type const * get() const
Definition: RefToBase.h:234
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
const TrackBaseRef getLeadTrackRef(const Jet &) const
double eta() const final
momentum pseudorapidity
Coll filterCandRefs(const Coll &refcoll, bool invert=false) const
Filter a ref vector of Candidates.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
const Track * getLeadTrack(const Jet &) const
#define nullptr
bool exists(std::string const &parameterName) const
checks if a parameter exists
unsigned long long EventNumber_t
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:684
double pt() const final
transverse momentum
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
Definition: weight.py:1
void setEvent(const edm::Event &evt)
Load the vertices from the event.
bool isNonnull() const
Checks for non-null.
Definition: RefToBase.h:340
std::vector< CandidatePtr > pfChargedCands(const Jet &jet, bool sort=true)
Extract all non-neutral candidates from a PFJet.
InputIterator leadCand(InputIterator begin, InputIterator end)
std::vector< reco::VertexRef > selectedVertices_
const CandidatePtr getLeadCand(const Jet &) const
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:690
edm::RefToBase< reco::Track > TrackBaseRef
persistent reference to a Track, using views
Definition: TrackFwd.h:36
Definition: Jet.py:1
void setPV(const reco::VertexRef &vtx) const
Update the primary vertex.
const CandidatePtr & leadChargedHadrCand() const
Definition: PFTau.cc:67
double pt() const
track transverse momentum
Definition: TrackBase.h:660
const JetBaseRef & jetRef() const
Definition: PFTau.cc:58
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
Definition: TrackBase.h:814
bool isNull() const
Checks for null.
Definition: Ptr.h:164
const double infinity
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
virtual const reco::Track & pseudoTrack() const
T min(T a, T b)
Definition: MathUtil.h:58
StringCutObjectSelector< reco::Vertex > * vertexSelector_
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:168
RecoTauVertexAssociator(const edm::ParameterSet &pset, edm::ConsumesCollector &&iC)
bool hasTrackDetails() const
Return true if a bestTrack can be extracted from this Candidate.
std::map< const reco::Jet *, reco::VertexRef > * jetToVertexAssociation_
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:648
const PFTauTagInfoRef & pfTauTagInfoRef() const
Definition: PFTau.cc:61
static const TrackGhostTrackState * getTrack(const BasicGhostTrackState *basic)
edm::Ptr< Candidate > CandidatePtr
persistent reference to an object in a collection of Candidate objects
Definition: CandidateFwd.h:25
ParameterSet const & getParameterSet(std::string const &) const
edm::EDGetTokenT< reco::VertexCollection > vxToken_
bool isNull() const
Checks for null.
Definition: RefToBase.h:331
edm::EventID id() const
Definition: EventBase.h:59
fixed size matrix
reco::VertexRef associatedVertex(const Jet &jet) const
double phi() const final
momentum azimuthal angle