CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationAgainstMuon2.cc

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; // CV: map into range 0..3
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);