CMS 3D CMS Logo
PFRecoTauDiscriminationAgainstMuonSimple.cc
Go to the documentation of this file.
1 
13 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
14 
15 #include "FWCore/Utilities/interface/Exception.h"
16 
17 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
18 #include "DataFormats/PatCandidates/interface/Muon.h"
19 #include "DataFormats/Common/interface/Handle.h"
20 #include "DataFormats/TrackReco/interface/HitPattern.h"
21 #include "DataFormats/TrackReco/interface/Track.h"
22 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
23 #include "DataFormats/Math/interface/deltaR.h"
24 
25 #include "RecoTauTag/RecoTau/interface/RecoTauMuonTools.h"
26 
27 #include <vector>
28 #include <string>
29 #include <iostream>
30 #include <atomic>
31 
32 using reco::tau::format_vint;
33 
34 namespace {
35 
36  class PFRecoTauDiscriminationAgainstMuonSimple final : public PFTauDiscriminationProducerBase {
37  public:
38  explicit PFRecoTauDiscriminationAgainstMuonSimple(const edm::ParameterSet& cfg)
39  : PFTauDiscriminationProducerBase(cfg), moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
40  hop_ = cfg.getParameter<double>("HoPMin");
41  doCaloMuonVeto_ = cfg.getParameter<bool>("doCaloMuonVeto");
42  srcPatMuons_ = cfg.getParameter<edm::InputTag>("srcPatMuons");
43  Muons_token = consumes<pat::MuonCollection>(srcPatMuons_);
44  dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
45  dRmuonMatchLimitedToJetArea_ = cfg.getParameter<bool>("dRmuonMatchLimitedToJetArea");
46  minPtMatchedMuon_ = cfg.getParameter<double>("minPtMatchedMuon");
47  maxNumberOfMatches_ = cfg.getParameter<int>("maxNumberOfMatches");
48  maxNumberOfHitsLast2Stations_ = cfg.getParameter<int>("maxNumberOfHitsLast2Stations");
49  typedef std::vector<int> vint;
50  maskMatchesDT_ = cfg.getParameter<vint>("maskMatchesDT");
51  maskMatchesCSC_ = cfg.getParameter<vint>("maskMatchesCSC");
52  maskMatchesRPC_ = cfg.getParameter<vint>("maskMatchesRPC");
53  maskHitsDT_ = cfg.getParameter<vint>("maskHitsDT");
54  maskHitsCSC_ = cfg.getParameter<vint>("maskHitsCSC");
55  maskHitsRPC_ = cfg.getParameter<vint>("maskHitsRPC");
56  maxNumberOfSTAMuons_ = cfg.getParameter<int>("maxNumberOfSTAMuons");
57  maxNumberOfRPCMuons_ = cfg.getParameter<int>("maxNumberOfRPCMuons");
58 
59  verbosity_ = cfg.exists("verbosity") ? cfg.getParameter<int>("verbosity") : 0;
60  }
61  ~PFRecoTauDiscriminationAgainstMuonSimple() override {}
62 
63  void beginEvent(const edm::Event&, const edm::EventSetup&) override;
64 
65  double discriminate(const reco::PFTauRef&) const override;
66 
67  private:
68  std::string moduleLabel_;
69  int discriminatorOption_;
70  double hop_;
71  bool doCaloMuonVeto_;
72  edm::InputTag srcPatMuons_;
73  edm::Handle<pat::MuonCollection> muons_;
74  edm::EDGetTokenT<pat::MuonCollection> Muons_token;
75  double dRmuonMatch_;
76  bool dRmuonMatchLimitedToJetArea_;
77  double minPtMatchedMuon_;
78  int maxNumberOfMatches_;
79  std::vector<int> maskMatchesDT_;
80  std::vector<int> maskMatchesCSC_;
81  std::vector<int> maskMatchesRPC_;
82  int maxNumberOfHitsLast2Stations_;
83  std::vector<int> maskHitsDT_;
84  std::vector<int> maskHitsCSC_;
85  std::vector<int> maskHitsRPC_;
86  int maxNumberOfSTAMuons_, maxNumberOfRPCMuons_;
87 
88  static std::atomic<unsigned int> numWarnings_;
89  static constexpr unsigned int maxWarnings_ = 3;
90  int verbosity_;
91  };
92 
93  std::atomic<unsigned int> PFRecoTauDiscriminationAgainstMuonSimple::numWarnings_{0};
94 
95  void PFRecoTauDiscriminationAgainstMuonSimple::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
96  evt.getByToken(Muons_token, muons_);
97  }
98 
99  double PFRecoTauDiscriminationAgainstMuonSimple::discriminate(const reco::PFTauRef& pfTau) const {
100  if (verbosity_) {
101  edm::LogPrint("PFTauAgainstMuonSimple") << "<PFRecoTauDiscriminationAgainstMuonSimple::discriminate>:";
102  edm::LogPrint("PFTauAgainstMuonSimple") << " moduleLabel = " << moduleLabel_;
103  edm::LogPrint("PFTauAgainstMuonSimple")
104  << "tau #" << pfTau.key() << ": Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
105  << ", phi = " << pfTau->phi() << ", decay mode = " << pfTau->decayMode();
106  }
107 
108  //if (pfTau->decayMode() >= 5) return true; //MB: accept all multi-prongs??
109 
110  const reco::CandidatePtr& pfLeadChargedHadron = pfTau->leadChargedHadrCand();
111  bool passesCaloMuonVeto = true;
112  if (pfLeadChargedHadron.isNonnull()) {
113  const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(pfLeadChargedHadron.get());
114  if (pCand != nullptr) {
115  double rawCaloEnergyFraction = pCand->rawCaloFraction();
116  //if ( !(rawCaloEnergyFraction > 0.) ) rawCaloEnergyFraction = 99; //MB: hack against cases when rawCaloEnergyFraction is not stored; it makes performance of the H/P cut rather poor
117  if (verbosity_) {
118  edm::LogPrint("PFTauAgainstMuonSimple")
119  << "decayMode = " << pfTau->decayMode()
120  << ", rawCaloEnergy(ECAL+HCAL)Fraction = " << rawCaloEnergyFraction
121  << ", leadPFChargedHadronP = " << pCand->p() << ", leadPFChargedHadron pdgId = " << pCand->pdgId();
122  }
123  if (pfTau->decayMode() == 0 && rawCaloEnergyFraction < hop_)
124  passesCaloMuonVeto = false;
125  }
126  }
127  if (doCaloMuonVeto_ && !passesCaloMuonVeto) {
128  if (verbosity_)
129  edm::LogPrint("PFTauAgainstMuonSimple") << "--> CaloMuonVeto failed, returning 0";
130  return 0.;
131  }
132 
133  //iterate over muons to find ones to use for discrimination
134  std::vector<const pat::Muon*> muonsToCheck;
135  size_t numMuons = muons_->size();
136  for (size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon) {
137  bool matched = false;
138  const pat::MuonRef muon(muons_, idxMuon);
139  if (verbosity_)
140  edm::LogPrint("PFTauAgainstMuonSimple") << "muon #" << muon.key() << ": Pt = " << muon->pt()
141  << ", eta = " << muon->eta() << ", phi = " << muon->phi();
142  //firsty check if the muon corrsponds with the leading tau particle
143  for (size_t iCand = 0; iCand < muon->numberOfSourceCandidatePtrs(); ++iCand) {
144  const reco::CandidatePtr& srcCand = muon->sourceCandidatePtr(iCand);
145  if (srcCand.isNonnull() && pfLeadChargedHadron.isNonnull() && srcCand.id() == pfLeadChargedHadron.id() &&
146  srcCand.key() == pfLeadChargedHadron.key()) {
147  muonsToCheck.push_back(muon.get());
148  matched = true;
149  if (verbosity_)
150  edm::LogPrint("PFTauAgainstMuonSimple") << " tau has muonRef.";
151  break;
152  }
153  }
154  if (matched)
155  continue;
156  //check dR matching
157  if (!(muon->pt() > minPtMatchedMuon_)) {
158  if (verbosity_) {
159  edm::LogPrint("PFTauAgainstMuonSimple") << " fails Pt cut --> skipping it.";
160  }
161  continue;
162  }
163  double dR = deltaR(muon->p4(), pfTau->p4());
164  double dRmatch = dRmuonMatch_;
165  if (dRmuonMatchLimitedToJetArea_) {
166  double jetArea = 0.;
167  if (pfTau->jetRef().isNonnull())
168  jetArea = pfTau->jetRef()->jetArea();
169  if (jetArea > 0.) {
170  dRmatch = TMath::Min(dRmatch, TMath::Sqrt(jetArea / TMath::Pi()));
171  } else {
172  if (numWarnings_ < maxWarnings_) {
173  edm::LogInfo("PFRecoTauDiscriminationAgainstMuonSimple::discriminate")
174  << "Jet associated to Tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
175  << ", phi = " << pfTau->phi() << " has area = " << jetArea << " !!";
176  ++numWarnings_;
177  }
178  dRmatch = pfTau->signalConeSize();
179  }
180  dRmatch = TMath::Max(dRmatch, pfTau->signalConeSize());
181  if (dR < dRmatch) {
182  if (verbosity_)
183  edm::LogPrint("PFTauAgainstMuonSimple") << " overlaps with tau, dR = " << dR;
184  muonsToCheck.push_back(muon.get());
185  }
186  }
187  }
188  //now examine selected muons
189  bool pass = true;
190  std::vector<int> numMatchesDT(4);
191  std::vector<int> numMatchesCSC(4);
192  std::vector<int> numMatchesRPC(4);
193  std::vector<int> numHitsDT(4);
194  std::vector<int> numHitsCSC(4);
195  std::vector<int> numHitsRPC(4);
196  //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.
197  for (int iStation = 0; iStation < 4; ++iStation) {
198  numMatchesDT[iStation] = 0;
199  numMatchesCSC[iStation] = 0;
200  numMatchesRPC[iStation] = 0;
201  numHitsDT[iStation] = 0;
202  numHitsCSC[iStation] = 0;
203  numHitsRPC[iStation] = 0;
204  }
205  int numSTAMuons = 0, numRPCMuons = 0;
206  if (verbosity_ && !muonsToCheck.empty())
207  edm::LogPrint("PFTauAgainstMuonSimple") << "Muons to check (" << muonsToCheck.size() << "):";
208  size_t iMu = 0;
209  for (const auto& mu : muonsToCheck) {
210  if (mu->isStandAloneMuon())
211  numSTAMuons++;
212  if (mu->muonID("RPCMuLoose"))
213  numRPCMuons++;
214  reco::tau::countMatches(*mu, numMatchesDT, numMatchesCSC, numMatchesRPC);
215  int numStationsWithMatches = 0;
216  for (int iStation = 0; iStation < 4; ++iStation) {
217  if (numMatchesDT[iStation] > 0 && !maskMatchesDT_[iStation])
218  ++numStationsWithMatches;
219  if (numMatchesCSC[iStation] > 0 && !maskMatchesCSC_[iStation])
220  ++numStationsWithMatches;
221  if (numMatchesRPC[iStation] > 0 && !maskMatchesRPC_[iStation])
222  ++numStationsWithMatches;
223  }
224  reco::tau::countHits(*mu, numHitsDT, numHitsCSC, numHitsRPC);
225  int numLast2StationsWithHits = 0;
226  for (int iStation = 2; iStation < 4; ++iStation) {
227  if (numHitsDT[iStation] > 0 && !maskHitsDT_[iStation])
228  ++numLast2StationsWithHits;
229  if (numHitsCSC[iStation] > 0 && !maskHitsCSC_[iStation])
230  ++numLast2StationsWithHits;
231  if (numHitsRPC[iStation] > 0 && !maskHitsRPC_[iStation])
232  ++numLast2StationsWithHits;
233  }
234  if (verbosity_)
235  edm::LogPrint("PFTauAgainstMuonSimple")
236  << "\t" << iMu << ": Pt = " << mu->pt() << ", eta = " << mu->eta() << ", phi = " << mu->phi() << "\n\t"
237  << " isSTA: " << mu->isStandAloneMuon() << ", isRPCLoose: " << mu->muonID("RPCMuLoose")
238  << "\n\t numMatchesDT = " << format_vint(numMatchesDT)
239  << "\n\t numMatchesCSC = " << format_vint(numMatchesCSC)
240  << "\n\t numMatchesRPC = " << format_vint(numMatchesRPC)
241  << "\n\t --> numStationsWithMatches = " << numStationsWithMatches
242  << "\n\t numHitsDT = " << format_vint(numHitsDT) << "\n\t numHitsCSC = " << format_vint(numHitsCSC)
243  << "\n\t numHitsRPC = " << format_vint(numHitsRPC)
244  << "\n\t --> numLast2StationsWithHits = " << numLast2StationsWithHits;
245  if (maxNumberOfMatches_ >= 0 && numStationsWithMatches > maxNumberOfMatches_) {
246  pass = false;
247  break;
248  }
249  if (maxNumberOfHitsLast2Stations_ >= 0 && numLast2StationsWithHits > maxNumberOfHitsLast2Stations_) {
250  pass = false;
251  break;
252  }
253  if (maxNumberOfSTAMuons_ >= 0 && numSTAMuons > maxNumberOfSTAMuons_) {
254  pass = false;
255  break;
256  }
257  if (maxNumberOfRPCMuons_ >= 0 && numRPCMuons > maxNumberOfRPCMuons_) {
258  pass = false;
259  break;
260  }
261  iMu++;
262  }
263 
264  double discriminatorValue = pass ? 1. : 0.;
265  if (verbosity_)
266  edm::LogPrint("PFTauAgainstMuonSimple") << "--> returning discriminatorValue = " << discriminatorValue;
267 
268  return discriminatorValue;
269  }
270 
271 } // namespace
272 
273 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuonSimple);
muonTagProbeFilters_cff.matched
matched
Definition: muonTagProbeFilters_cff.py:62
Handle.h
pat::PackedCandidate::rawCaloFraction
float rawCaloFraction() const
Definition: PackedCandidate.h:920
TauDiscriminationProducerBase.h
vint
double vint[400]
Definition: JetMatchingHook.cc:18
muon
Definition: MuonCocktails.h:17
pat::PackedCandidate::p
double p() const override
magnitude of momentum vector
Definition: PackedCandidate.h:463
pat::PackedCandidate::pdgId
int pdgId() const override
PDG identifier.
Definition: PackedCandidate.h:832
amptDefaultParameters_cff.mu
mu
Definition: amptDefaultParameters_cff.py:16
edm::EDGetTokenT< pat::MuonCollection >
Muon.h
edm::LogInfo
Definition: MessageLogger.h:254
RecoTauMuonTools.h
edm::Ptr::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:139
HLT_2018_cff.muon
muon
Definition: HLT_2018_cff.py:10349
edm::Ptr::key
key_type key() const
Definition: Ptr.h:163
edm::Handle
Definition: AssociativeIterator.h:50
edm::Ref< PFTauCollection >
reco::tau::format_vint
std::string format_vint(const std::vector< int > &vi)
Definition: RecoTauMuonTools.cc:33
deltaR.h
Track.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
reco::tau::countMatches
void countMatches(const reco::Muon &muon, std::vector< int > &numMatchesDT, std::vector< int > &numMatchesCSC, std::vector< int > &numMatchesRPC)
Definition: RecoTauMuonTools.cc:46
TauDiscriminationProducerBase::discriminate
virtual TauDiscriminatorDataType discriminate(const TauRef &tau) const =0
TauDiscriminationProducerBase
Definition: TauDiscriminationProducerBase.h:55
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:528
PbPb_ZMuSkimMuonDPG_cff.deltaR
deltaR
Definition: PbPb_ZMuSkimMuonDPG_cff.py:63
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::ParameterSet
Definition: ParameterSet.h:36
edm::Ptr::id
ProductID id() const
Accessor for product ID.
Definition: Ptr.h:158
TauDiscriminationProducerBase::moduleLabel_
std::string moduleLabel_
Definition: TauDiscriminationProducerBase.h:106
pat::PackedCandidate
Definition: PackedCandidate.h:22
PackedCandidate.h
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
GsfTrack.h
edm::LogPrint
Definition: MessageLogger.h:342
Max
T Max(T a, T b)
Definition: MathUtil.h:44
edm::EventSetup
Definition: EventSetup.h:57
HitPattern.h
TauDiscriminationProducerBase::beginEvent
virtual void beginEvent(const edm::Event &, const edm::EventSetup &)
Definition: TauDiscriminationProducerBase.h:74
edm::Ptr< Candidate >
looper.cfg
cfg
Definition: looper.py:297
std
Definition: JetResolutionObject.h:76
reco::tau::countHits
void countHits(const reco::Muon &muon, std::vector< int > &numHitsDT, std::vector< int > &numHitsCSC, std::vector< int > &numHitsRPC)
Definition: RecoTauMuonTools.cc:7
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
Exception.h
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
edm::Ref::key
key_type key() const
Accessor for product key.
Definition: Ref.h:250
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
Min
T Min(T a, T b)
Definition: MathUtil.h:39
edm::Event
Definition: Event.h:73
patZpeak.numMuons
numMuons
Definition: patZpeak.py:41
edm::InputTag
Definition: InputTag.h:15
Generated for CMSSW Reference Manual by doxygen 1.8.16