diff --git a/interface/PatUtilities.h b/interface/PatUtilities.h index 6fbb906d..4fa4030f 100644 --- a/interface/PatUtilities.h +++ b/interface/PatUtilities.h @@ -6,6 +6,7 @@ #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h" #include "DataFormats/PatCandidates/interface/Electron.h" #include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/MET.h" #include "DataFormats/PatCandidates/interface/Jet.h" #include "JetMETCorrections/Objects/interface/JetCorrector.h" #include "FWCore/Framework/interface/Event.h" @@ -30,4 +31,63 @@ pat::JetCollection applyNewJec( pat::JetCollection jets, const JetCorrector* cor float getJECForJet(const pat::Jet jet, const JetCorrector* corrector, edm::Event& iEvent, const edm::EventSetup& iSetup ); +namespace ntp { +namespace utils { +/** + * @param jets + * @return pT sum of all jets + */ +double ht(const pat::JetCollection& jets); +/** + * + * @param jets + * @param met + * @param lepton + * @return ht + MET.et() + lepton.pt() + */ +double st(const pat::JetCollection& jets, const pat::MET& met, const reco::Candidate* lepton); +/** + * + * @param met + * @param lepton + * @return transverse momentum of the leptonically decaying W boson + */ +double wpt(const pat::MET& met, const reco::Candidate* lepton); +/** + * + * @param met + * @param lepton + * @return transverse mass of the leptonically decaying W boson + */ +double mt(const pat::MET& met, const reco::Candidate* lepton); +/** + * Calculates the invariant mass of three jets where the triplet of jets is + * chosen by having the largest (vector) sum of pT + * @param jets + * @return invariant mass of three jets + */ +double m3(const pat::JetCollection& jets); +/** + * @param jets + * @param lepton + * @return angle between the lepton and the closest b-tagged jet + */ +double angle_bl(const pat::JetCollection& jets, const reco::Candidate* lepton); +/** + * @param jets + * @param lepton + * @return invariant mass of the lepton and the closest b-tagged jet + */ +double m_bl(const pat::JetCollection& jets, const reco::Candidate* lepton); + +/** + * + * @param jets + * @param lepton + * @return closest jet to the given lepton + */ +pat::Jet getClosestJet(const pat::JetCollection& jets, const reco::Candidate* lepton); +} +} + #endif diff --git a/plugins/BuildFile.xml b/plugins/BuildFile.xml index b3924b97..8cf4fc8c 100644 --- a/plugins/BuildFile.xml +++ b/plugins/BuildFile.xml @@ -14,4 +14,5 @@ + diff --git a/plugins/UnfoldingAnalyser.cc b/plugins/UnfoldingAnalyser.cc index fc07edca..70acbaf6 100644 --- a/plugins/UnfoldingAnalyser.cc +++ b/plugins/UnfoldingAnalyser.cc @@ -1,3 +1,7 @@ +/** + * @deprecated + */ + #include "BristolAnalysis/NTupleTools/plugins/UnfoldingAnalyser.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" diff --git a/plugins/UnfoldingProducer.cc b/plugins/UnfoldingProducer.cc index 5fe33e58..d1a94372 100644 --- a/plugins/UnfoldingProducer.cc +++ b/plugins/UnfoldingProducer.cc @@ -1,3 +1,7 @@ +/** + * + * @deprecated + */ #include "BristolAnalysis/NTupleTools/plugins/UnfoldingProducer.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" diff --git a/plugins/userdata/EventUserData.cc b/plugins/userdata/EventUserData.cc new file mode 100644 index 00000000..e6b39161 --- /dev/null +++ b/plugins/userdata/EventUserData.cc @@ -0,0 +1,182 @@ +/** + * Produces global variables and similar: + * - HT: scalar sum of pT of all jets above a threshold + * - ST: HT + MET + signal lepton pT + * - + * + * Why do we not include this in BristolNTuple_Event.cc? + * BristolNTuple_Event is meant to run once per event, + * while this module can run N times (once for each selection, jet energy variation, etc). + */ +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "FWCore/Utilities/interface/StreamID.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/MET.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +#include "BristolAnalysis/NTupleTools/interface/PatUtilities.h" + +namespace ntp { +namespace userdata { + +class EventUserData: public edm::stream::EDProducer<> { +public: + typedef std::vector IndexCollection; + explicit EventUserData(const edm::ParameterSet&); + ~EventUserData(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + virtual void beginStream(edm::StreamID) override; + virtual void produce(edm::Event&, const edm::EventSetup&) override; + virtual void endStream() override; + + pat::JetCollection getCleanedJets(const pat::JetCollection& jets, const IndexCollection& indices); + // inputs + const edm::EDGetTokenT > vtxInputTag_; + edm::EDGetToken electronInputTag_, muonInputTag_, jetInputTag_, bJetInputTag_, metInputTag_; + bool isElectron_; + + std::string prefix_; + std::string suffix_; + +}; + +EventUserData::EventUserData(const edm::ParameterSet& iConfig) : + vtxInputTag_( + consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("vertexCollection"))), + electronInputTag_( + consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("electronCollection"))), + muonInputTag_(consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("muonCollection"))), + jetInputTag_(consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("jetCollection"))), // + bJetInputTag_(consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("bJetCollection"))), // + metInputTag_(consumes < std::vector > (iConfig.getParameter < edm::InputTag > ("metInputTag"))), // + isElectron_(iConfig.getParameter("isElectron")), // + prefix_(iConfig.getParameter < std::string > ("prefix")), // + suffix_(iConfig.getParameter < std::string > ("suffix")) +{ + //register your products + produces(prefix_ + "HT" + suffix_); + produces(prefix_ + "ST" + suffix_); + produces(prefix_ + "M3" + suffix_); + produces(prefix_ + "Mbl" + suffix_); + produces(prefix_ + "anglebl" + suffix_); + produces(prefix_ + "MT" + suffix_); + produces(prefix_ + "WPT" + suffix_); +} + +EventUserData::~EventUserData() +{ + + // do anything here that needs to be done at destruction time + // (e.g. close files, deallocate resources etc.) + +} + +// +// member functions +// + +// ------------ method called to produce the data ------------ +void EventUserData::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::Handle < std::vector > electrons; + iEvent.getByToken(electronInputTag_, electrons); + + edm::Handle < std::vector > muons; + iEvent.getByToken(muonInputTag_, muons); + + edm::Handle < std::vector > jetHandle; + iEvent.getByToken(jetInputTag_, jetHandle); + const pat::JetCollection& jets = *jetHandle.product(); + + edm::Handle < std::vector > bJetHandle; + iEvent.getByToken(bJetInputTag_, bJetHandle); + const pat::JetCollection& bjets = *bJetHandle.product(); + + edm::Handle < std::vector > mets; + iEvent.getByToken(metInputTag_, mets); + + const pat::MET & met(mets->at(0)); + + const double DEFAULT_VALUE=-999; + std::auto_ptr ht(new double(utils::ht(jets))); + std::auto_ptr st(new double(DEFAULT_VALUE)); + std::auto_ptr m3(new double(DEFAULT_VALUE)); + std::auto_ptr m_bl(new double(DEFAULT_VALUE)); + std::auto_ptr angle_bl(new double(-9)); + std::auto_ptr mt(new double(DEFAULT_VALUE)); + std::auto_ptr wpt(new double(DEFAULT_VALUE)); + + if (jets.size() >= 3) + *m3 = utils::m3(jets); + + if ((electrons.isValid() && isElectron_) || (muons.isValid() && !isElectron_)) { + const reco::Candidate* lepton = + isElectron_ ? + (reco::Candidate*) &electrons->front() : + (reco::Candidate*) &muons->front(); + if (lepton) { + if (bjets.size() > 0) { + *m_bl = utils::m_bl(bjets, lepton); + *angle_bl = utils::angle_bl(bjets, lepton); + } + *st = utils::st(jets, met, lepton); + *mt = utils::mt(met, lepton); + *wpt = utils::wpt(met, lepton); + } + + } + + iEvent.put(ht, prefix_ + "HT" + suffix_); + iEvent.put(st, prefix_ + "ST" + suffix_); + iEvent.put(m3, prefix_ + "M3" + suffix_); + iEvent.put(m_bl, prefix_ + "Mbl" + suffix_); + iEvent.put(angle_bl, prefix_ + "anglebl" + suffix_); + iEvent.put(mt, prefix_ + "MT" + suffix_); + iEvent.put(wpt, prefix_ + "WPT" + suffix_); + +} + +// ------------ method called once each stream before processing any runs, lumis or events ------------ +void EventUserData::beginStream(edm::StreamID) +{ +} + +// ------------ method called once each stream after processing all runs, lumis and events ------------ +void EventUserData::endStream() +{ +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void EventUserData::fillDescriptions(edm::ConfigurationDescriptions& descriptions) +{ + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE (EventUserData); +} //userdata +} //ntp diff --git a/python/run/miniAODToNTuple_cfg.py b/python/run/miniAODToNTuple_cfg.py index 1f712d57..beeb2092 100644 --- a/python/run/miniAODToNTuple_cfg.py +++ b/python/run/miniAODToNTuple_cfg.py @@ -129,6 +129,9 @@ # mapping between MiniAOD collections and our object selections process.load('BristolAnalysis.NTupleTools.indices_cff') +# adds process.eventUserDataSequence +process.load('BristolAnalysis.NTupleTools.userdata.EventUserData_cff') + if isTTbarMC: process.makingNTuples = cms.Path( @@ -143,6 +146,7 @@ process.selectionCriteriaAnalyzer * process.makePseudoTop * process.indexSequence * + process.eventUserDataSequence * process.printEventContent * process.nTuples * process.nTupleTree @@ -158,6 +162,7 @@ process.qcdElectronSelectionAnalyzerSequence * process.selectionCriteriaAnalyzer * process.indexSequence * + process.eventUserDataSequence * process.printEventContent * process.nTuples * process.nTupleTree @@ -178,6 +183,7 @@ 'keep bool_topPairEPlusJetsQCDSelectionTagging_*FullSelection*_*', 'keep bool_topPairEPlusJetsConversionSelectionTagging_*FullSelection*_*', 'keep uint*_*Indices*_*_*', + 'keep double_eventUserData*_*_*', ] ) @@ -311,6 +317,7 @@ process.load('BristolAnalysis.NTupleTools.userdata.ElectronUserData_cfi') process.load('BristolAnalysis.NTupleTools.userdata.MuonUserData_cfi') process.load('BristolAnalysis.NTupleTools.userdata.JetUserData_cfi') + if options.useJECFromFile: process.jetUserData.jetCollection=cms.InputTag("patJetsReapplyJEC") diff --git a/python/userdata/EventUserData_cff.py b/python/userdata/EventUserData_cff.py new file mode 100644 index 00000000..5683f0ab --- /dev/null +++ b/python/userdata/EventUserData_cff.py @@ -0,0 +1,51 @@ +import FWCore.ParameterSet.Config as cms + +from BristolAnalysis.NTupleTools.userdata.EventUserData_cfi import eventUserData + +# electron signal selection +eventUserDataTopPairElectronPlusJetsSelection = eventUserData.clone() +# electron main QCD control region +eventUserDataTopPairElectronPlusJetsConversionSelection = eventUserData.clone( + prefix='TopPairElectronPlusJetsConversionSelection.', + electronCollection='goodConversionElectrons', + jetCollection='goodJetsEConversionRegion', + bJetCollection='goodBJetsEConversionRegion', +) +# electron secondary QCD control region +eventUserDataTopPairElectronPlusJetsNonIsoSelection = eventUserData.clone( + prefix='TopPairElectronPlusJetsQCDSelection.', + electronCollection='goodNonIsoElectrons', + jetCollection='goodJetsENonIsoRegion', + bJetCollection='goodBJetsENonIsoRegion', +) + +# muon signal selection +eventUserDataTopPairMuonPlusJetsSelection = eventUserData.clone( + prefix='TopPairMuonPlusJetsSelection.', + isElectron=False +) +# muon main QCD control region +eventUserDataTopPairMuonPlusJetsQCD1Selection = eventUserData.clone( + prefix='TopPairMuonPlusJetsQCDSelection3toInf.', + muonCollection='goodNonIsoR1Muons', + jetCollection='goodJetsMuNonIsoR1Region', + bJetCollection='goodBJetsMuNonIsoR1Region', +) +# muon secondary QCD control region +eventUserDataTopPairMuonPlusJetsQCD2Selection = eventUserData.clone( + prefix='TopPairMuonPlusJetsQCDSelection1p5to3.', + muonCollection='goodNonIsoR2Muons', + jetCollection='goodJetsMuNonIsoR2Region', + bJetCollection='goodBJetsMuNonIsoR2Region', +) + +eventUserDataSequence = cms.Sequence( + # electrons + eventUserDataTopPairElectronPlusJetsSelection + + eventUserDataTopPairElectronPlusJetsConversionSelection + + eventUserDataTopPairElectronPlusJetsNonIsoSelection + + # muons + eventUserDataTopPairMuonPlusJetsSelection + + eventUserDataTopPairMuonPlusJetsQCD1Selection + + eventUserDataTopPairMuonPlusJetsQCD2Selection +) diff --git a/python/userdata/EventUserData_cfi.py b/python/userdata/EventUserData_cfi.py new file mode 100644 index 00000000..3f0fb977 --- /dev/null +++ b/python/userdata/EventUserData_cfi.py @@ -0,0 +1,17 @@ +import FWCore.ParameterSet.Config as cms + +eventUserData = cms.EDProducer( + 'EventUserData', + vertexCollection=cms.InputTag('offlineSlimmedPrimaryVertices'), + electronCollection=cms.InputTag("goodElectrons"), + muonCollection=cms.InputTag("goodMuons"), + jetCollection=cms.InputTag("goodJets"), + bJetCollection=cms.InputTag("goodBJets"), + metInputTag=cms.InputTag('slimmedMETs'), + + # selection and similar + isElectron=cms.bool(True), + # + prefix=cms.string('TopPairElectronPlusJetsSelection.'), + suffix=cms.string(''), +) diff --git a/src/BristolNTuple_GlobalEventVars.cc b/src/BristolNTuple_GlobalEventVars.cc index abc8f540..183b51ed 100644 --- a/src/BristolNTuple_GlobalEventVars.cc +++ b/src/BristolNTuple_GlobalEventVars.cc @@ -1,3 +1,8 @@ +/** + * + * @deprecated + */ + #include "BristolAnalysis/NTupleTools/interface/BristolNTuple_GlobalEventVars.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" diff --git a/src/PatUtilities.cc b/src/PatUtilities.cc index b295c2a7..625e2d90 100644 --- a/src/PatUtilities.cc +++ b/src/PatUtilities.cc @@ -5,6 +5,7 @@ #include "DataFormats/PatCandidates/interface/Isolation.h" #include "EgammaAnalysis/ElectronTools/interface/ElectronEffectiveArea.h" #include "FWCore/Framework/interface/EDConsumerBase.h" +#include "TLorentzVector.h" using namespace edm; using namespace std; @@ -90,58 +91,6 @@ unsigned int findTrigger(const std::string& triggerWildCard, const HLTConfigProv return found; } -//double getRelativeIsolation(const pat::Electron& electron, double cone, double rho, bool isRealData, bool useDeltaBetaCorrections, -// bool useRhoActiveAreaCorrections) { -// //code from: https://twiki.cern.ch/twiki/bin/view/CMS/PfIsolation -// float AEff = 0.00; -// //so far this is only for 0.3. 0.4 is available but not 0.5 -// -//// std::cout << "ElectronEffectiveArea::kEleEAData2011: " << ElectronEffectiveArea::kEleEAData2011 << std::endl; -//// std::cout << "ElectronEffectiveArea::kEleEAData2012: " << ElectronEffectiveArea::kEleEAData2012<< std::endl; -// -//// std::cout << "ElectronEffectiveArea::kEleEAFall1MC: " << ElectronEffectiveArea::kEleEAFall11MC << std::endl; -//// std::cout << "ElectronEffectiveArea::kEleEASummer12MC: " << ElectronEffectiveArea::kEleEASummer12MC<< std::endl; -// -// if (isRealData) { -// AEff = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, -// electron.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2012); -// } else { -// AEff = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, -// electron.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2012); -// } -// -// AbsVetos vetos_ch; -// AbsVetos vetos_nh; -// AbsVetos vetos_ph; -// -// Direction Dir = Direction(electron.superCluster()->eta(), electron.superCluster()->phi()); -// -// //pf isolation veto setup EGM recommendation -// if (fabs(electron.superCluster()->eta()) > 1.479) { -// vetos_ch.push_back(new ConeVeto(Dir, 0.015)); -// vetos_ph.push_back(new ConeVeto(Dir, 0.08)); -// } -// -// const double chIso = electron.isoDeposit(pat::PfChargedHadronIso)->depositAndCountWithin(cone, vetos_ch).first; -// const double nhIso = electron.isoDeposit(pat::PfNeutralHadronIso)->depositAndCountWithin(cone, vetos_nh).first; -// const double phIso = electron.isoDeposit(pat::PfGammaIso)->depositAndCountWithin(cone, vetos_ph).first; -// -// const double puChIso = electron.isoDeposit(pat::PfPUChargedHadronIso)->depositAndCountWithin(cone, vetos_ch).first; -// -// const double relIso = (chIso + nhIso + phIso) / electron.pt(); -// const double relIsodb = (chIso + max(0.0, nhIso + phIso - 0.5 * puChIso)) / electron.pt(); -// const double relIsorho = (chIso + max(0.0, nhIso + phIso - rho * AEff)) / electron.pt(); -// -// if (useRhoActiveAreaCorrections) -// return relIsorho; -// -// if (useDeltaBetaCorrections) -// return relIsodb; -// -// -// return relIso; -//} - double getSmearedJetPtScale(const pat::Jet& jet, int jet_smearing_systematic) { // Get the jet energy resolution scale factors, depending on the jet eta, from // https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetResolution#Recommendations_for_7_and_8_TeV @@ -245,4 +194,110 @@ float getJECForJet(const pat::Jet jet, const JetCorrector* corrector, edm::Event } +namespace ntp { +namespace utils { +double ht(const pat::JetCollection& jets) +{ + double result = 0; + for (auto jet : jets) { + result += jet.pt(); + } + return result; +} + +double st(const pat::JetCollection& jets, const pat::MET& met, const reco::Candidate* lepton) +{ + return ht(jets) + met.pt() + lepton->pt(); +} + +double wpt(const pat::MET& met, const reco::Candidate* lepton) +{ + return sqrt(pow(lepton->px() + met.px(), 2) + pow(lepton->py() + met.py(), 2)); +} + +double mt(const pat::MET& met, const reco::Candidate* lepton) +{ + double energy_squared = pow(lepton->et() + met.pt(), 2); + double momentum_squared = pow(lepton->px() + met.px(), 2) + pow(lepton->py() + met.py(), 2); + double MT_squared = energy_squared - momentum_squared; + + if (MT_squared > 0) + return sqrt(MT_squared); + else + return -1; +} + +double m3(const pat::JetCollection& jets) +{ + double m3(-1), max_pt(0); + if (jets.size() >= 3) { + for (unsigned int index1 = 0; index1 < jets.size() - 2; ++index1) { + for (unsigned int index2 = index1 + 1; index2 < jets.size() - 1; ++index2) { + for (unsigned int index3 = index2 + 1; index3 < jets.size(); ++index3) { + TLorentzVector jet1(jets.at(index1).px(), jets.at(index1).py(), jets.at(index1).pz(), + jets.at(index1).energy()); + TLorentzVector jet2(jets.at(index2).px(), jets.at(index2).py(), jets.at(index2).pz(), + jets.at(index2).energy()); + TLorentzVector jet3(jets.at(index3).px(), jets.at(index3).py(), jets.at(index3).pz(), + jets.at(index3).energy()); + + TLorentzVector m3Vector(jet1 + jet2 + jet3); + double currentPt = m3Vector.Pt(); + if (currentPt > max_pt) { + max_pt = currentPt; + m3 = m3Vector.M(); + } + } + } + } + } + + return m3; +} + +double angle_bl(const pat::JetCollection& jets, const reco::Candidate* lepton) +{ + double angle(-1); + const pat::Jet closestJet(getClosestJet(jets, lepton)); + + TLorentzVector b(closestJet.px(), closestJet.py(), closestJet.pz(), closestJet.energy()); + TLorentzVector l(lepton->px(), lepton->py(), lepton->pz(), lepton->energy()); + + angle = b.Angle(l.Vect()); + + return angle; +} + +pat::Jet getClosestJet(const pat::JetCollection& jets, const reco::Candidate* lepton) +{ + double closestDR = 9999; + pat::Jet closestJet; + + for (auto jet : jets) { + double DR = deltaR(lepton->eta(), lepton->phi(), jet.eta(), jet.phi()); + if (DR < closestDR) { + closestDR = DR; + closestJet = pat::Jet(jet); + } + } + + return closestJet; +} + +double m_bl(const pat::JetCollection& jets, const reco::Candidate* lepton) +{ + double mass(-1); + const pat::Jet closestJet(getClosestJet(jets, lepton)); + + TLorentzVector b(closestJet.px(), closestJet.py(), closestJet.pz(), closestJet.energy()); + TLorentzVector l(lepton->px(), lepton->py(), lepton->pz(), lepton->energy()); + + mass = (b + l).M(); + + return mass; +} +} +} + +