CMS 3D CMS Logo

TauDiscriminationAgainstElectronMVA6.cc
Go to the documentation of this file.
1 /* class TauDiscriminationAgainstElectronMVA6
2  * created : Nov 2 2015,
3  * revised : May 29 2020,
4  * Authors : Fabio Colombo (KIT)
5  * Anne-Catherine Le Bihan (IPHC),
6  * Michal Bluj (NCBJ)
7  */
8 
13 
14 template <class TauType, class TauDiscriminator, class ElectronType>
16  reco::TauDiscriminatorContainer,
17  reco::SingleTauDiscriminatorContainer,
18  TauDiscriminator> {
19 public:
20  typedef std::vector<TauType> TauCollection;
22  typedef std::vector<ElectronType> ElectronCollection;
23 
27  reco::SingleTauDiscriminatorContainer,
28  TauDiscriminator>::TauDiscriminationProducerBase(cfg),
29  moduleLabel_(cfg.getParameter<std::string>("@module_label")),
30  mva_(
31  std::make_unique<AntiElectronIDMVA6<TauType, ElectronType>>(cfg, edm::EDConsumerBase::consumesCollector())),
32  Electron_token(edm::EDConsumerBase::consumes<ElectronCollection>(
33  cfg.getParameter<edm::InputTag>("srcElectrons"))), // MB: full specification with prefix mandatory
34  positionAtECalEntrance_(PositionAtECalEntranceComputer(edm::EDConsumerBase::consumesCollector(),
35  cfg.getParameter<bool>("isPhase2"))),
36  vetoEcalCracks_(cfg.getParameter<bool>("vetoEcalCracks")),
37  isPhase2_(cfg.getParameter<bool>("isPhase2")),
38  verbosity_(cfg.getParameter<int>("verbosity")) {
39  deltaREleTauMax_ = (isPhase2_ ? 0.2 : 0.3);
40  }
41 
42  void beginEvent(const edm::Event& evt, const edm::EventSetup& es) override {
43  mva_->beginEvent(evt, es);
45  evt.getByToken(this->Tau_token, taus_);
47  }
48 
50 
52 
53  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
54 
55 private:
56  bool isInEcalCrack(double) const;
57 
58  // Overloaded method with explicit type specification to avoid partial
59  //implementation of full class
60  std::pair<float, float> getTauEtaAtECalEntrance(const reco::PFTauRef& theTauRef) const;
61  std::pair<float, float> getTauEtaAtECalEntrance(const pat::TauRef& theTauRef) const;
62 
64  std::unique_ptr<AntiElectronIDMVA6<TauType, ElectronType>> mva_;
65 
69 
71 
72  static constexpr float ecalBarrelEndcapEtaBorder_ = 1.479;
73  static constexpr float ecalEndcapVFEndcapEtaBorder_ = 2.4;
74 
76 
77  bool isPhase2_;
79 
81 };
82 
83 template <class TauType, class TauDiscriminator, class ElectronType>
86  const TauRef& theTauRef) const {
88  result.rawValues = {1., -1.};
89  double category = -1.;
90  bool isGsfElectronMatched = false;
91 
92  double deltaRDummy = 9.9;
93 
94  std::pair<float, float> tauEtaAtECalEntrance;
96  tauEtaAtECalEntrance = getTauEtaAtECalEntrance(theTauRef);
97  else
98  throw cms::Exception("TauDiscriminationAgainstElectronMVA6")
99  << "Unsupported TauType used. You must use either reco::PFTau or pat::Tau.";
100 
101  if ((*theTauRef).leadChargedHadrCand().isNonnull()) {
102  int numSignalGammaCandsInSigCone = 0;
103  double signalrad = std::clamp(3.0 / std::max(1.0, theTauRef->pt()), 0.05, 0.10);
104  for (const auto& gamma : theTauRef->signalGammaCands()) {
105  double dR = deltaR(gamma->p4(), theTauRef->leadChargedHadrCand()->p4());
106  // pfGammas inside the tau signal cone
107  if (dR < signalrad) {
108  numSignalGammaCandsInSigCone += 1;
109  }
110  }
111 
112  bool hasGsfTrack = false;
113  const reco::CandidatePtr& leadChCand = theTauRef->leadChargedHadrCand();
114  if (leadChCand.isNonnull()) {
115  if (isPhase2_) {
116  //MB: for phase-2 has gsf-track reads lead charged cand is pf-electron
117  hasGsfTrack = (std::abs(leadChCand->pdgId()) == 11);
118  } else {
119  const pat::PackedCandidate* packedLeadChCand = dynamic_cast<const pat::PackedCandidate*>(leadChCand.get());
120  if (packedLeadChCand != nullptr) {
121  hasGsfTrack = (std::abs(packedLeadChCand->pdgId()) == 11);
122  } else {
123  const reco::PFCandidate* pfLeadChCand = dynamic_cast<const reco::PFCandidate*>(leadChCand.get());
124  //pfLeadChCand can not be a nullptr here as it would be imply taus not built either with PFCandidates or PackedCandidates
125  hasGsfTrack = pfLeadChCand->gsfTrackRef().isNonnull();
126  }
127  }
128  }
129 
130  // loop over the electrons
131  size_t iElec = 0;
132  for (const auto& theElectron : *electrons_) {
133  edm::Ref<ElectronCollection> theElecRef(electrons_, iElec);
134  iElec++;
135  if (theElectron.pt() > 10.) { // CV: only take electrons above some minimal energy/Pt into account...
136  double deltaREleTau = deltaR(theElectron.p4(), theTauRef->p4());
137  deltaRDummy = std::min(deltaREleTau, deltaRDummy);
138  if (deltaREleTau < deltaREleTauMax_) {
139  double mva_match = mva_->mvaValue(*theTauRef, theElecRef);
140  if (!hasGsfTrack)
141  hasGsfTrack = theElectron.gsfTrack().isNonnull();
142 
143  // veto taus that go to ECal crack
144  if (vetoEcalCracks_ &&
145  (isInEcalCrack(tauEtaAtECalEntrance.first) || isInEcalCrack(tauEtaAtECalEntrance.second))) {
146  // add category index
147  result.rawValues.at(1) = category;
148  // return MVA output value
149  result.rawValues.at(0) = -99.;
150  return result;
151  }
152  // veto taus that go to ECal crack
153 
154  if (std::abs(tauEtaAtECalEntrance.first) < ecalBarrelEndcapEtaBorder_) { // Barrel
155  if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
156  category = 5.;
157  } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
158  category = 7.;
159  }
160  } else if (!isPhase2_ || std::abs(tauEtaAtECalEntrance.first) < ecalEndcapVFEndcapEtaBorder_) { // Endcap
161  if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
162  category = 13.;
163  } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
164  category = 15.;
165  }
166  } else { // VeryForwardEndcap
167  if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
168  category = 14.;
169  } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
170  category = 16.;
171  }
172  }
173 
174  result.rawValues.at(0) = std::min(result.rawValues.at(0), float(mva_match));
175  isGsfElectronMatched = true;
176  } // deltaR < deltaREleTauMax_
177  } // electron pt > 10
178  } // end of loop over electrons
179 
180  if (!isGsfElectronMatched) {
181  double mva_nomatch = mva_->mvaValue(*theTauRef);
182 
183  // veto taus that go to ECal crack
184  if (vetoEcalCracks_ &&
185  (isInEcalCrack(tauEtaAtECalEntrance.first) || isInEcalCrack(tauEtaAtECalEntrance.second))) {
186  // add category index
187  result.rawValues.at(1) = category;
188  // return MVA output value
189  result.rawValues.at(0) = -99.;
190  return result;
191  }
192  // veto taus that go to ECal crack
193 
194  if (std::abs(tauEtaAtECalEntrance.first) < ecalBarrelEndcapEtaBorder_) { // Barrel
195  if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
196  category = 0.;
197  } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
198  category = 2.;
199  }
200  } else if (!isPhase2_ || std::abs(tauEtaAtECalEntrance.first) < ecalEndcapVFEndcapEtaBorder_) { // Endcap
201  if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
202  category = 8.;
203  } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
204  category = 10.;
205  }
206  } else { // VeryForwardEndcap
207  if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
208  category = 9.;
209  } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
210  category = 11.;
211  }
212  }
213 
214  result.rawValues.at(0) = std::min(result.rawValues.at(0), float(mva_nomatch));
215  }
216  }
217 
218  if (verbosity_) {
219  edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
220  << "<" + this->getTauTypeString() + "AgainstElectronMVA6::discriminate>:";
221  edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
222  << " tau: Pt = " << theTauRef->pt() << ", eta = " << theTauRef->eta() << ", phi = " << theTauRef->phi();
223  edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
224  << " deltaREleTau = " << deltaRDummy << ", isGsfElectronMatched = " << isGsfElectronMatched;
225  edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
226  << " #Prongs = " << theTauRef->signalChargedHadrCands().size();
227  edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
228  << " MVA = " << result.rawValues.at(0) << ", category = " << category;
229  }
230 
231  // add category index
232  result.rawValues.at(1) = category;
233  // return MVA output value
234  return result;
235 }
236 
237 template <class TauType, class TauDiscriminator, class ElectronType>
239  double absEta = std::abs(eta);
240  return (absEta > 1.460 && absEta < 1.558);
241 }
242 
243 template <class TauType, class TauDiscriminator, class ElectronType>
244 std::pair<float, float>
246  const reco::PFTauRef& theTauRef) const {
247  float tauEtaAtECalEntrance = -99;
248  float leadChargedCandEtaAtECalEntrance = -99;
249  float sumEtaTimesEnergy = 0;
250  float sumEnergy = 0;
251  float leadChargedCandPt = -99;
252 
253  for (const auto& candidate : theTauRef->signalCands()) {
254  float etaAtECalEntrance = candidate->eta();
255  const reco::Track* track = nullptr;
256  const reco::PFCandidate* pfCandidate = dynamic_cast<const reco::PFCandidate*>(candidate.get());
257  if (pfCandidate != nullptr) {
258  if (!isPhase2_ || std::abs(theTauRef->eta()) < ecalBarrelEndcapEtaBorder_) { // ECal
259  etaAtECalEntrance = pfCandidate->positionAtECALEntrance().eta();
260  } else { // HGCal
261  bool success = false;
262  reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
263  if (success) {
264  etaAtECalEntrance = posAtECal.eta();
265  }
266  }
267  if (pfCandidate->trackRef().isNonnull())
268  track = pfCandidate->trackRef().get();
269  else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull())
270  track = pfCandidate->muonRef()->innerTrack().get();
271  else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull())
272  track = pfCandidate->muonRef()->globalTrack().get();
273  else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull())
274  track = pfCandidate->muonRef()->outerTrack().get();
275  else if (pfCandidate->gsfTrackRef().isNonnull())
276  track = pfCandidate->gsfTrackRef().get();
277  } else {
278  bool success = false;
279  reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
280  if (success) {
281  etaAtECalEntrance = posAtECal.eta();
282  }
283  track = candidate->bestTrack();
284  }
285  if (track != nullptr) {
286  if (track->pt() > leadChargedCandPt) {
287  leadChargedCandEtaAtECalEntrance = etaAtECalEntrance;
288  leadChargedCandPt = track->pt();
289  }
290  }
291  sumEtaTimesEnergy += etaAtECalEntrance * candidate->energy();
292  sumEnergy += candidate->energy();
293  }
294  if (sumEnergy > 0.) {
295  tauEtaAtECalEntrance = sumEtaTimesEnergy / sumEnergy;
296  }
297  return std::pair<float, float>(tauEtaAtECalEntrance, leadChargedCandEtaAtECalEntrance);
298 }
299 
300 template <class TauType, class TauDiscriminator, class ElectronType>
301 std::pair<float, float>
303  const pat::TauRef& theTauRef) const {
304  if (!isPhase2_ || std::abs(theTauRef->eta()) < ecalBarrelEndcapEtaBorder_) { // ECal
305  return std::pair<float, float>(theTauRef->etaAtEcalEntrance(), theTauRef->etaAtEcalEntranceLeadChargedCand());
306  } else { // HGCal
307  float tauEtaAtECalEntrance = -99;
308  float leadChargedCandEtaAtECalEntrance = -99;
309  float sumEtaTimesEnergy = 0.;
310  float sumEnergy = 0.;
311  float leadChargedCandPt = -99;
312 
313  for (const auto& candidate : theTauRef->signalCands()) {
314  float etaAtECalEntrance = candidate->eta();
315  bool success = false;
316  reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
317  if (success) {
318  etaAtECalEntrance = posAtECal.eta();
319  }
320  const reco::Track* track = candidate->bestTrack();
321  if (track != nullptr) {
322  if (track->pt() > leadChargedCandPt) {
323  leadChargedCandEtaAtECalEntrance = etaAtECalEntrance;
324  leadChargedCandPt = track->pt();
325  }
326  }
327  sumEtaTimesEnergy += etaAtECalEntrance * candidate->energy();
328  sumEnergy += candidate->energy();
329  }
330  if (sumEnergy > 0.) {
331  tauEtaAtECalEntrance = sumEtaTimesEnergy / sumEnergy;
332  }
333  return std::pair<float, float>(tauEtaAtECalEntrance, leadChargedCandEtaAtECalEntrance);
334  }
335 }
336 
337 template <class TauType, class TauDiscriminator, class ElectronType>
339  edm::ConfigurationDescriptions& descriptions) {
340  // {pfReco,pat}TauDiscriminationAgainstElectronMVA6
342 
343  desc.add<std::string>("method", "BDTG");
344  desc.add<bool>("loadMVAfromDB", true);
345  desc.add<bool>("returnMVA", true);
346 
347  desc.add<std::string>("mvaName_NoEleMatch_woGwoGSF_BL", "gbr_NoEleMatch_woGwoGSF_BL");
348  desc.add<std::string>("mvaName_NoEleMatch_wGwoGSF_BL", "gbr_NoEleMatch_wGwoGSF_BL");
349  desc.add<std::string>("mvaName_woGwGSF_BL", "gbr_woGwGSF_BL");
350  desc.add<std::string>("mvaName_wGwGSF_BL", "gbr_wGwGSF_BL");
351  desc.add<std::string>("mvaName_NoEleMatch_woGwoGSF_EC", "gbr_NoEleMatch_woGwoGSF_EC");
352  desc.add<std::string>("mvaName_NoEleMatch_wGwoGSF_EC", "gbr_NoEleMatch_wGwoGSF_EC");
353  desc.add<std::string>("mvaName_woGwGSF_EC", "gbr_woGwGSF_EC");
354  desc.add<std::string>("mvaName_wGwGSF_EC", "gbr_wGwGSF_EC");
355 
356  desc.add<double>("minMVANoEleMatchWOgWOgsfBL", 0.0);
357  desc.add<double>("minMVANoEleMatchWgWOgsfBL", 0.0);
358  desc.add<double>("minMVAWOgWgsfBL", 0.0);
359  desc.add<double>("minMVAWgWgsfBL", 0.0);
360  desc.add<double>("minMVANoEleMatchWOgWOgsfEC", 0.0);
361  desc.add<double>("minMVANoEleMatchWgWOgsfEC", 0.0);
362  desc.add<double>("minMVAWOgWgsfEC", 0.0);
363  desc.add<double>("minMVAWgWgsfEC", 0.0);
364 
365  desc.ifValue(
366  edm::ParameterDescription<bool>("isPhase2", false, true),
367  // MB: "srcElectrons" present for both phase-2 and non-phase2 to have a non-empy case for default, i.e. isPhase2=false
368  false >> (edm::ParameterDescription<edm::InputTag>("srcElectrons", edm::InputTag("fixme"), true)) or
369  // The following used only for Phase2
370  true >> (edm::ParameterDescription<edm::InputTag>("srcElectrons", edm::InputTag("fixme"), true) and
371  edm::ParameterDescription<std::string>("mvaName_wGwGSF_VFEC", "gbr_wGwGSF_VFEC", true) and
372  edm::ParameterDescription<std::string>("mvaName_woGwGSF_VFEC", "gbr_woGwGSF_VFEC", true) and
374  "mvaName_NoEleMatch_wGwoGSF_VFEC", "gbr_NoEleMatch_wGwoGSF_VFEC", true) and
376  "mvaName_NoEleMatch_woGwoGSF_VFEC", "gbr_NoEleMatch_woGwoGSF_VFEC", true) and
377  edm::ParameterDescription<double>("minMVAWOgWgsfVFEC", 0.0, true) and
378  edm::ParameterDescription<double>("minMVAWgWgsfVFEC", 0.0, true) and
379  edm::ParameterDescription<double>("minMVANoEleMatchWgWOgsfVFEC", 0.0, true) and
380  edm::ParameterDescription<double>("minMVANoEleMatchWOgWOgsfVFEC", 0.0, true)));
381 
382  // Relevant only for gsfElectrons for Phase2
384  desc.add<std::vector<edm::InputTag>>("hgcalElectronIDs", std::vector<edm::InputTag>())
385  ->setComment("Relevant only for Phase-2");
386  }
387  desc.add<bool>("vetoEcalCracks", true);
388  desc.add<bool>("usePhiAtEcalEntranceExtrapolation", false);
389  desc.add<int>("verbosity", 0);
390 
394  TauDiscriminator>::fillProducerDescriptions(desc); // inherited from the base-class
395 
396  descriptions.addWithDefaultLabel(desc);
397 }
398 
403 
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
TauDiscriminationAgainstElectronMVA6(const edm::ParameterSet &cfg)
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
reco::GsfTrackRef gsfTrackRef() const
Definition: PFCandidate.cc:469
void beginEvent(const edm::Event &evt, const edm::EventSetup &es) override
void beginEvent(const edm::EventSetup &)
TauDiscriminationAgainstElectronMVA6< reco::PFTau, reco::PFTauDiscriminator, reco::GsfElectron > PFRecoTauDiscriminationAgainstElectronMVA6
std::pair< float, float > getTauEtaAtECalEntrance(const reco::PFTauRef &theTauRef) const
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
const math::XYZPointF & positionAtECALEntrance() const
Definition: PFCandidate.h:388
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:528
int pdgId() const override
PDG identifier.
std::unique_ptr< AntiElectronIDMVA6< TauType, ElectronType > > mva_
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
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
reco::SingleTauDiscriminatorContainer discriminate(const TauRef &) const override
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
TauDiscriminationAgainstElectronMVA6< pat::Tau, pat::PATTauDiscriminator, pat::Electron > PATTauDiscriminationAgainstElectronMVA6
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Log< level::Warning, true > LogPrint
reco::MuonRef muonRef() const
Definition: PFCandidate.cc:450
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:139
edm::ValueMap< SingleTauDiscriminatorContainer > TauDiscriminatorContainer
edm::EDGetTokenT< ElectronCollection > Electron_token
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
fixed size matrix
HLT enums.
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
reco::TrackRef trackRef() const
Definition: PFCandidate.cc:437
math::XYZPoint Point
point in the space
Definition: Candidate.h:40