Go to the documentation of this file.00001
00014 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
00015
00016 #include "FWCore/Utilities/interface/Exception.h"
00017
00018 #include "DataFormats/MuonReco/interface/Muon.h"
00019 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00020 #include "DataFormats/Common/interface/Handle.h"
00021 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00022 #include "DataFormats/TrackReco/interface/HitPattern.h"
00023 #include "DataFormats/TrackReco/interface/Track.h"
00024 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00025 #include "DataFormats/Math/interface/deltaR.h"
00026
00027 #include <vector>
00028 #include <string>
00029 #include <iostream>
00030
00031 class PFRecoTauDiscriminationAgainstMuon2 : public PFTauDiscriminationProducerBase
00032 {
00033 enum { kLoose, kMedium, kTight, kCustom };
00034 public:
00035 explicit PFRecoTauDiscriminationAgainstMuon2(const edm::ParameterSet& cfg)
00036 : PFTauDiscriminationProducerBase(cfg),
00037 moduleLabel_(cfg.getParameter<std::string>("@module_label"))
00038 {
00039 std::string discriminatorOption_string = cfg.getParameter<std::string>("discriminatorOption");
00040 if ( discriminatorOption_string == "loose" ) discriminatorOption_ = kLoose;
00041 else if ( discriminatorOption_string == "medium" ) discriminatorOption_ = kMedium;
00042 else if ( discriminatorOption_string == "tight" ) discriminatorOption_ = kTight;
00043 else if ( discriminatorOption_string == "custom" ) discriminatorOption_ = kCustom;
00044 else throw edm::Exception(edm::errors::UnimplementedFeature)
00045 << " Invalid Configuration parameter 'discriminatorOption' = " << discriminatorOption_string << " !!\n";
00046 hop_ = cfg.getParameter<double>("HoPMin");
00047 maxNumberOfMatches_ = cfg.exists("maxNumberOfMatches") ? cfg.getParameter<int>("maxNumberOfMatches"): 0;
00048 doCaloMuonVeto_ = cfg.exists("doCaloMuonVeto") ? cfg.getParameter<bool>("doCaloMuonVeto"): false;
00049 maxNumberOfHitsLast2Stations_ = cfg.exists("maxNumberOfHitsLast2Stations") ? cfg.getParameter<int>("maxNumberOfHitsLast2Stations"): 0;
00050 if ( cfg.exists("srcMuons") ) {
00051 srcMuons_ = cfg.getParameter<edm::InputTag>("srcMuons");
00052 dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
00053 }
00054 verbosity_ = cfg.exists("verbosity") ? cfg.getParameter<int>("verbosity") : 0;
00055 }
00056 ~PFRecoTauDiscriminationAgainstMuon2() {}
00057
00058 void beginEvent(const edm::Event&, const edm::EventSetup&);
00059
00060 double discriminate(const reco::PFTauRef&);
00061
00062 private:
00063 std::string moduleLabel_;
00064 int discriminatorOption_;
00065 double hop_;
00066 int maxNumberOfMatches_;
00067 bool doCaloMuonVeto_;
00068 int maxNumberOfHitsLast2Stations_;
00069 edm::InputTag srcMuons_;
00070 edm::Handle<reco::MuonCollection> muons_;
00071 double dRmuonMatch_;
00072 int verbosity_;
00073 };
00074
00075 void PFRecoTauDiscriminationAgainstMuon2::beginEvent(const edm::Event& evt, const edm::EventSetup& es)
00076 {
00077 if ( srcMuons_.label() != "" ) {
00078 evt.getByLabel(srcMuons_, muons_);
00079 }
00080 }
00081
00082 namespace
00083 {
00084 void countHits(const reco::Muon& muon, std::vector<int>& numHitsDT, std::vector<int>& numHitsCSC, std::vector<int>& numHitsRPC)
00085 {
00086 if ( muon.outerTrack().isNonnull() ) {
00087 const reco::HitPattern& muonHitPattern = muon.outerTrack()->hitPattern();
00088 for ( int iHit = 0; iHit < muonHitPattern.numberOfHits(); ++iHit ) {
00089 uint32_t hit = muonHitPattern.getHitPattern(iHit);
00090 if ( hit == 0 ) break;
00091 if ( muonHitPattern.muonHitFilter(hit) && (muonHitPattern.getHitType(hit) == TrackingRecHit::valid || muonHitPattern.getHitType(hit) == TrackingRecHit::bad) ) {
00092 int muonStation = muonHitPattern.getMuonStation(hit) - 1;
00093 if ( muonStation >= 0 && muonStation < 4 ) {
00094 if ( muonHitPattern.muonDTHitFilter(hit) ) ++numHitsDT[muonStation];
00095 else if ( muonHitPattern.muonCSCHitFilter(hit) ) ++numHitsCSC[muonStation];
00096 else if ( muonHitPattern.muonRPCHitFilter(hit) ) ++numHitsRPC[muonStation];
00097 }
00098 }
00099 }
00100 }
00101 }
00102
00103 std::string format_vint(const std::vector<int>& vi)
00104 {
00105 std::ostringstream os;
00106 os << "{ ";
00107 unsigned numEntries = vi.size();
00108 for ( unsigned iEntry = 0; iEntry < numEntries; ++iEntry ) {
00109 os << vi[iEntry];
00110 if ( iEntry < (numEntries - 1) ) os << ", ";
00111 }
00112 os << " }";
00113 return os.str();
00114 }
00115 }
00116
00117 double PFRecoTauDiscriminationAgainstMuon2::discriminate(const reco::PFTauRef& pfTau)
00118 {
00119 if ( verbosity_ ) {
00120 std::cout << "<PFRecoTauDiscriminationAgainstMuon2::discriminate>:" << std::endl;
00121 std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
00122 std::cout << "tau #" << pfTau.key() << ": Pt = " << pfTau->pt() << ", eta = " << pfTau->eta() << ", phi = " << pfTau->phi() << std::endl;
00123 }
00124
00125 int numMatches = 0;
00126
00127 std::vector<int> numHitsDT(4);
00128 std::vector<int> numHitsCSC(4);
00129 std::vector<int> numHitsRPC(4);
00130 for ( int iStation = 0; iStation < 4; ++iStation ) {
00131 numHitsDT[iStation] = 0;
00132 numHitsCSC[iStation] = 0;
00133 numHitsRPC[iStation] = 0;
00134 }
00135
00136 const reco::PFCandidateRef& pfLeadChargedHadron = pfTau->leadPFChargedHadrCand();
00137 if ( pfLeadChargedHadron.isNonnull() ) {
00138 reco::MuonRef muonRef = pfLeadChargedHadron->muonRef();
00139 if ( muonRef.isNonnull() ) {
00140 if ( verbosity_ ) std::cout << " has muonRef." << std::endl;
00141 numMatches = muonRef->numberOfMatches(reco::Muon::NoArbitration);
00142 countHits(*muonRef, numHitsDT, numHitsCSC, numHitsRPC);
00143 }
00144 }
00145
00146 if ( srcMuons_.label() != "" ) {
00147 size_t numMuons = muons_->size();
00148 for ( size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon ) {
00149 reco::MuonRef muon(muons_, idxMuon);
00150 if ( verbosity_ ) std::cout << "muon #" << muon.key() << ": Pt = " << muon->pt() << ", eta = " << muon->eta() << ", phi = " << muon->phi() << std::endl;
00151 if ( pfLeadChargedHadron.isNonnull() && pfLeadChargedHadron->muonRef().isNonnull() && muon == pfLeadChargedHadron->muonRef() ) {
00152 if ( verbosity_ ) std::cout << " matches muonRef of tau --> skipping it." << std::endl;
00153 continue;
00154 }
00155 double dR = deltaR(muon->p4(), pfTau->p4());
00156 if ( dR < dRmuonMatch_ ) {
00157 if ( verbosity_ ) std::cout << " overlaps with tau, dR = " << dR << std::endl;
00158 numMatches += muon->numberOfMatches(reco::Muon::NoArbitration);
00159 countHits(*muon, numHitsDT, numHitsCSC, numHitsRPC);
00160 }
00161 }
00162 }
00163
00164 if ( verbosity_ ) {
00165 std::cout << "numMatches = " << numMatches << std::endl;
00166 std::cout << "numHitsDT = " << format_vint(numHitsDT) << std::endl;
00167 std::cout << "numHitsCSC = " << format_vint(numHitsCSC) << std::endl;
00168 std::cout << "numHitsRPC = " << format_vint(numHitsRPC) << std::endl;
00169 }
00170
00171 int numLast2StationsWithHits = 0;
00172 for ( int iStation = 2; iStation < 4; ++iStation ) {
00173 if ( numHitsDT[iStation] > 0 ) ++numLast2StationsWithHits;
00174 if ( numHitsCSC[iStation] > 0 ) ++numLast2StationsWithHits;
00175 if ( numHitsRPC[iStation] > 0 ) ++numLast2StationsWithHits;
00176 }
00177
00178
00179 bool passesCaloMuonVeto = true;
00180 if ( pfLeadChargedHadron.isNonnull() ) {
00181 double energyECALplusHCAL = pfLeadChargedHadron->ecalEnergy() + pfLeadChargedHadron->hcalEnergy();
00182 if ( verbosity_ ) {
00183 if ( pfLeadChargedHadron->trackRef().isNonnull() ) std::cout << "decayMode = " << pfTau->decayMode() << ", energy(ECAL+HCAL) = " << energyECALplusHCAL << ", leadPFChargedHadronP = " << pfLeadChargedHadron->trackRef()->p() << std::endl;
00184 else if ( pfLeadChargedHadron->gsfTrackRef().isNonnull() )
00185 std::cout << "decayMode = " << pfTau->decayMode() << ", energy(ECAL+HCAL) = " << energyECALplusHCAL << ", leadPFChargedHadronP = " << pfLeadChargedHadron->gsfTrackRef()->p() << std::endl;
00186
00187 }
00188 const reco::Track* leadTrack = 0;
00189 if ( pfLeadChargedHadron->trackRef().isNonnull() ) leadTrack = pfLeadChargedHadron->trackRef().get();
00190 else if ( pfLeadChargedHadron->gsfTrackRef().isNonnull() ) leadTrack = pfLeadChargedHadron->gsfTrackRef().get();
00191 if ( pfTau->decayMode() == 0 && leadTrack && energyECALplusHCAL < (hop_*leadTrack->p()) ) passesCaloMuonVeto = false;
00192 }
00193
00194 double discriminatorValue = 0.;
00195 if ( discriminatorOption_ == kLoose && numMatches <= maxNumberOfMatches_ ) discriminatorValue = 1.;
00196 else if ( discriminatorOption_ == kMedium && numMatches <= maxNumberOfMatches_ && numLast2StationsWithHits <= maxNumberOfHitsLast2Stations_ ) discriminatorValue = 1.;
00197 else if ( discriminatorOption_ == kTight && numMatches <= maxNumberOfMatches_ && numLast2StationsWithHits <= maxNumberOfHitsLast2Stations_ && passesCaloMuonVeto ) discriminatorValue = 1.;
00198 else if ( discriminatorOption_ == kCustom ) {
00199 bool pass = true;
00200 if ( maxNumberOfMatches_ >= 0 && numMatches > maxNumberOfMatches_ ) pass = false;
00201 if ( maxNumberOfHitsLast2Stations_ >= 0 && numLast2StationsWithHits > maxNumberOfHitsLast2Stations_ ) pass = false;
00202 if ( doCaloMuonVeto_ && !passesCaloMuonVeto ) pass = false;
00203 discriminatorValue = pass ? 1.: 0.;
00204 }
00205 if ( verbosity_ ) std::cout << "--> returning discriminatorValue = " << discriminatorValue << std::endl;
00206
00207 return discriminatorValue;
00208 }
00209
00210 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuon2);