CMS 3D CMS Logo

PFRecoTauDiscriminationAgainstMuonSimple.cc
Go to the documentation of this file.
1 
14 
16 
24 
26 
27 #include <vector>
28 #include <string>
29 #include <iostream>
30 #include <atomic>
31 
33 
34 namespace {
35 
36  struct PFRecoTauDiscriminationAgainstMuonSimpleConfigSet {
37  PFRecoTauDiscriminationAgainstMuonSimpleConfigSet(double hop, int mNOM, bool doCMV, int mNHL2S, int mNSTA, int mNRPC)
38  : hop(hop),
39  maxNumberOfMatches(mNOM),
40  doCaloMuonVeto(doCMV),
42  maxNumberOfSTAMuons(mNSTA),
43  maxNumberOfRPCMuons(mNRPC) {}
44 
45  double hop;
47  bool doCaloMuonVeto;
50  };
51 
52  class PFRecoTauDiscriminationAgainstMuonSimple final : public PFTauDiscriminationContainerProducerBase {
53  public:
54  explicit PFRecoTauDiscriminationAgainstMuonSimple(const edm::ParameterSet& cfg)
55  : PFTauDiscriminationContainerProducerBase(cfg), moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
56  auto const wpDefs = cfg.getParameter<std::vector<edm::ParameterSet>>("IDWPdefinitions");
57  for (auto& wpDef : wpDefs) {
58  wpDefs_.push_back(
59  PFRecoTauDiscriminationAgainstMuonSimpleConfigSet(wpDef.getParameter<double>("HoPMin"),
60  wpDef.getParameter<int>("maxNumberOfMatches"),
61  wpDef.getParameter<bool>("doCaloMuonVeto"),
62  wpDef.getParameter<int>("maxNumberOfHitsLast2Stations"),
63  wpDef.getParameter<int>("maxNumberOfSTAMuons"),
64  wpDef.getParameter<int>("maxNumberOfRPCMuons")));
65  }
66  srcPatMuons_ = cfg.getParameter<edm::InputTag>("srcPatMuons");
67  muons_token = consumes<pat::MuonCollection>(srcPatMuons_);
68  dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
69  dRmuonMatchLimitedToJetArea_ = cfg.getParameter<bool>("dRmuonMatchLimitedToJetArea");
70  minPtMatchedMuon_ = cfg.getParameter<double>("minPtMatchedMuon");
71  typedef std::vector<int> vint;
72  maskMatchesDT_ = cfg.getParameter<vint>("maskMatchesDT");
73  maskMatchesCSC_ = cfg.getParameter<vint>("maskMatchesCSC");
74  maskMatchesRPC_ = cfg.getParameter<vint>("maskMatchesRPC");
75  maskHitsDT_ = cfg.getParameter<vint>("maskHitsDT");
76  maskHitsCSC_ = cfg.getParameter<vint>("maskHitsCSC");
77  maskHitsRPC_ = cfg.getParameter<vint>("maskHitsRPC");
78  verbosity_ = cfg.getParameter<int>("verbosity");
79  }
80  ~PFRecoTauDiscriminationAgainstMuonSimple() override {}
81 
82  void beginEvent(const edm::Event&, const edm::EventSetup&) override;
83 
85 
86  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
87 
88  private:
90  std::vector<PFRecoTauDiscriminationAgainstMuonSimpleConfigSet> wpDefs_;
91  edm::InputTag srcPatMuons_;
94  double dRmuonMatch_;
95  bool dRmuonMatchLimitedToJetArea_;
96  double minPtMatchedMuon_;
97  std::vector<int> maskMatchesDT_;
98  std::vector<int> maskMatchesCSC_;
99  std::vector<int> maskMatchesRPC_;
100  std::vector<int> maskHitsDT_;
101  std::vector<int> maskHitsCSC_;
102  std::vector<int> maskHitsRPC_;
103 
104  static std::atomic<unsigned int> numWarnings_;
105  static constexpr unsigned int maxWarnings_ = 3;
106  int verbosity_;
107  };
108 
109  std::atomic<unsigned int> PFRecoTauDiscriminationAgainstMuonSimple::numWarnings_{0};
110 
112  evt.getByToken(muons_token, muons_);
113  }
114 
116  const reco::PFTauRef& pfTau) const {
117  if (verbosity_) {
118  edm::LogPrint("PFTauAgainstMuonSimple") << "<PFRecoTauDiscriminationAgainstMuonSimple::discriminate>:";
119  edm::LogPrint("PFTauAgainstMuonSimple") << " moduleLabel = " << moduleLabel_;
120  edm::LogPrint("PFTauAgainstMuonSimple")
121  << "tau #" << pfTau.key() << ": Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
122  << ", phi = " << pfTau->phi() << ", decay mode = " << pfTau->decayMode();
123  }
124 
125  //energy fraction carried by leading charged hadron
126  const reco::CandidatePtr& pfLeadChargedHadron = pfTau->leadChargedHadrCand();
127  double caloEnergyFraction = 99;
128  if (pfLeadChargedHadron.isNonnull()) {
129  const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(pfLeadChargedHadron.get());
130  if (pCand != nullptr) {
131  caloEnergyFraction = pCand->caloFraction();
132  if (pCand->hasTrackDetails()) //MB: relate energy fraction to track momentum rather than to energy of candidate
133  caloEnergyFraction *= pCand->energy() / pCand->bestTrack()->p();
134  if (verbosity_) {
135  edm::LogPrint("PFTauAgainstMuonSimple")
136  << "decayMode = " << pfTau->decayMode() << ", caloEnergy(ECAL+HCAL)Fraction = " << caloEnergyFraction
137  << ", leadPFChargedHadronP = " << pCand->p() << ", leadPFChargedHadron pdgId = " << pCand->pdgId();
138  }
139  }
140  }
141  //iterate over muons to find ones to use for discrimination
142  std::vector<const pat::Muon*> muonsToCheck;
143  size_t numMuons = muons_->size();
144  for (size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon) {
145  bool matched = false;
146  const pat::MuonRef muon(muons_, idxMuon);
147  if (verbosity_)
148  edm::LogPrint("PFTauAgainstMuonSimple") << "muon #" << muon.key() << ": Pt = " << muon->pt()
149  << ", eta = " << muon->eta() << ", phi = " << muon->phi();
150  //firsty check if the muon corrsponds with the leading tau particle
151  for (size_t iCand = 0; iCand < muon->numberOfSourceCandidatePtrs(); ++iCand) {
152  const reco::CandidatePtr& srcCand = muon->sourceCandidatePtr(iCand);
153  if (srcCand.isNonnull() && pfLeadChargedHadron.isNonnull() && srcCand.id() == pfLeadChargedHadron.id() &&
154  srcCand.key() == pfLeadChargedHadron.key()) {
155  muonsToCheck.push_back(muon.get());
156  matched = true;
157  if (verbosity_)
158  edm::LogPrint("PFTauAgainstMuonSimple") << " tau has muonRef.";
159  break;
160  }
161  }
162  if (matched)
163  continue;
164  //check dR matching
165  if (!(muon->pt() > minPtMatchedMuon_)) {
166  if (verbosity_) {
167  edm::LogPrint("PFTauAgainstMuonSimple") << " fails Pt cut --> skipping it.";
168  }
169  continue;
170  }
171  double dR = deltaR(muon->p4(), pfTau->p4());
172  double dRmatch = dRmuonMatch_;
173  if (dRmuonMatchLimitedToJetArea_) {
174  double jetArea = 0.;
175  if (pfTau->jetRef().isNonnull())
176  jetArea = pfTau->jetRef()->jetArea();
177  if (jetArea > 0.) {
178  dRmatch = std::min(dRmatch, std::sqrt(jetArea / M_PI));
179  } else {
180  if (numWarnings_ < maxWarnings_) {
181  edm::LogInfo("PFRecoTauDiscriminationAgainstMuonSimple::discriminate")
182  << "Jet associated to Tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
183  << ", phi = " << pfTau->phi() << " has area = " << jetArea << " !!";
184  ++numWarnings_;
185  }
186  dRmatch = pfTau->signalConeSize();
187  }
188  dRmatch = std::max(dRmatch, pfTau->signalConeSize());
189  if (dR < dRmatch) {
190  if (verbosity_)
191  edm::LogPrint("PFTauAgainstMuonSimple") << " overlaps with tau, dR = " << dR;
192  muonsToCheck.push_back(muon.get());
193  }
194  }
195  }
196  //now examine selected muons
197  std::vector<int> numMatchesDT(4);
198  std::vector<int> numMatchesCSC(4);
199  std::vector<int> numMatchesRPC(4);
200  std::vector<int> numHitsDT(4);
201  std::vector<int> numHitsCSC(4);
202  std::vector<int> numHitsRPC(4);
203  //MB: clear counters of matched segments and hits globally as in the AgainstMuon2 discriminant, but note that they will sum for all matched muons. However it is not likely that there is more than one matched muon.
204  for (int iStation = 0; iStation < 4; ++iStation) {
205  numMatchesDT[iStation] = 0;
206  numMatchesCSC[iStation] = 0;
207  numMatchesRPC[iStation] = 0;
208  numHitsDT[iStation] = 0;
209  numHitsCSC[iStation] = 0;
210  numHitsRPC[iStation] = 0;
211  }
212  int numSTAMuons = 0, numRPCMuons = 0;
213  int numStationsWithMatches = 0;
214  int numLast2StationsWithHits = 0;
215  if (verbosity_ && !muonsToCheck.empty())
216  edm::LogPrint("PFTauAgainstMuonSimple") << "Muons to check (" << muonsToCheck.size() << "):";
217  size_t iMu = 0;
218  for (const auto& mu : muonsToCheck) {
219  if (mu->isStandAloneMuon())
220  numSTAMuons++;
221  if (mu->muonID("RPCMuLoose"))
222  numRPCMuons++;
223  reco::tau::countMatches(*mu, numMatchesDT, numMatchesCSC, numMatchesRPC);
224  for (int iStation = 0; iStation < 4; ++iStation) {
225  if (numMatchesDT[iStation] > 0 && !maskMatchesDT_[iStation])
226  ++numStationsWithMatches;
227  if (numMatchesCSC[iStation] > 0 && !maskMatchesCSC_[iStation])
228  ++numStationsWithMatches;
229  if (numMatchesRPC[iStation] > 0 && !maskMatchesRPC_[iStation])
230  ++numStationsWithMatches;
231  }
232  reco::tau::countHits(*mu, numHitsDT, numHitsCSC, numHitsRPC);
233  for (int iStation = 2; iStation < 4; ++iStation) {
234  if (numHitsDT[iStation] > 0 && !maskHitsDT_[iStation])
235  ++numLast2StationsWithHits;
236  if (numHitsCSC[iStation] > 0 && !maskHitsCSC_[iStation])
237  ++numLast2StationsWithHits;
238  if (numHitsRPC[iStation] > 0 && !maskHitsRPC_[iStation])
239  ++numLast2StationsWithHits;
240  }
241  if (verbosity_)
242  edm::LogPrint("PFTauAgainstMuonSimple")
243  << "\t" << iMu << ": Pt = " << mu->pt() << ", eta = " << mu->eta() << ", phi = " << mu->phi() << "\n\t"
244  << " isSTA: " << mu->isStandAloneMuon() << ", isRPCLoose: " << mu->muonID("RPCMuLoose")
245  << "\n\t numMatchesDT = " << format_vint(numMatchesDT)
246  << "\n\t numMatchesCSC = " << format_vint(numMatchesCSC)
247  << "\n\t numMatchesRPC = " << format_vint(numMatchesRPC)
248  << "\n\t --> numStationsWithMatches = " << numStationsWithMatches
249  << "\n\t numHitsDT = " << format_vint(numHitsDT) << "\n\t numHitsCSC = " << format_vint(numHitsCSC)
250  << "\n\t numHitsRPC = " << format_vint(numHitsRPC)
251  << "\n\t --> numLast2StationsWithHits = " << numLast2StationsWithHits;
252  iMu++;
253  }
254 
256  for (auto const& wpDef : wpDefs_) {
257  bool discriminatorValue = true;
258  if (wpDef.maxNumberOfMatches >= 0 && numStationsWithMatches > wpDef.maxNumberOfMatches)
259  discriminatorValue = false;
260  if (wpDef.maxNumberOfHitsLast2Stations >= 0 && numLast2StationsWithHits > wpDef.maxNumberOfHitsLast2Stations)
261  discriminatorValue = false;
262  if (wpDef.maxNumberOfSTAMuons >= 0 && numSTAMuons > wpDef.maxNumberOfSTAMuons) {
263  discriminatorValue = false;
264  }
265  if (wpDef.maxNumberOfRPCMuons >= 0 && numRPCMuons > wpDef.maxNumberOfRPCMuons) {
266  discriminatorValue = false;
267  }
268  bool passesCaloMuonVeto = true;
269  if (pfTau->decayMode() == 0 && caloEnergyFraction < wpDef.hop) {
270  passesCaloMuonVeto = false;
271  }
272  if (wpDef.doCaloMuonVeto && !passesCaloMuonVeto) {
273  discriminatorValue = false;
274  }
275  result.workingPoints.push_back(discriminatorValue);
276  if (verbosity_)
277  edm::LogPrint("PFTauAgainstMuonSimple") << "--> returning discriminatorValue = " << result.workingPoints.back();
278  }
279  return result;
280  }
281 
283  // pfRecoTauDiscriminationAgainstMuonSimple
285  desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfRecoTauProducer"));
286  {
288  psd0.add<std::string>("BooleanOperator", "and");
289  {
291  psd1.add<double>("cut"); //, 0.5);
292  psd1.add<edm::InputTag>("Producer"); //, edm::InputTag("pfRecoTauDiscriminationByLeadingTrackFinding"));
293  psd0.addOptional<edm::ParameterSetDescription>("leadTrack", psd1); // optional with default?
294  }
295  // Prediscriminants can be
296  // Prediscriminants = noPrediscriminants,
297  // as in RecoTauTag/Configuration/python/HPSPFTaus_cff.py
298  //
299  // and the definition is:
300  // RecoTauTag/RecoTau/python/TauDiscriminatorTools.py
301  // # Require no prediscriminants
302  // noPrediscriminants = cms.PSet(
303  // BooleanOperator = cms.string("and"),
304  // )
305  // -- so this is the minimum required definition
306  // otherwise it inserts the leadTrack with Producer = InpuTag(...) and does not find the corresponding output during the run
307  desc.add<edm::ParameterSetDescription>("Prediscriminants", psd0);
308  }
309  {
311  desc_wp.add<std::string>("IDname");
312  desc_wp.add<double>("HoPMin");
313  desc_wp.add<int>("maxNumberOfMatches")->setComment("negative value would turn off this cut");
314  desc_wp.add<bool>("doCaloMuonVeto");
315  desc_wp.add<int>("maxNumberOfHitsLast2Stations")->setComment("negative value would turn off this cut");
316  desc_wp.add<int>("maxNumberOfSTAMuons")->setComment("negative value would turn off this cut");
317  desc_wp.add<int>("maxNumberOfRPCMuons")->setComment("negative value would turn off this cut");
318  edm::ParameterSet pset_wp;
319  pset_wp.addParameter<std::string>("IDname", "ByLooseMuonRejectionSimple");
320  pset_wp.addParameter<double>("HoPMin", 0.2);
321  pset_wp.addParameter<int>("maxNumberOfMatches", 1);
322  pset_wp.addParameter<bool>("doCaloMuonVeto", true);
323  pset_wp.addParameter<int>("maxNumberOfHitsLast2Stations", -1);
324  pset_wp.addParameter<int>("maxNumberOfSTAMuons", -1);
325  pset_wp.addParameter<int>("maxNumberOfRPCMuons", -1);
326  std::vector<edm::ParameterSet> vpsd_wp;
327  vpsd_wp.push_back(pset_wp);
328  desc.addVPSet("IDWPdefinitions", desc_wp, vpsd_wp);
329  }
330  //add empty raw value config to simplify subsequent provenance searches
331  desc.addVPSet("IDdefinitions", edm::ParameterSetDescription(), {});
332 
333  desc.add<edm::InputTag>("srcPatMuons", edm::InputTag("slimmedMuons"));
334  desc.add<double>("dRmuonMatch", 0.3);
335  desc.add<bool>("dRmuonMatchLimitedToJetArea", false);
336  desc.add<double>("minPtMatchedMuon", 5.0);
337  desc.add<std::vector<int>>("maskMatchesDT", {0, 0, 0, 0});
338  desc.add<std::vector<int>>("maskMatchesCSC", {1, 0, 0, 0})
339  ->setComment(
340  "flags to mask/unmask DT, CSC and RPC chambers in individual muon stations. Segments and hits that are "
341  "present in that muon station are ignored in case the 'mask' is set to 1. Per default only the innermost "
342  "CSC "
343  "chamber is ignored, as it is affected by spurious hits in high pile-up events.");
344  desc.add<std::vector<int>>("maskMatchesRPC", {0, 0, 0, 0});
345  desc.add<std::vector<int>>("maskHitsDT", {0, 0, 0, 0});
346  desc.add<std::vector<int>>("maskHitsCSC", {0, 0, 0, 0});
347  desc.add<std::vector<int>>("maskHitsRPC", {0, 0, 0, 0});
348  desc.add<int>("verbosity", 0);
349  descriptions.add("pfRecoTauDiscriminationAgainstMuonSimple", desc);
350  }
351 
352 } // namespace
353 
354 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuonSimple);
ParameterDescriptionBase * addOptional(U const &iLabel, T const &value)
ProductID id() const
Accessor for product ID.
Definition: Ptr.h:160
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:540
void countMatches(const reco::Muon &muon, std::vector< int > &numMatchesDT, std::vector< int > &numMatchesCSC, std::vector< int > &numMatchesRPC)
float caloFraction() const
Set the fraction of ECAL+HCAL energy over candidate energy.
virtual TauDiscriminatorDataType discriminate(const TauRef &tau) const =0
int pdgId() const override
PDG identifier.
key_type key() const
Accessor for product key.
Definition: Ref.h:250
double vint[400]
double p() const override
magnitude of momentum vector
T sqrt(T t)
Definition: SSEVec.h:19
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
void addParameter(std::string const &name, T const &value)
Definition: ParameterSet.h:135
virtual void beginEvent(const edm::Event &, const edm::EventSetup &)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Log< level::Warning, true > LogPrint
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:139
#define M_PI
const reco::Track * bestTrack() const override
return a pointer to the track if present. otherwise, return a null pointer
Log< level::Info, false > LogInfo
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void add(std::string const &label, ParameterSetDescription const &psetDescription)
key_type key() const
Definition: Ptr.h:165
double energy() const override
energy
std::string format_vint(const std::vector< int > &vi)
void countHits(const reco::Muon &muon, std::vector< int > &numHitsDT, std::vector< int > &numHitsCSC, std::vector< int > &numHitsRPC)
bool hasTrackDetails() const
Return true if a bestTrack can be extracted from this Candidate.