CMS 3D CMS Logo

RecoTauQualityCuts.cc
Go to the documentation of this file.
6 
7 #include <boost/bind.hpp>
8 
9 namespace reco { namespace tau {
10 
11 namespace {
12  // Get the KF track if it exists. Otherwise, see if PFCandidate has a GSF track.
13  const reco::TrackBaseRef getTrackRef(const PFCandidate& cand)
14  {
15  if ( cand.trackRef().isNonnull() ) return reco::TrackBaseRef(cand.trackRef());
16  else if ( cand.gsfTrackRef().isNonnull() ) return reco::TrackBaseRef(cand.gsfTrackRef());
17  else return reco::TrackBaseRef();
18  }
19 
20  // Translate GsfTrackRef to TrackBaseRef
21  template <typename T>
22  reco::TrackBaseRef convertRef(const T& ref) {
23  return reco::TrackBaseRef(ref);
24  }
25 }
26 
27 // Quality cut implementations
28 namespace qcuts {
29 
30 bool ptMin(const TrackBaseRef& track, double cut)
31 {
32  LogDebug("TauQCuts") << "<ptMin>: Pt = " << track->pt() << ", cut = " << cut ;
33  return (track->pt() > cut);
34 }
35 
36 bool ptMin_cand(const PFCandidate& cand, double cut)
37 {
38  LogDebug("TauQCuts") << "<ptMin_cand>: Pt = " << cand.pt() << ", cut = " << cut ;
39  return (cand.pt() > cut);
40 }
41 
42 bool etMin_cand(const PFCandidate& cand, double cut)
43 {
44  LogDebug("TauQCuts") << "<etMin_cand>: Et = " << cand.et() << ", cut = " << cut ;
45  return (cand.et() > cut);
46 }
47 
48 bool trkPixelHits(const TrackBaseRef& track, int cut)
49 {
50  // For some reason, the number of hits is signed
51  LogDebug("TauQCuts") << "<trkPixelHits>: #Pxl hits = " << track->hitPattern().numberOfValidPixelHits() << ", cut = " << cut ;
52  return (track->hitPattern().numberOfValidPixelHits() >= cut);
53 }
54 
55 bool trkPixelHits_cand(const PFCandidate& cand, int cut)
56 {
57  // For some reason, the number of hits is signed
58  auto track = getTrackRef(cand);
59  if ( track.isNonnull() ) {
60  LogDebug("TauQCuts") << "<trkPixelHits_cand>: #Pxl hits = " << trkPixelHits(track, cut) << ", cut = " << cut ;
61  return trkPixelHits(track, cut);
62  } else {
63  LogDebug("TauQCuts") << "<trkPixelHits_cand>: #Pxl hits = N/A, cut = " << cut ;
64  return false;
65  }
66 }
67 
69 {
70  LogDebug("TauQCuts") << "<trkTrackerHits>: #Trk hits = " << track->hitPattern().numberOfValidHits() << ", cut = " << cut ;
71  return (track->hitPattern().numberOfValidHits() >= cut);
72 }
73 
74 bool trkTrackerHits_cand(const PFCandidate& cand, int cut)
75 {
76  auto track = getTrackRef(cand);
77  if ( track.isNonnull() ) {
78  LogDebug("TauQCuts") << "<trkTrackerHits>: #Trk hits = " << track->hitPattern().numberOfValidHits() << ", cut = " << cut ;
79  return trkTrackerHits(track, cut);
80  } else {
81  LogDebug("TauQCuts") << "<trkTrackerHits>: #Trk hits = N/A, cut = " << cut ;
82  return false;
83  }
84 }
85 
87 {
88  if ( pv->isNull() ) {
89  edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
90  "RecoTauQualityCuts is invalid. - trkTransverseImpactParameter";
91  return false;
92  }
93  LogDebug("TauQCuts") << " track: Pt = " << track->pt() << ", eta = " << track->eta() << ", phi = " << track->phi() ;
94  LogDebug("TauQCuts") << " vertex: x = " << (*pv)->position().x() << ", y = " << (*pv)->position().y() << ", z = " << (*pv)->position().z() ;
95  LogDebug("TauQCuts") << "--> dxy = " << std::fabs(track->dxy((*pv)->position())) << " (cut = " << cut << ")" ;
96  return (std::fabs(track->dxy((*pv)->position())) <= cut);
97 }
98 
100 {
101  auto track = getTrackRef(cand);
102  if ( track.isNonnull() ) {
103  return trkTransverseImpactParameter(track, pv, cut);
104  } else {
105  LogDebug("TauQCuts") << "<trkTransverseImpactParameter_cand>: dXY = N/A, cut = " << cut ;
106  return false;
107  }
108 }
109 
111 {
112  if ( pv->isNull() ) {
113  edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
114  "RecoTauQualityCuts is invalid. - trkLongitudinalImpactParameter";
115  return false;
116  }
117  LogDebug("TauQCuts") << " track: Pt = " << track->pt() << ", eta = " << track->eta() << ", phi = " << track->phi() ;
118  LogDebug("TauQCuts") << " vertex: x = " << (*pv)->position().x() << ", y = " << (*pv)->position().y() << ", z = " << (*pv)->position().z() ;
119  LogDebug("TauQCuts") << "--> dz = " << std::fabs(track->dz((*pv)->position())) << " (cut = " << cut << ")" ;
120  return (std::fabs(track->dz((*pv)->position())) <= cut);
121 }
122 
124 {
125  auto track = getTrackRef(cand);
126  if ( track.isNonnull() ) {
127  return trkLongitudinalImpactParameter(track, pv, cut);
128  } else {
129  LogDebug("TauQCuts") << "<trkLongitudinalImpactParameter_cand>: dZ = N/A, cut = " << cut ;
130  return false;
131  }
132 }
133 
136 {
137  if ( leadTrack->isNull()) {
138  edm::LogError("QCutsNoValidLeadTrack") << "Lead track Ref in " <<
139  "RecoTauQualityCuts is invalid. - trkLongitudinalImpactParameterWrtTrack";
140  return false;
141  }
142  return (std::fabs(track->dz((*pv)->position()) - (*leadTrack)->dz((*pv)->position())) <= cut);
143 }
144 
146 {
147  auto track = getTrackRef(cand);
148  if ( track.isNonnull() ) return trkLongitudinalImpactParameterWrtTrack(track, leadTrack, pv, cut);
149  else return false;
150 }
151 
153 {
154  if ( pv->isNull() ) {
155  edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
156  "RecoTauQualityCuts is invalid. - minTrackVertexWeight";
157  return false;
158  }
159  LogDebug("TauQCuts") << " track: Pt = " << track->pt() << ", eta = " << track->eta() << ", phi = " << track->phi() ;
160  LogDebug("TauQCuts") << " vertex: x = " << (*pv)->position().x() << ", y = " << (*pv)->position().y() << ", z = " << (*pv)->position().z() ;
161  LogDebug("TauQCuts") << "--> trackWeight = " << (*pv)->trackWeight(track) << " (cut = " << cut << ")" ;
162  return ((*pv)->trackWeight(track) >= cut);
163 }
164 
165 bool minTrackVertexWeight_cand(const PFCandidate& cand, const reco::VertexRef* pv, double cut)
166 {
167  auto track = getTrackRef(cand);
168  if ( track.isNonnull() ) {
169  return minTrackVertexWeight(track, pv, cut);
170  } else {
171  LogDebug("TauQCuts") << "<minTrackVertexWeight_cand>: weight = N/A, cut = " << cut ;
172  return false;
173  }
174 }
175 
176 bool trkChi2(const TrackBaseRef& track, double cut)
177 {
178  LogDebug("TauQCuts") << "<trkChi2>: chi^2 = " << track->normalizedChi2() << ", cut = " << cut ;
179  return (track->normalizedChi2() <= cut);
180 }
181 
182 bool trkChi2_cand(const PFCandidate& cand, double cut)
183 {
184  auto track = getTrackRef(cand);
185  if ( track.isNonnull() ) {
186  LogDebug("TauQCuts") << "<trkChi2_cand>: chi^2 = " << track->normalizedChi2() << ", cut = " << cut ;
187  return trkChi2(track, cut);
188  } else {
189  LogDebug("TauQCuts") << "<trkChi2_cand>: chi^2 = N/A, cut = " << cut ;
190  return false;
191  }
192 }
193 
194 // And a set of qcuts
196 {
197  BOOST_FOREACH( const RecoTauQualityCuts::TrackQCutFunc& func, cuts ) {
198  if ( !func(track) ) return false;
199  }
200  return true;
201 }
202 
204 {
205  BOOST_FOREACH( const RecoTauQualityCuts::CandQCutFunc& func, cuts ) {
206  if ( !func(cand) ) return false;
207  }
208  return true;
209 }
210 
211 // Get the set of Q cuts for a given type (i.e. gamma)
213 {
214  // Find the cuts that for this particle type
215  RecoTauQualityCuts::CandQCutFuncMap::const_iterator cuts = funcMap.find(cand.particleId());
216  // Return false if we dont' know how to deal with this particle type
217  if ( cuts == funcMap.end() ) return false;
218  return AND_cand(cand, cuts->second); // Otherwise AND all the cuts
219 }
220 
221 } // end qcuts implementation namespace
222 
224 {
225  // Setup all of our predicates
226  CandQCutFuncCollection chargedHadronCuts;
227  CandQCutFuncCollection gammaCuts;
228  CandQCutFuncCollection neutralHadronCuts;
229 
230  // Make sure there are no extra passed options
231  std::set<std::string> passedOptionSet;
232  std::vector<std::string> passedOptions = qcuts.getParameterNames();
233 
234  BOOST_FOREACH(const std::string& option, passedOptions) {
235  passedOptionSet.insert(option);
236  }
237 
238  unsigned int nCuts = 0;
239  auto getDouble = [&qcuts, &passedOptionSet, &nCuts](const std::string& name) {
240  if(qcuts.exists(name)) {
241  ++nCuts;
242  passedOptionSet.erase(name);
243  return qcuts.getParameter<double>(name);
244  }
245  return -1.0;
246  };
247  auto getUint = [&qcuts, &passedOptionSet, &nCuts](const std::string& name) -> unsigned int {
248  if(qcuts.exists(name)) {
249  ++nCuts;
250  passedOptionSet.erase(name);
251  return qcuts.getParameter<unsigned int>(name);
252  }
253  return 0;
254  };
255 
256  // Build all the QCuts for tracks
257  minTrackPt_ = getDouble("minTrackPt");
258  maxTrackChi2_ = getDouble("maxTrackChi2");
259  minTrackPixelHits_ = getUint("minTrackPixelHits");
260  minTrackHits_ = getUint("minTrackHits");
261  maxTransverseImpactParameter_ = getDouble("maxTransverseImpactParameter");
262  maxDeltaZ_ = getDouble("maxDeltaZ");
263  maxDeltaZToLeadTrack_ = getDouble("maxDeltaZToLeadTrack");
264  // Require tracks to contribute a minimum weight to the associated vertex.
265  minTrackVertexWeight_ = getDouble("minTrackVertexWeight");
266 
267  // Use bit-wise & to avoid conditional code
268  checkHitPattern_ = (minTrackHits_ > 0) || (minTrackPixelHits_ > 0);
269  checkPV_ = (maxTransverseImpactParameter_ >= 0) ||
270  (maxDeltaZ_ >= 0) ||
271  (maxDeltaZToLeadTrack_ >= 0) ||
272  (minTrackVertexWeight_ >= 0);
273 
274  // Build the QCuts for gammas
275  minGammaEt_ = getDouble("minGammaEt");
276 
277  // Build QCuts for netural hadrons
278  minNeutralHadronEt_ = getDouble("minNeutralHadronEt");
279 
280  // Check if there are any remaining unparsed QCuts
281  if ( !passedOptionSet.empty() ) {
282  std::string unParsedOptions;
283  bool thereIsABadParameter = false;
284  BOOST_FOREACH( const std::string& option, passedOptionSet ) {
285  // Workaround for HLT - TODO FIXME
286  if ( option == "useTracksInsteadOfPFHadrons" ) {
287  // Crash if true - no one should have this option enabled.
288  if ( qcuts.getParameter<bool>("useTracksInsteadOfPFHadrons") ) {
289  throw cms::Exception("DontUseTracksInQcuts")
290  << "The obsolete exception useTracksInsteadOfPFHadrons "
291  << "is set to true in the quality cut config." << std::endl;
292  }
293  continue;
294  }
295 
296  // If we get to this point, there is a real unknown parameter
297  thereIsABadParameter = true;
298 
299  unParsedOptions += option;
300  unParsedOptions += "\n";
301  }
302  if ( thereIsABadParameter ) {
303  throw cms::Exception("BadQualityCutConfig")
304  << " The PSet passed to the RecoTauQualityCuts class had"
305  << " the following unrecognized options: " << std::endl
306  << unParsedOptions;
307  }
308  }
309 
310  // Make sure there are at least some quality cuts
311  if ( !nCuts ) {
312  throw cms::Exception("BadQualityCutConfig")
313  << " No options were passed to the quality cut class!" << std::endl;
314  }
315 }
316 
317 std::pair<edm::ParameterSet, edm::ParameterSet> factorizePUQCuts(const edm::ParameterSet& input)
318 {
319  edm::ParameterSet puCuts;
320  edm::ParameterSet nonPUCuts;
321 
322  std::vector<std::string> inputNames = input.getParameterNames();
323  BOOST_FOREACH( const std::string& cut, inputNames ) {
324  if ( cut == "minTrackVertexWeight" ||
325  cut == "maxDeltaZ" ||
326  cut == "maxDeltaZToLeadTrack" ) {
327  puCuts.copyFrom(input, cut);
328  } else {
329  nonPUCuts.copyFrom(input, cut);
330  }
331  }
332  return std::make_pair(puCuts, nonPUCuts);
333 }
334 
336 {
337  return filterTrack_(track);
338 }
339 
341 {
342  return filterTrack_(track);
343 }
344 
345 template <typename T>
346 bool RecoTauQualityCuts::filterTrack_(const T& trackRef) const
347 {
348  const Track *track = trackRef.get();
349  if(minTrackPt_ >= 0 && !(track->pt() > minTrackPt_)) return false;
350  if(maxTrackChi2_ >= 0 && !(track->normalizedChi2() <= maxTrackChi2_)) return false;
351  if(checkHitPattern_) {
352  const reco::HitPattern &hitPattern = track->hitPattern();
353  if(minTrackPixelHits_ > 0 && !(hitPattern.numberOfValidPixelHits() >= minTrackPixelHits_)) return false;
354  if(minTrackHits_ > 0 && !(hitPattern.numberOfValidHits() >= minTrackHits_)) return false;
355  }
356  if(checkPV_ && pv_.isNull()) {
357  edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
358  "RecoTauQualityCuts is invalid. - filterTrack";
359  return false;
360  }
361 
362  if(maxTransverseImpactParameter_ >= 0 &&
363  !(std::fabs(track->dxy(pv_->position())) <= maxTransverseImpactParameter_))
364  return false;
365  if(maxDeltaZ_ >= 0 && !(std::fabs(track->dz(pv_->position())) <= maxDeltaZ_)) return false;
366  if(maxDeltaZToLeadTrack_ >= 0) {
367  if ( leadTrack_.isNull()) {
368  edm::LogError("QCutsNoValidLeadTrack") << "Lead track Ref in " <<
369  "RecoTauQualityCuts is invalid. - filterTrack";
370  return false;
371  }
372 
373  if(!(std::fabs(track->dz(pv_->position()) - leadTrack_->dz(pv_->position())) <= maxDeltaZToLeadTrack_))
374  return false;
375  }
376  if(minTrackVertexWeight_ > -1.0 && !(pv_->trackWeight(convertRef(trackRef)) >= minTrackVertexWeight_)) return false;
377 
378  return true;
379 }
380 
382  if(minGammaEt_ >= 0 && !(cand.et() > minGammaEt_)) return false;
383  return true;
384 }
385 
387  if(minNeutralHadronEt_ >= 0 && !(cand.et() > minNeutralHadronEt_)) return false;
388  return true;
389 }
390 
392  switch(cand.particleId()) {
393  case PFCandidate::gamma:
394  return filterGammaCand(cand);
395  case PFCandidate::h0:
396  return filterNeutralHadronCand(cand);
397  // We use the same qcuts for muons/electrons and charged hadrons.
398  case PFCandidate::h:
399  case PFCandidate::e:
400  case PFCandidate::mu:
401  // no cuts ATM (track cuts applied in filterCand)
402  return true;
403  // Return false if we dont' know how to deal with this particle type
404  default:
405  return false;
406  };
407  return false;
408 }
409 
411 {
412  auto trackRef = cand.trackRef();
413  bool result = true;
414  if(trackRef.isNonnull()) {
415  result = filterTrack_(trackRef);
416  }
417  else {
418  auto gsfTrackRef = cand.gsfTrackRef();
419  if(gsfTrackRef.isNonnull()) {
420  result = filterTrack_(gsfTrackRef);
421  }
422  }
423  if(result)
424  result = filterCandByType(cand);
425  return result;
426 }
427 
429 {
430  leadTrack_ = reco::TrackBaseRef(leadTrack);
431 }
432 
434 {
435  leadTrack_ = getTrackRef(leadCand);
436 }
437 
439 {
440  if ( leadCand.isNonnull() ) {
441  leadTrack_ = getTrackRef(*leadCand);
442  } else {
443  // Set null
444  leadTrack_ = reco::TrackBaseRef();
445  }
446 }
447 
448 }} // end namespace reco::tau
#define LogDebug(id)
T getParameter(std::string const &) const
bool filterCandByType(const reco::PFCandidate &cand) const
bool trkLongitudinalImpactParameter_cand(const PFCandidate &cand, const reco::VertexRef *pv, double cut)
bool minTrackVertexWeight(const TrackBaseRef &track, const reco::VertexRef *pv, double cut)
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:253
bool trkTransverseImpactParameter_cand(const PFCandidate &cand, const reco::VertexRef *pv, double cut)
RecoTauQualityCuts(const edm::ParameterSet &qcuts)
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:561
bool minTrackVertexWeight_cand(const PFCandidate &cand, const reco::VertexRef *pv, double cut)
int numberOfValidHits() const
Definition: HitPattern.h:823
bool exists(std::string const &parameterName) const
checks if a parameter exists
bool etMin_cand(const PFCandidate &cand, double cut)
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:645
double pt() const final
transverse momentum
bool trkPixelHits(const TrackBaseRef &track, int cut)
bool filterTrack_(const T &trackRef) const
bool filterNeutralHadronCand(const reco::PFCandidate &cand) const
bool trkTransverseImpactParameter(const TrackBaseRef &track, const reco::VertexRef *pv, double cut)
static std::string const input
Definition: EdmProvDump.cc:44
std::vector< CandQCutFunc > CandQCutFuncCollection
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:438
bool filterGammaCand(const reco::PFCandidate &cand) const
void copyFrom(ParameterSet const &from, std::string const &name)
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:651
bool trkLongitudinalImpactParameterWrtTrack_cand(const PFCandidate &cand, const reco::TrackBaseRef *leadTrack, const reco::VertexRef *pv, double cut)
edm::RefToBase< reco::Track > TrackBaseRef
persistent reference to a Track, using views
Definition: TrackFwd.h:32
double et() const final
transverse energy
std::map< PFCandidate::ParticleType, CandQCutFuncCollection > CandQCutFuncMap
double pt() const
track transverse momentum
Definition: TrackBase.h:621
def pv(vc)
Definition: MetAnalyzer.py:6
std::pair< edm::ParameterSet, edm::ParameterSet > factorizePUQCuts(const edm::ParameterSet &inputSet)
std::vector< TrackQCutFunc > TrackQCutFuncCollection
std::vector< std::string > getParameterNames() const
bool ptMin(const TrackBaseRef &track, double cut)
bool isNull() const
Checks for null.
Definition: Ref.h:250
bool trkTrackerHits(const TrackBaseRef &track, int cut)
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:609
bool trkTrackerHits_cand(const PFCandidate &cand, int cut)
bool mapAndCutByType(const PFCandidate &cand, const RecoTauQualityCuts::CandQCutFuncMap &funcMap)
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:446
boost::function< bool(const TrackBaseRef &)> TrackQCutFunc
bool ptMin_cand(const PFCandidate &cand, double cut)
bool trkChi2_cand(const PFCandidate &cand, double cut)
boost::function< bool(const PFCandidate &)> CandQCutFunc
bool trkChi2(const TrackBaseRef &track, double cut)
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:39
fixed size matrix
bool AND(const TrackBaseRef &track, const RecoTauQualityCuts::TrackQCutFuncCollection &cuts)
bool trkLongitudinalImpactParameter(const TrackBaseRef &track, const reco::VertexRef *pv, double cut)
reco::GsfTrackRef gsfTrackRef() const
Definition: PFCandidate.cc:476
bool trkPixelHits_cand(const PFCandidate &cand, int cut)
int numberOfValidPixelHits() const
Definition: HitPattern.h:838
void setLeadTrack(const reco::TrackRef &leadTrack) const
Update the leading track.
bool trkLongitudinalImpactParameterWrtTrack(const TrackBaseRef &track, const reco::TrackBaseRef *leadTrack, const reco::VertexRef *pv, double cut)
DZ cut, with respect to the current lead rack.
bool filterCand(const reco::PFCandidate &cand) const
Filter a single PFCandidate.
virtual ParticleType particleId() const
Definition: PFCandidate.h:373
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:591
long double T
bool AND_cand(const PFCandidate &cand, const RecoTauQualityCuts::CandQCutFuncCollection &cuts)
bool filterTrack(const reco::TrackBaseRef &track) const
Filter a single Track.