Skimming/src/BTagSkimLeptonJet.cc

0001 /** \class BTaSkimLeptonJet
0002  *
0003  *
0004  *
0005  * \author Francisco Yumiceva, FERMILAB
0006  *
0007  */
0008
0009 #include <cmath>
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014
0015 #include "DataFormats/JetReco/interface/CaloJet.h"
0016 #include "DataFormats/MuonReco/interface/Muon.h"
0017 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0018
0019 //#include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectron.h"
0020 //#include "DataFormats/EgammaCandidates/interface/PixelMatchGsfElectronFwd.h"
0021 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0022 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0023
0024 #include "RecoBTag/Skimming/interface/BTagSkimLeptonJet.h"
0025
0026 #include "TVector3.h"
0027 #include "Math/GenVector/VectorUtil.h"
0028 #include "Math/GenVector/PxPyPzE4D.h"
0029
0030 using namespace edm;
0031 using namespace std;
0032 using namespace reco;
0033
0034 BTagSkimLeptonJet::BTagSkimLeptonJet(const edm::ParameterSet& iConfig) : nEvents_(0), nAccepted_(0) {
0035   CaloJetInput_ = iConfig.getParameter<InputTag>("CaloJet");
0036   MinCaloJetPt_ = iConfig.getParameter<double>("MinimumCaloJetPt");
0037   MaxCaloJetEta_ = iConfig.getParameter<double>("MaximumCaloJetEta");
0038
0039   LeptonType_ = iConfig.getParameter<std::string>("LeptonType");
0040   if (LeptonType_ != "electron" && LeptonType_ != "muon")
0041     edm::LogError("BTagSkimLeptonJet") << "unknown lepton type !!";
0042   LeptonInput_ = iConfig.getParameter<InputTag>("Lepton");
0043   MinLeptonPt_ = iConfig.getParameter<double>("MinimumLeptonPt");
0044   MaxLeptonEta_ = iConfig.getParameter<double>("MaximumLeptonEta");
0045   //MinNLepton_ = iConfig.getParameter<int>( "MinimumNLepton" );
0046
0047   MaxDeltaR_ = iConfig.getParameter<double>("MaximumDeltaR");
0048
0049   MinPtRel_ = iConfig.getParameter<double>("MinimumPtRel");
0050
0051   MinNLeptonJet_ = iConfig.getParameter<int>("MinimumNLeptonJet");
0052   if (MinNLeptonJet_ < 1)
0053     edm::LogError("BTagSkimLeptonJet") << "MinimunNCaloJet < 1 !!";
0054 }
0055
0056 /*------------------------------------------------------------------------*/
0057
0058 BTagSkimLeptonJet::~BTagSkimLeptonJet() {}
0059
0060 /*------------------------------------------------------------------------*/
0061
0062 bool BTagSkimLeptonJet::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0063   nEvents_++;
0064
0065   Handle<CaloJetCollection> CaloJetsHandle;
0066   //try {
0067   iEvent.getByLabel(CaloJetInput_, CaloJetsHandle);
0068   //}
0069   //catch ( cms::Exception& ex ) {
0070   //edm::LogError( "BTagSkimLeptonJet" )
0071   //  << "Unable to get CaloJet collection "
0072   //  << CaloJetInput_.label();
0073   //return false;
0074   //}
0075   if (CaloJetsHandle->empty())
0076     return false;
0077   CaloJetCollection TheCaloJets = *CaloJetsHandle;
0078   std::stable_sort(TheCaloJets.begin(), TheCaloJets.end(), PtSorter());
0079
0080   //int nLepton = 0;
0081
0082   Handle<MuonCollection> MuonHandle;
0083   MuonCollection TheMuons;
0084   if (LeptonType_ == "muon") {
0085     //try{
0086     iEvent.getByLabel(LeptonInput_, MuonHandle);
0087     //  }
0088     //    catch ( cms::Exception& ex ) {
0089     //edm::LogError( "BTagSkimLeptonJet" )
0090     //    << "Unable to get muon collection "
0091     //    << LeptonInput_.label();
0092     //return false;
0093     //}
0094     TheMuons = *MuonHandle;
0095     std::stable_sort(TheMuons.begin(), TheMuons.end(), PtSorter());
0096   }
0097
0098   Handle<GsfElectronCollection> ElectronHandle;
0099   GsfElectronCollection TheElectrons;
0100   if (LeptonType_ == "electron") {
0101     //try{
0102     iEvent.getByLabel(LeptonInput_, ElectronHandle);
0103     //  }
0104     // catch ( cms::Exception& ex ) {
0105     // edm::LogError( "BTagSkimLeptonJet" )
0106     //        << "Unable to get electron collection "
0107     //        << LeptonInput_.label();
0108     // return false;
0109     // }
0110     TheElectrons = *ElectronHandle;
0111     std::stable_sort(TheElectrons.begin(), TheElectrons.end(), PtSorter());
0112   }
0113
0114   // Jet cuts
0115   int nJet = 0;
0116   for (CaloJetCollection::const_iterator ajet = TheCaloJets.begin(); ajet != TheCaloJets.end(); ++ajet) {
0117     if ((fabs(ajet->eta()) < MaxCaloJetEta_) && (ajet->pt() > MinCaloJetPt_)) {
0118       if (LeptonType_ == "muon") {
0119         for (MuonCollection::const_iterator amuon = TheMuons.begin(); amuon != TheMuons.end(); ++amuon) {
0120           // select good muon
0121           if ((amuon->pt() > MinLeptonPt_) && (fabs(amuon->eta()) < MaxLeptonEta_)) {
0122             double deltar = ROOT::Math::VectorUtil::DeltaR(ajet->p4().Vect(), amuon->momentum());
0123
0124             TVector3 jetvec(ajet->p4().Vect().X(), ajet->p4().Vect().Y(), ajet->p4().Vect().Z());
0125             TVector3 muvec(amuon->momentum().X(), amuon->momentum().Y(), amuon->momentum().Z());
0126             jetvec += muvec;
0127             double ptrel = muvec.Perp(jetvec);
0128
0129             if ((deltar < MaxDeltaR_) && (ptrel > MinPtRel_))
0130               nJet++;
0131             if (nJet >= MinNLeptonJet_)
0132               break;
0133           }
0134         }
0135       }
0136
0137       if (LeptonType_ == "electron") {
0138         for (GsfElectronCollection::const_iterator anelectron = TheElectrons.begin(); anelectron != TheElectrons.end();
0139              anelectron++) {
0140           if ((anelectron->pt() > MinLeptonPt_) && (fabs(anelectron->eta()) < MaxLeptonEta_)) {
0141             double deltar = ROOT::Math::VectorUtil::DeltaR(ajet->p4().Vect(), anelectron->momentum());
0142
0143             TVector3 jetvec(ajet->p4().Vect().X(), ajet->p4().Vect().Y(), ajet->p4().Vect().Z());
0144             TVector3 evec(anelectron->momentum().X(), anelectron->momentum().Y(), anelectron->momentum().Z());
0145             jetvec += evec;
0146             double ptrel = evec.Perp(jetvec);
0147
0148             if ((deltar < MaxDeltaR_) && (ptrel > MinPtRel_))
0149               nJet++;
0150             if (nJet >= MinNLeptonJet_)
0151               break;
0152           }
0153         }
0154       }
0155
0156     }  // close jet selection
0157   }
0158
0159   if (nJet < MinNLeptonJet_)
0160     return false;
0161
0162   nAccepted_++;
0163
0164   return true;
0165 }
0166
0167 /*------------------------------------------------------------------------*/
0168
0169 void BTagSkimLeptonJet::endJob() {
0170   edm::LogVerbatim("BTagSkimLeptonJet")
0171       << "=============================================================================\n"
0172       << " Events read: " << nEvents_ << "\n Events accepted by (" << LeptonType_
0173       << ") BTagSkimLeptonJet: " << nAccepted_ << "\n Efficiency: " << (double)(nAccepted_) / (double)(nEvents_)
0174       << "\n===========================================================================" << endl;
0175 }
0176
0177 //define this as a plug-in
0178 DEFINE_FWK_MODULE(BTagSkimLeptonJet);