CMS 3D CMS Logo

PFRecoTauDiscriminationByIsolationContainer.cc
Go to the documentation of this file.
1 #include <functional>
2 #include <memory>
3 
14 
15 #include "TMath.h"
16 #include "TFormula.h"
17 
18 /* class PFRecoTauDiscriminationByIsolationContainer
19  * created : Jul 23 2007,
20  * revised : Thu Aug 13 14:44:40 PDT 2009
21  * contributors : Ludovic Houchu (IPHC, Strasbourg),
22  * Christian Veelken (UC Davis),
23  * Evan K. Friis (UC Davis)
24  * Michalis Bachtis (UW Madison)
25  */
26 
27 using namespace reco;
28 using namespace std;
29 
31 public:
32  enum StoredRawType { None, SumPt, PUsumPt, Occupancy, FootPrintCorrection, PhotonSumPt };
35  moduleLabel_(pset.getParameter<std::string>("@module_label")),
36  qualityCutsPSet_(pset.getParameter<edm::ParameterSet>("qualityCuts")) {
37  // RIC: multiply neutral isolation by a flat factor.
38  // Useful, for instance, to combine charged and neutral isolations
39  // with different relative weights
40  weightGammas_ = pset.getParameter<double>("WeightECALIsolation");
41 
42  // RIC: allow to relax the isolation completely beyond a given tau pt
43  minPtForNoIso_ = pset.getParameter<double>("minTauPtForNoIso");
44 
45  // Get configs for raw values
46  bool storeRawFootprintCorrection = false;
47  deltaBetaNeeded_ = false;
48  weightsNeeded_ = false;
49  tracksNeeded_ = false;
50  gammasNeeded_ = false;
51  storeRawValue_.clear();
52  auto const& rawDefs = pset.getParameter<std::vector<edm::ParameterSet>>("IDdefinitions");
53  std::vector<std::string> idnames;
54  for (auto const& rawDefsEntry : rawDefs) {
55  idnames.push_back(rawDefsEntry.getParameter<std::string>("IDname"));
56  // Can only store one type
57  int numStoreOptions = 0;
58  if (rawDefsEntry.getParameter<bool>("storeRawSumPt")) {
59  storeRawValue_.push_back(SumPt);
60  ++numStoreOptions;
61  }
62  if (rawDefsEntry.getParameter<bool>("storeRawOccupancy")) {
63  storeRawValue_.push_back(Occupancy);
64  ++numStoreOptions;
65  }
66  if (rawDefsEntry.getParameter<bool>("storeRawPUsumPt")) {
67  storeRawValue_.push_back(PUsumPt);
68  ++numStoreOptions;
69  }
70  if (rawDefsEntry.getParameter<bool>("storeRawFootprintCorrection")) {
71  storeRawValue_.push_back(FootPrintCorrection);
73  ++numStoreOptions;
74  }
75  if (rawDefsEntry.getParameter<bool>("storeRawPhotonSumPt_outsideSignalCone")) {
76  storeRawValue_.push_back(PhotonSumPt);
77  ++numStoreOptions;
78  }
79  if (numStoreOptions != 1) {
80  throw cms::Exception("BadIsoConfig")
81  << "Multiple or none of 'store sum pt' and/or 'store occupancy' options are set."
82  << " These options are mutually exclusive.";
83  }
84 
85  includeGammas_.push_back(rawDefsEntry.getParameter<bool>("ApplyDiscriminationByECALIsolation"));
86  if (includeGammas_.back())
87  gammasNeeded_ = true;
88  calculateWeights_.push_back(rawDefsEntry.getParameter<bool>("ApplyDiscriminationByWeightedECALIsolation"));
89  if (calculateWeights_.back())
90  weightsNeeded_ = true;
91  includeTracks_.push_back(rawDefsEntry.getParameter<bool>("ApplyDiscriminationByTrackerIsolation"));
92  if (includeTracks_.back())
93  tracksNeeded_ = true;
94  applyDeltaBetaCorrection_.push_back(rawDefsEntry.getParameter<bool>("applyDeltaBetaCorrection"));
95  if (applyDeltaBetaCorrection_.back())
96  deltaBetaNeeded_ = true;
97  useAllPFCandsForWeights_.push_back(rawDefsEntry.getParameter<bool>("UseAllPFCandsForWeights"));
98 
99  // sanity check2 - can't use weighted and unweighted iso at the same time
100  if (includeGammas_.back() && calculateWeights_.back()) {
101  throw cms::Exception("BadIsoConfig")
102  << "Both 'ApplyDiscriminationByECALIsolation' and 'ApplyDiscriminationByWeightedECALIsolation' "
103  << "have been set to true. These options are mutually exclusive.";
104  }
105  }
106 
107  // Get configs for WPs - negative cut values are used to switch of the condition
108  std::vector<edm::ParameterSet> wpDefs = pset.getParameter<std::vector<edm::ParameterSet>>("IDWPdefinitions");
109  for (std::vector<edm::ParameterSet>::iterator wpDefsEntry = wpDefs.begin(); wpDefsEntry != wpDefs.end();
110  ++wpDefsEntry) {
111  maxAbsValue_.push_back(wpDefsEntry->getParameter<std::vector<double>>("maximumAbsoluteValues"));
112  maxRelValue_.push_back(wpDefsEntry->getParameter<std::vector<double>>("maximumRelativeValues"));
113  offsetRelValue_.push_back(wpDefsEntry->getParameter<std::vector<double>>("relativeValueOffsets"));
114  auto refRawIDNames = wpDefsEntry->getParameter<std::vector<std::string>>("referenceRawIDNames");
115  if (!maxAbsValue_.back().empty() && maxAbsValue_.back().size() != refRawIDNames.size())
116  throw cms::Exception("BadIsoConfig")
117  << "WP configuration: Length of 'maximumAbsoluteValues' does not match length of 'referenceRawIDNames'!";
118  if (!maxRelValue_.back().empty() && maxRelValue_.back().size() != refRawIDNames.size())
119  throw cms::Exception("BadIsoConfig")
120  << "WP configuration: Length of 'maximumRelativeValues' does not match length of 'referenceRawIDNames'!";
121  if (!offsetRelValue_.back().empty() && offsetRelValue_.back().size() != refRawIDNames.size())
122  throw cms::Exception("BadIsoConfig")
123  << "WP configuration: Length of 'relativeValueOffsets' does not match length of 'referenceRawIDNames'!";
124  else if (offsetRelValue_.back().empty())
125  offsetRelValue_.back().assign(refRawIDNames.size(), 0.0);
126  rawValue_reference_.push_back(std::vector<int>(refRawIDNames.size()));
127  for (size_t i = 0; i < refRawIDNames.size(); i++) {
128  bool found = false;
129  for (size_t j = 0; j < idnames.size(); j++) {
130  if (refRawIDNames[i] == idnames[j]) {
131  rawValue_reference_.back()[i] = j;
132  found = true;
133  break;
134  }
135  }
136  if (!found)
137  throw cms::Exception("BadIsoConfig")
138  << "WP configuration: Requested raw ID '" << refRawIDNames[i] << "' not defined!";
139  }
140  }
141 
142  customIsoCone_ = pset.getParameter<double>("customOuterCone");
143 
144  applyFootprintCorrection_ = pset.getParameter<bool>("applyFootprintCorrection");
145  if (applyFootprintCorrection_ || storeRawFootprintCorrection) {
146  edm::VParameterSet cfgFootprintCorrections = pset.getParameter<edm::VParameterSet>("footprintCorrections");
147  for (edm::VParameterSet::const_iterator cfgFootprintCorrection = cfgFootprintCorrections.begin();
148  cfgFootprintCorrection != cfgFootprintCorrections.end();
149  ++cfgFootprintCorrection) {
150  std::string selection = cfgFootprintCorrection->getParameter<std::string>("selection");
151  std::string offset = cfgFootprintCorrection->getParameter<std::string>("offset");
152  std::unique_ptr<FootprintCorrection> footprintCorrection(new FootprintCorrection(selection, offset));
153  footprintCorrections_.push_back(std::move(footprintCorrection));
154  }
155  }
156 
157  // Get the quality cuts specific to the isolation region
158  edm::ParameterSet isolationQCuts = qualityCutsPSet_.getParameterSet("isolationQualityCuts");
159 
160  qcuts_ = std::make_unique<tau::RecoTauQualityCuts>(isolationQCuts);
161 
162  vertexAssociator_ = std::make_unique<tau::RecoTauVertexAssociator>(qualityCutsPSet_, consumesCollector());
163 
164  if (deltaBetaNeeded_ || weightsNeeded_) {
165  // Factorize the isolation QCuts into those that are used to
166  // select PU and those that are not.
167  std::pair<edm::ParameterSet, edm::ParameterSet> puFactorizedIsoQCuts =
168  reco::tau::factorizePUQCuts(isolationQCuts);
169 
170  // Determine the pt threshold for the PU tracks
171  // First check if the user specifies explicitly the cut.
172  // For that the user has to provide a >= 0 value for the PtCutOverride.
173  bool deltaBetaPUTrackPtCutOverride = pset.getParameter<bool>("deltaBetaPUTrackPtCutOverride");
175  double deltaBetaPUTrackPtCutOverride_val = pset.getParameter<double>("deltaBetaPUTrackPtCutOverride_val");
176  puFactorizedIsoQCuts.second.addParameter<double>("minTrackPt", deltaBetaPUTrackPtCutOverride_val);
177  } else {
178  // Secondly take it from the minGammaEt
179  puFactorizedIsoQCuts.second.addParameter<double>("minTrackPt",
180  isolationQCuts.getParameter<double>("minGammaEt"));
181  }
182 
183  pileupQcutsPUTrackSelection_ = std::make_unique<tau::RecoTauQualityCuts>(puFactorizedIsoQCuts.first);
184 
185  pileupQcutsGeneralQCuts_ = std::make_unique<tau::RecoTauQualityCuts>(puFactorizedIsoQCuts.second);
186 
187  pfCandSrc_ = pset.getParameter<edm::InputTag>("particleFlowSrc");
188  pfCand_token = consumes<edm::View<reco::Candidate>>(pfCandSrc_);
189  vertexSrc_ = pset.getParameter<edm::InputTag>("vertexSrc");
190  vertex_token = consumes<reco::VertexCollection>(vertexSrc_);
191  deltaBetaCollectionCone_ = pset.getParameter<double>("isoConeSizeForDeltaBeta");
192  std::string deltaBetaFactorFormula = pset.getParameter<string>("deltaBetaFactor");
193  deltaBetaFormula_ = std::make_unique<TFormula>("DB_corr", deltaBetaFactorFormula.c_str());
194  }
195 
196  applyRhoCorrection_ = pset.getParameter<bool>("applyRhoCorrection");
197  if (applyRhoCorrection_) {
198  rhoProducer_ = pset.getParameter<edm::InputTag>("rhoProducer");
199  rho_token = consumes<double>(rhoProducer_);
200  rhoConeSize_ = pset.getParameter<double>("rhoConeSize");
201  rhoUEOffsetCorrection_ = pset.getParameter<double>("rhoUEOffsetCorrection");
202  }
203 
204  verbosity_ = pset.getParameter<int>("verbosity");
205  }
206 
208 
209  void beginEvent(const edm::Event& evt, const edm::EventSetup& evtSetup) override;
210  reco::SingleTauDiscriminatorContainer discriminate(const PFTauRef& pfTau) const override;
211 
212  inline double weightedSum(const std::vector<CandidatePtr>& inColl_, double eta, double phi) const {
213  double out = 1.0;
214  for (auto const& inObj_ : inColl_) {
215  double sum = (inObj_->pt() * inObj_->pt()) / (deltaR2(eta, phi, inObj_->eta(), inObj_->phi()));
216  if (sum > 1.0)
217  out *= sum;
218  }
219  return out;
220  }
221 
222  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
223 
224 private:
226 
228  std::unique_ptr<tau::RecoTauQualityCuts> qcuts_;
229 
230  // Inverted QCut which selects tracks with bad DZ/trackWeight
231  std::unique_ptr<tau::RecoTauQualityCuts> pileupQcutsPUTrackSelection_;
232  std::unique_ptr<tau::RecoTauQualityCuts> pileupQcutsGeneralQCuts_;
233 
234  std::unique_ptr<tau::RecoTauVertexAssociator> vertexAssociator_;
235 
241  // RIC:
243 
247  : selection_(selection), offset_(offset) {}
251  };
252  std::vector<std::unique_ptr<FootprintCorrection>> footprintCorrections_;
253 
254  // Options to store the raw value in the discriminator instead of boolean pass/fail flag
255  std::vector<StoredRawType> storeRawValue_;
256  // Options to store the boolean pass/fail flag
257  std::vector<std::vector<int>> rawValue_reference_;
258  std::vector<std::vector<double>> maxAbsValue_;
259  std::vector<std::vector<double>> maxRelValue_;
260  std::vector<std::vector<double>> offsetRelValue_;
261  // Options used for both raw and WP definitions
262  std::vector<bool> includeGammas_;
263  std::vector<bool> calculateWeights_;
264  std::vector<bool> includeTracks_;
265  std::vector<bool> applyDeltaBetaCorrection_;
266  std::vector<bool> useAllPFCandsForWeights_;
267 
268  /* **********************************************************************
269  **** Pileup Subtraction Parameters ***********************************
270  **********************************************************************/
271 
272  // Delta Beta correction
276  // Keep track of how many vertices are in the event
279  std::vector<reco::CandidatePtr> chargedPFCandidatesInEvent_;
280  // Size of cone used to collect PU tracks
282  std::unique_ptr<TFormula> deltaBetaFormula_;
284 
285  // Rho correction
289  double rhoConeSize_;
293 
294  // Flag to enable/disable debug output
296 };
297 
299  const edm::EventSetup& eventSetup) {
300  // NB: The use of the PV in this context is necessitated by its use in
301  // applying quality cuts to the different objects in the isolation cone
302  // The vertex associator contains the logic to select the appropriate vertex
303  // We need to pass it the event so it can load the vertices.
304  vertexAssociator_->setEvent(event);
305 
306  // If we are applying the delta beta correction, we need to get the PF
307  // candidates from the event so we can find the PU tracks.
308  if (deltaBetaNeeded_ || weightsNeeded_) {
309  // Collect all the PF pile up tracks
311  event.getByToken(pfCand_token, pfCandidates);
312  chargedPFCandidatesInEvent_.clear();
313  chargedPFCandidatesInEvent_.reserve(pfCandidates->size());
314  size_t numPFCandidates = pfCandidates->size();
315  for (size_t i = 0; i < numPFCandidates; ++i) {
316  reco::CandidatePtr pfCandidate(pfCandidates, i);
317  if (pfCandidate->charge() != 0) {
318  chargedPFCandidatesInEvent_.push_back(pfCandidate);
319  }
320  }
321  // Count all the vertices in the event, to parameterize the DB
322  // correction factor
324  event.getByToken(vertex_token, vertices);
325  size_t nVtxThisEvent = vertices->size();
326  deltaBetaFactorThisEvent_ = deltaBetaFormula_->Eval(nVtxThisEvent);
327  }
328 
329  if (applyRhoCorrection_) {
330  edm::Handle<double> rhoHandle_;
331  event.getByToken(rho_token, rhoHandle_);
332  rhoThisEvent_ = (*rhoHandle_ - rhoUEOffsetCorrection_) * (3.14159) * rhoConeSize_ * rhoConeSize_;
333  }
334 }
335 
337  const PFTauRef& pfTau) const {
338  LogDebug("discriminate") << " tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta() << ", phi = " << pfTau->phi();
339  LogDebug("discriminate") << *pfTau;
340 
341  // collect the objects we are working with (ie tracks, tracks+gammas, etc)
342  std::vector<CandidatePtr> isoCharged_;
343  std::vector<CandidatePtr> isoNeutral_;
344  std::vector<CandidatePtr> isoPU_;
345  std::vector<CandidatePtr> isoPUall_;
346  CandidateCollection isoNeutralWeight_;
347  CandidateCollection isoNeutralWeight_UseAllPFCands_;
348  std::vector<CandidatePtr> chPV_;
349  std::vector<CandidatePtr> chPVall_;
350  isoCharged_.reserve(pfTau->isolationChargedHadrCands().size());
351  isoNeutral_.reserve(pfTau->isolationGammaCands().size());
352  isoPU_.reserve(std::min(100UL, chargedPFCandidatesInEvent_.size()));
353  isoPUall_.reserve(std::min(100UL, chargedPFCandidatesInEvent_.size()));
354  isoNeutralWeight_.reserve(pfTau->isolationGammaCands().size());
355  isoNeutralWeight_UseAllPFCands_.reserve(pfTau->isolationGammaCands().size());
356 
357  chPV_.reserve(std::min(50UL, chargedPFCandidatesInEvent_.size()));
358  chPVall_.reserve(std::min(50UL, chargedPFCandidatesInEvent_.size()));
359 
360  // Get the primary vertex associated to this tau
361  reco::VertexRef pv = vertexAssociator_->associatedVertex(*pfTau);
362  // Let the quality cuts know which the vertex to use when applying selections
363  // on dz, etc.
364  if (verbosity_) {
365  if (pv.isNonnull()) {
366  LogTrace("discriminate") << "pv: x = " << pv->position().x() << ", y = " << pv->position().y()
367  << ", z = " << pv->position().z();
368  } else {
369  LogTrace("discriminate") << "pv: N/A";
370  }
371  if (pfTau->leadChargedHadrCand().isNonnull()) {
372  LogTrace("discriminate") << "leadPFChargedHadron:"
373  << " Pt = " << pfTau->leadChargedHadrCand()->pt() << ","
374  << " eta = " << pfTau->leadChargedHadrCand()->eta() << ","
375  << " phi = " << pfTau->leadChargedHadrCand()->phi();
376  } else {
377  LogTrace("discriminate") << "leadPFChargedHadron: N/A";
378  }
379  }
380 
381  // CV: isolation is not well defined in case primary vertex or leading charged hadron do not exist
382  if (!(pv.isNonnull() && pfTau->leadChargedHadrCand().isNonnull()))
383  return 0.;
384 
385  qcuts_->setPV(pv);
386  qcuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
387 
388  if (deltaBetaNeeded_ || weightsNeeded_) {
389  pileupQcutsGeneralQCuts_->setPV(pv);
390  pileupQcutsGeneralQCuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
391  pileupQcutsPUTrackSelection_->setPV(pv);
392  pileupQcutsPUTrackSelection_->setLeadTrack(*pfTau->leadChargedHadrCand());
393  }
394 
395  // Load the tracks if they are being used.
396  if (tracksNeeded_) {
397  for (auto const& cand : pfTau->isolationChargedHadrCands()) {
398  if (qcuts_->filterCandRef(cand)) {
399  LogTrace("discriminate") << "adding charged iso cand with pt " << cand->pt();
400  isoCharged_.push_back(cand);
401  }
402  }
403  }
404  if (gammasNeeded_ || weightsNeeded_) {
405  for (auto const& cand : pfTau->isolationGammaCands()) {
406  if (qcuts_->filterCandRef(cand)) {
407  LogTrace("discriminate") << "adding neutral iso cand with pt " << cand->pt();
408  isoNeutral_.push_back(cand);
409  }
410  }
411  }
412 
415 
416  // If desired, get PU tracks.
417  if (deltaBetaNeeded_ || weightsNeeded_) {
418  // First select by inverted the DZ/track weight cuts. True = invert
419  if (verbosity_) {
420  std::cout << "Initial PFCands: " << chargedPFCandidatesInEvent_.size() << std::endl;
421  }
422 
423  std::vector<CandidatePtr> allPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_, true);
424 
425  std::vector<CandidatePtr> allNPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_);
426  LogTrace("discriminate") << "After track cuts: " << allPU.size();
427 
428  // Now apply the rest of the cuts, like pt, and TIP, tracker hits, etc
429  std::vector<CandidatePtr> cleanPU = pileupQcutsGeneralQCuts_->filterCandRefs(allPU);
430 
431  std::vector<CandidatePtr> cleanNPU = pileupQcutsGeneralQCuts_->filterCandRefs(allNPU);
432 
433  LogTrace("discriminate") << "After cleaning cuts: " << cleanPU.size();
434 
435  // Only select PU tracks inside the isolation cone.
436  DRFilter deltaBetaFilter(pfTau->p4(), 0, deltaBetaCollectionCone_);
437  for (auto const& cand : cleanPU) {
438  if (deltaBetaFilter(cand))
439  isoPU_.push_back(cand);
440  }
441 
442  for (auto const& cand : cleanNPU) {
443  if (deltaBetaFilter(cand))
444  chPV_.push_back(cand);
445  }
446  LogTrace("discriminate") << "After cone cuts: " << isoPU_.size() << " " << chPV_.size();
447  isoPUall_ = std::move(allPU);
448  chPVall_ = std::move(allNPU);
449  }
450 
451  if (weightsNeeded_) {
452  for (auto const& isoObject : isoNeutral_) {
453  if (isoObject->charge() != 0) {
454  // weight only neutral objects
455  isoNeutralWeight_.push_back(*isoObject);
456  isoNeutralWeight_UseAllPFCands_.push_back(*isoObject);
457  continue;
458  }
459 
460  double eta = isoObject->eta();
461  double phi = isoObject->phi();
462  {
463  double sumNPU = 0.5 * log(weightedSum(chPV_, eta, phi));
464 
465  double sumPU = 0.5 * log(weightedSum(isoPU_, eta, phi));
466  LeafCandidate neutral(*isoObject);
467  if ((sumNPU + sumPU) > 0)
468  neutral.setP4(((sumNPU) / (sumNPU + sumPU)) * neutral.p4());
469 
470  isoNeutralWeight_.push_back(neutral);
471  }
472  {
473  double sumNPU = 0.5 * log(weightedSum(chPVall_, eta, phi));
474 
475  double sumPU = 0.5 * log(weightedSum(isoPUall_, eta, phi));
476  LeafCandidate neutral(*isoObject);
477  if ((sumNPU + sumPU) > 0)
478  neutral.setP4(((sumNPU) / (sumNPU + sumPU)) * neutral.p4());
479 
480  isoNeutralWeight_UseAllPFCands_.push_back(neutral);
481  }
482  }
483  }
484 
485  // Check if we want a custom iso cone
486  if (customIsoCone_ >= 0.) {
487  DRFilter filter(pfTau->p4(), 0, customIsoCone_);
488  DRFilter2 filter2(pfTau->p4(), 0, customIsoCone_);
489  std::vector<CandidatePtr> isoCharged_filter;
490  std::vector<CandidatePtr> isoNeutral_filter;
491  // Remove all the objects not in our iso cone
492  for (auto const& isoObject : isoCharged_) {
493  if (filter(isoObject))
494  isoCharged_filter.push_back(isoObject);
495  }
496  isoCharged_ = isoCharged_filter;
497  for (auto const& isoObject : isoNeutral_) {
498  if (filter(isoObject))
499  isoNeutral_filter.push_back(isoObject);
500  }
501  isoNeutral_ = isoNeutral_filter;
502  {
503  CandidateCollection isoNeutralWeight_filter;
504  for (auto const& isoObject : isoNeutralWeight_) {
505  if (filter2(isoObject))
506  isoNeutralWeight_filter.push_back(isoObject);
507  }
508  isoNeutralWeight_ = isoNeutralWeight_filter;
509  }
510  {
511  CandidateCollection isoNeutralWeight_filter;
512  for (auto const& isoObject : isoNeutralWeight_UseAllPFCands_) {
513  if (filter2(isoObject))
514  isoNeutralWeight_filter.push_back(isoObject);
515  }
516  isoNeutralWeight_UseAllPFCands_ = isoNeutralWeight_filter;
517  }
518  }
519 
520  //Now all needed incredients are ready. Loop over all ID configurations and produce output
522  for (size_t i = 0; i < includeGammas_.size(); i++) {
523  //--- nObjects requirement
524  int neutrals = isoNeutral_.size();
525 
526  if (applyDeltaBetaCorrection_.at(i)) {
527  neutrals -= TMath::Nint(deltaBetaFactorThisEvent_ * isoPU_.size());
528  }
529  if (neutrals < 0) {
530  neutrals = 0;
531  }
532 
533  int nOccupants = isoCharged_.size() + neutrals;
534 
535  double footprintCorrection_value = 0.;
536  if (applyFootprintCorrection_ || storeRawValue_.at(i) == FootPrintCorrection) {
537  for (std::vector<std::unique_ptr<FootprintCorrection>>::const_iterator footprintCorrection =
538  footprintCorrections_.begin();
539  footprintCorrection != footprintCorrections_.end();
540  ++footprintCorrection) {
541  if ((*footprintCorrection)->selection_(*pfTau)) {
542  footprintCorrection_value = (*footprintCorrection)->offset_(*pfTau);
543  }
544  }
545  }
546 
547  double totalPt = 0.;
548  double puPt = 0.;
549  //--- Sum PT requirement
550  if (storeRawValue_.at(i) == SumPt || storeRawValue_.at(i) == PUsumPt) {
551  double chargedPt = 0.;
552  double neutralPt = 0.;
553  double weightedNeutralPt = 0.;
554  if (includeTracks_.at(i)) {
555  for (auto const& isoObject : isoCharged_) {
556  chargedPt += isoObject->pt();
557  }
558  }
559 
560  if (calculateWeights_.at(i)) {
561  if (useAllPFCandsForWeights_.at(i)) {
562  for (auto const& isoObject : isoNeutralWeight_UseAllPFCands_) {
563  weightedNeutralPt += isoObject.pt();
564  }
565  } else {
566  for (auto const& isoObject : isoNeutralWeight_) {
567  weightedNeutralPt += isoObject.pt();
568  }
569  }
570  } else if (includeGammas_.at(i)) {
571  for (auto const& isoObject : isoNeutral_) {
572  neutralPt += isoObject->pt();
573  }
574  }
575  for (auto const& isoObject : isoPU_) {
576  puPt += isoObject->pt();
577  }
578  LogTrace("discriminate") << "chargedPt = " << chargedPt;
579  LogTrace("discriminate") << "neutralPt = " << neutralPt;
580  LogTrace("discriminate") << "weighted neutral Pt = " << weightedNeutralPt;
581  LogTrace("discriminate") << "puPt = " << puPt << " (delta-beta corr. = " << (deltaBetaFactorThisEvent_ * puPt)
582  << ")";
583 
584  if (calculateWeights_.at(i)) {
585  neutralPt = weightedNeutralPt;
586  }
587 
588  if (applyDeltaBetaCorrection_.at(i)) {
589  neutralPt -= (deltaBetaFactorThisEvent_ * puPt);
590  }
591 
592  if (applyFootprintCorrection_) {
593  neutralPt -= footprintCorrection_value;
594  }
595 
596  if (applyRhoCorrection_) {
597  neutralPt -= rhoThisEvent_;
598  }
599 
600  if (neutralPt < 0.) {
601  neutralPt = 0.;
602  }
603 
604  totalPt = chargedPt + weightGammas_ * neutralPt;
605  }
606 
607  double photonSumPt_outsideSignalCone = 0.;
608  if (storeRawValue_.at(i) == PhotonSumPt) {
609  const std::vector<reco::CandidatePtr>& signalGammas = pfTau->signalGammaCands();
610  for (std::vector<reco::CandidatePtr>::const_iterator signalGamma = signalGammas.begin();
611  signalGamma != signalGammas.end();
612  ++signalGamma) {
613  double dR = deltaR(pfTau->eta(), pfTau->phi(), (*signalGamma)->eta(), (*signalGamma)->phi());
614  if (dR > pfTau->signalConeSize())
615  photonSumPt_outsideSignalCone += (*signalGamma)->pt();
616  }
617  }
618 
619  // We did error checking in the constructor, so this is safe.
620  if (storeRawValue_.at(i) == SumPt) {
621  result.rawValues.push_back(totalPt);
622  } else if (storeRawValue_.at(i) == PUsumPt) {
623  if (applyDeltaBetaCorrection_.at(i))
624  result.rawValues.push_back(puPt);
625  else if (applyRhoCorrection_)
626  result.rawValues.push_back(rhoThisEvent_);
627  else
628  result.rawValues.push_back(0.);
629  } else if (storeRawValue_.at(i) == Occupancy) {
630  result.rawValues.push_back(nOccupants);
631  } else if (storeRawValue_.at(i) == FootPrintCorrection) {
632  result.rawValues.push_back(footprintCorrection_value);
633  } else if (storeRawValue_.at(i) == PhotonSumPt) {
634  result.rawValues.push_back(photonSumPt_outsideSignalCone);
635  }
636  }
637  for (size_t i = 0; i < rawValue_reference_.size(); i++) {
638  bool pass = true;
639  if (minPtForNoIso_ > 0. && pfTau->pt() > minPtForNoIso_)
640  LogDebug("discriminate") << "tau pt = " << pfTau->pt() << "\t min cutoff pt = " << minPtForNoIso_;
641  else {
642  for (size_t j = 0; j < rawValue_reference_[i].size(); j++) {
643  double rawValue = result.rawValues[rawValue_reference_[i][j]];
644  LogTrace("discriminate") << "Iso sum = " << rawValue << " (max_abs = " << maxAbsValue_[i][j]
645  << ", max_rel = " << maxRelValue_[i][j] << ", offset_rel = " << offsetRelValue_[i][j]
646  << ")";
647  if (!maxAbsValue_[i].empty() && maxAbsValue_[i][j] >= 0.0)
648  pass = rawValue <= maxAbsValue_[i][j];
649  if (!maxRelValue_[i].empty() && maxRelValue_[i][j] >= 0.0)
650  pass = rawValue <= maxRelValue_[i][j] * (pfTau->pt() - offsetRelValue_[i][j]);
651  if (!pass)
652  break; // do not pass if one of the conditions in the j list fails
653  }
654  }
655  result.workingPoints.push_back(pass);
656  }
657  return result;
658 }
659 
661  // pfRecoTauDiscriminationByIsolationContainer
663  desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfRecoTauProducer"));
664 
665  edm::ParameterSetDescription desc_qualityCuts;
667  desc.add<edm::ParameterSetDescription>("qualityCuts", desc_qualityCuts);
668 
669  desc.add<double>("minTauPtForNoIso", -99.0);
670  desc.add<edm::InputTag>("vertexSrc", edm::InputTag("offlinePrimaryVertices"));
671  desc.add<double>("rhoConeSize", 0.5);
672  desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAll"));
673 
674  {
676  vpsd1.add<std::string>("selection");
677  vpsd1.add<std::string>("offset");
678  desc.addVPSet("footprintCorrections", vpsd1, {});
679  }
680 
681  desc.add<std::string>("deltaBetaFactor", "0.38");
682  desc.add<bool>("applyFootprintCorrection", false);
683  {
684  edm::ParameterSetDescription pset_Prediscriminants;
685  pset_Prediscriminants.add<std::string>("BooleanOperator", "and");
686  {
688  psd1.add<double>("cut", 0.5);
689  psd1.add<edm::InputTag>("Producer", edm::InputTag("pfRecoTauDiscriminationByLeadingTrackFinding"));
690  pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("leadTrack", psd1);
691  }
692  {
694  psd1.add<double>("cut", 0.5);
695  psd1.add<edm::InputTag>("Producer", edm::InputTag("hpsPFTauDiscriminationByDecayModeFindingNewDMs"));
696  pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("decayMode", psd1);
697  }
698  {
700  psd1.add<double>("cut", 0.5);
701  psd1.add<edm::InputTag>("Producer", edm::InputTag("hpsPFTauDiscriminationByLooseChargedIsolation"));
702  pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("preIso", psd1);
703  }
704  desc.add<edm::ParameterSetDescription>("Prediscriminants", pset_Prediscriminants);
705  }
706 
707  desc.add<int>("verbosity", 0);
708 
709  desc.add<bool>("deltaBetaPUTrackPtCutOverride", false);
710  desc.add<bool>("applyRhoCorrection", false);
711 
712  desc.add<double>("WeightECALIsolation", 1.0);
713  desc.add<double>("rhoUEOffsetCorrection", 1.0);
714  desc.add<double>("deltaBetaPUTrackPtCutOverride_val", -1.5);
715  desc.add<double>("isoConeSizeForDeltaBeta", 0.5);
716  desc.add<double>("customOuterCone", -1.0);
717  desc.add<edm::InputTag>("particleFlowSrc", edm::InputTag("particleFlow"));
718 
719  // options for various stored ID raw values
720  edm::ParameterSetDescription desc_idlist;
721  desc_idlist.add<string>("IDname"); //not needed by producer but required for mapping at PAT level
722  desc_idlist.add<bool>("storeRawSumPt", false);
723  desc_idlist.add<bool>("storeRawPUsumPt", false);
724  desc_idlist.add<bool>("storeRawOccupancy", false);
725  desc_idlist.add<bool>("storeRawFootprintCorrection", false);
726  desc_idlist.add<bool>("storeRawPhotonSumPt_outsideSignalCone", false);
727  desc_idlist.add<bool>("ApplyDiscriminationByECALIsolation", false);
728  desc_idlist.add<bool>("ApplyDiscriminationByWeightedECALIsolation", false);
729  desc_idlist.add<bool>("ApplyDiscriminationByTrackerIsolation", false);
730  desc_idlist.add<bool>("applyDeltaBetaCorrection", false);
731  desc_idlist.add<bool>("UseAllPFCandsForWeights", false);
732  desc.addVPSet("IDdefinitions", desc_idlist, {});
733  // options for various stored ID WPs
734  edm::ParameterSetDescription desc_idwplist;
735  desc_idwplist.add<string>("IDname"); //not needed by producer but required for mapping at PAT level
736  desc_idwplist.add<std::vector<string>>("referenceRawIDNames")
737  ->setComment(
738  "List of raw IDs defined in 'IDdefinitions' to pass all respective conditions defined in "
739  "'maximumAbsoluteValues', 'maximumRelativeValues' , and 'relativeValueOffsets'");
740  desc_idwplist.add<std::vector<double>>("maximumAbsoluteValues", {});
741  desc_idwplist.add<std::vector<double>>("maximumRelativeValues", {});
742  desc_idwplist.add<std::vector<double>>("relativeValueOffsets", {});
743  desc.addVPSet("IDWPdefinitions", desc_idwplist, {});
744 
745  descriptions.add("pfRecoTauDiscriminationByIsolationContainer", desc);
746 }
747 
double weightedSum(const std::vector< CandidatePtr > &inColl_, double eta, double phi) const
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
ParameterDescriptionBase * addOptional(U const &iLabel, T const &value)
static void fillDescriptions(edm::ParameterSetDescription &descriptions)
Declare all parameters read from python config file.
std::vector< ParameterSet > VParameterSet
Definition: ParameterSet.h:35
selection
main part
Definition: corrVsCorr.py:100
ParameterSet const & getParameterSet(std::string const &) const
std::unique_ptr< tau::RecoTauQualityCuts > pileupQcutsPUTrackSelection_
std::unique_ptr< tau::RecoTauVertexAssociator > vertexAssociator_
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
const LorentzVector & p4() const final
four-momentum Lorentz vector
#define LogTrace(id)
std::vector< std::unique_ptr< FootprintCorrection > > footprintCorrections_
void push_back(D *&d)
Definition: OwnVector.h:326
reco::SingleTauDiscriminatorContainer discriminate(const PFTauRef &pfTau) const override
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
virtual int charge() const =0
electric charge
std::pair< edm::ParameterSet, edm::ParameterSet > factorizePUQCuts(const edm::ParameterSet &inputSet)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
constexpr auto deltaR2(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:16
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void beginEvent(const edm::Event &evt, const edm::EventSetup &evtSetup) override
edm::EDGetTokenT< edm::View< reco::Candidate > > pfCand_token
void add(std::string const &label, ParameterSetDescription const &psetDescription)
fixed size matrix
HLT enums.
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void setP4(const LorentzVector &p4) final
set 4-momentum
def move(src, dest)
Definition: eostools.py:511
void reserve(size_t)
Definition: OwnVector.h:320
Definition: event.py:1
#define LogDebug(id)