diff --git a/DataFormats/METReco/interface/MET.h b/DataFormats/METReco/interface/MET.h index 642df0e0b0ac1..a1293376fb05f 100755 --- a/DataFormats/METReco/interface/MET.h +++ b/DataFormats/METReco/interface/MET.h @@ -71,6 +71,9 @@ namespace reco //________________________________________________________________________|| void setSignificanceMatrix(const reco::METCovMatrix& matrix); + void setSumEt(const double & sumEt){ + sumet = sumEt; + }; reco::METCovMatrix getSignificanceMatrix(void) const; private: diff --git a/RecoMET/METPUSubtraction/BuildFile.xml b/RecoMET/METPUSubtraction/BuildFile.xml index 092ffdb18113f..7289d95abf318 100644 --- a/RecoMET/METPUSubtraction/BuildFile.xml +++ b/RecoMET/METPUSubtraction/BuildFile.xml @@ -1,21 +1,50 @@ - - - - + + + + + + + + + + + + - + - - + + + + + + + + + + + + + + - - - + + + + + + + + + + + + diff --git a/RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h b/RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h deleted file mode 100644 index 96e51563305d5..0000000000000 --- a/RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h +++ /dev/null @@ -1,89 +0,0 @@ -#ifndef RecoMET_METPUSubtraction_mvaMEtUtilities_h -#define RecoMET_METPUSubtraction_mvaMEtUtilities_h - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "DataFormats/Candidate/interface/Candidate.h" -#include "DataFormats/METReco/interface/CommonMETData.h" -#include "DataFormats/METReco/interface/PUSubMETData.h" - -#include -#include - -class MvaMEtUtilities -{ - public: - - enum {kPFCands=0,kLeptons,kJets}; - enum {kPF=0, kChHS, kHS, kPU, kHSMinusNeutralPU}; - - private: - - CommonMETData leptonsSum_; - CommonMETData leptonsChSum_; - CommonMETData pfCandSum_; - CommonMETData pfCandChHSSum_; - CommonMETData pfCandChPUSum_; - CommonMETData neutralJetHSSum_; - CommonMETData neutralJetPUSum_; - - std::vector cleanedJets_; - - double dzCut_; - double ptThreshold_; - - public: - - MvaMEtUtilities(const edm::ParameterSet& cfg); - virtual ~MvaMEtUtilities(); - - reco::Candidate::LorentzVector leadJetP4(const std::vector&); - reco::Candidate::LorentzVector subleadJetP4(const std::vector&); - unsigned numJetsAboveThreshold(const std::vector&, double); - - const std::vector& getCleanedJets() const; - - //access functions for lepton suns ============ - double getLeptonsSumMEX() const; - double getLeptonsSumMEY() const; - - double getLeptonsChSumMEX() const; - double getLeptonsChSumMEY() const; - - //recoil and sum computing functions ======== - void computeAllSums(const std::vector& jets, - const std::vector& leptons, - const std::vector& pfCandidates); - - CommonMETData computeRecoil(int metType); - - - protected: - - reco::Candidate::LorentzVector jetP4(const std::vector&, unsigned); - - // cuts on jet Id. MVA output in bins of jet Pt and eta - double mvaCut_[3][4][4]; - - private: - - //utilities functions for jets =============== - bool passesMVA(const reco::Candidate::LorentzVector&, double); - - std::vector cleanJets(const std::vector&, - const std::vector&, double, double); - - //utilities functions for pf candidate ====== - std::vector cleanPFCands(const std::vector&, - const std::vector&, double, bool); - - CommonMETData computeCandSum( int compKey, double dZmax, int dZflag, - bool iCharged, bool mvaPassFlag, - const std::vector& objects ); - - - void finalize(CommonMETData& metData); - -}; - -#endif diff --git a/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h b/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h deleted file mode 100644 index 140546bf052ae..0000000000000 --- a/RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h +++ /dev/null @@ -1,121 +0,0 @@ -#ifndef RecoMET_METPUSubtraction_PFMETAlgorithmMVA_h -#define RecoMET_METPUSubtraction_PFMETAlgorithmMVA_h - -/** \class PFMETAlgorithmMVA - * - * MVA based algorithm for computing the particle-flow missing Et - * - * \authors Phil Harris, CERN - * Christian Veelken, LLR - * - */ - -#include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "CondFormats/EgammaObjects/interface/GBRForest.h" - -#include "DataFormats/Candidate/interface/Candidate.h" -#include "DataFormats/METReco/interface/MET.h" -#include "DataFormats/VertexReco/interface/Vertex.h" - -#include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h" - -//#include -#include - -#include -#include -#include - -class PFMETAlgorithmMVA -{ - - public: - - PFMETAlgorithmMVA(const edm::ParameterSet& cfg); - ~PFMETAlgorithmMVA(); - - void initialize(const edm::EventSetup&); - - void setHasPhotons(bool hasPhotons) { hasPhotons_ = hasPhotons; } - - void setInput(const std::vector&, - const std::vector&, - const std::vector&, - const std::vector&); - - void evaluateMVA(); - - reco::Candidate::LorentzVector getMEt() const { return mvaMEt_; } - const reco::METCovMatrix& getMEtCov() const { return mvaMEtCov_; } - - double getU() const { return mvaOutputU_; } - double getDPhi() const { return mvaOutputDPhi_; } - double getCovU1() const { return mvaOutputCovU1_; } - double getCovU2() const { return mvaOutputCovU2_; } - - void print(std::ostream&) const; - - private: - const std::string updateVariableNames(std::string input); - const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName); - const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName); - - const float evaluateU(); - const float evaluateDPhi(); - const float evaluateCovU1(); - const float evaluateCovU2(); - - MvaMEtUtilities utils_; - - std::string mvaNameU_; - std::string mvaNameDPhi_; - std::string mvaNameCovU1_; - std::string mvaNameCovU2_; - - int mvaType_; - bool hasPhotons_; - - double dZcut_; - std::unique_ptr createFloatVector(std::vector variableNames); - const float GetResponse(const GBRForest *Reader, std::vector &variableNames); - void computeMET(); - std::map var_; - - - float* mvaInputU_; - float* mvaInputDPhi_; - float* mvaInputCovU1_; - float* mvaInputCovU2_; - - float mvaOutputU_; - float mvaOutputDPhi_; - float mvaOutputCovU1_; - float mvaOutputCovU2_; - - std::vector varForU_; - std::vector varForDPhi_; - std::vector varForCovU1_; - std::vector varForCovU2_; - - - double sumLeptonPx_; - double sumLeptonPy_; - double chargedSumLeptonPx_; - double chargedSumLeptonPy_; - - reco::Candidate::LorentzVector mvaMEt_; - //TMatrixD mvaMEtCov_; - reco::METCovMatrix mvaMEtCov_; - - const GBRForest* mvaReaderU_; - const GBRForest* mvaReaderDPhi_; - const GBRForest* mvaReaderCovU1_; - const GBRForest* mvaReaderCovU2_; - - bool loadMVAfromDB_; - - edm::ParameterSet cfg_; -}; -#endif diff --git a/RecoMET/METPUSubtraction/interface/PFMEtSignInterfaceBase.h b/RecoMET/METPUSubtraction/interface/PFMEtSignInterfaceBase.h old mode 100755 new mode 100644 diff --git a/RecoMET/METPUSubtraction/plugins/MVAMET.cc b/RecoMET/METPUSubtraction/plugins/MVAMET.cc new file mode 100644 index 0000000000000..17398a78f0652 --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/MVAMET.cc @@ -0,0 +1,549 @@ +#include + +#include "RecoMET/METPUSubtraction/plugins/MVAMET.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +MVAMET::MVAMET(const edm::ParameterSet& cfg){ + + saveMap_ = cfg.getParameter("saveMap"); + if(saveMap_) + { + produces>(); + produces>(); + } + // get tokens for input METs and prepare for saving the corresponding recoils to the event + srcMETTags_ = cfg.getParameter("srcMETs"); + for(vInputTag::const_iterator it=srcMETTags_.begin();it!=srcMETTags_.end();it++) { + srcMETs_.push_back( consumes( *it ) ); + } + + + // take flags for the met + srcMETFlags_ = cfg.getParameter>("inputMETFlags"); + + if(srcMETFlags_.size() != srcMETTags_.size()) + throw cms::Exception("MVAMET::MVAMET") << " Failed to load MET flags !!\n" << "Expected " << srcMETTags_.size() << " but got " << srcMETFlags_.size() << std::endl; + + //get leptons to calculate Z vector + vInputTag srcLeptonsTags = cfg.getParameter("srcLeptons"); + for(vInputTag::const_iterator it=srcLeptonsTags.begin();it!=srcLeptonsTags.end();it++) { + srcLeptons_.push_back( consumes( *it ) ); + } + + debug_ = (cfg.existsAs("debug")) ? cfg.getParameter("debug") : false; + combineNLeptons_ = cfg.getParameter("combineNLeptons"); + requireOS_ = cfg.getParameter("requireOS"); + // check config: only when exactly two leptons are combined, OS requirement makes sense + if(requireOS_) + assert( combineNLeptons_ == 2 ); + useTauSig_ = cfg.getParameter("useTauSig"); + + srcVertices_ = consumes(cfg.getParameter("srcVertices")); + srcJets_ = consumes(cfg.getParameter("srcJets")); + srcTaus_ = consumes(cfg.getParameter("srcTaus")); + srcMuons_ = consumes(cfg.getParameter("srcMuons")); + if(useTauSig_) + srcTausSignificance_ = consumes::type>(cfg.getParameter("tausSignificance")); + + permuteLeptonsWithinPlugin_ = cfg.getParameter("permuteLeptonsWithinPlugin"); + if (!permuteLeptonsWithinPlugin_) + leptonPermutationsHandle_ = consumes(cfg.getParameter("leptonPermutations")); + + // load weight files + edm::FileInPath weightFile; + weightFile = cfg.getParameter("weightFile"); + mvaReaderPhiCorrection_ = loadMVAfromFile(weightFile, variablesForPhiTraining_, "PhiCorrectedRecoil"); + mvaReaderRecoilCorrection_ = loadMVAfromFile(weightFile, variablesForRecoilTraining_, "LongZCorrectedRecoil"); + mvaReaderCovU1_ = loadMVAfromFile(weightFile, variablesForCovU1_, "CovU1"); + mvaReaderCovU2_ = loadMVAfromFile(weightFile, variablesForCovU2_, "CovU2"); + + // prepare for saving the final mvaMET to the event + mvaMETLabel_ = cfg.getParameter("MVAMETLabel"); + produces(mvaMETLabel_); +} + +MVAMET::~MVAMET(){ + delete mvaReaderPhiCorrection_; + delete mvaReaderRecoilCorrection_; + delete mvaReaderCovU1_; + delete mvaReaderCovU2_; +} + +metPlus MVAMET::calculateRecoil(metPlus* MET, const recoilingBoson &Z, edm::Event& evt) +{ + reco::METCovMatrix rotateToZFrame; + auto Zp4 = Z.p4vec(); + rotateToZFrame(0,0) = rotateToZFrame(1,1) = std::cos(- Zp4.Phi()); + rotateToZFrame(0,1) = std::sin(- Zp4.Phi()); + rotateToZFrame(1,0) = - std::sin(- Zp4.Phi()); + + metPlus Recoil((*MET)); + Recoil.setP4( - Recoil.p4() ); + if( MET->containsCharged() ) + { + Recoil.setP4(Recoil.p4() - Z.chargedP4()); + Recoil.setSumEt(Recoil.sumEt()-Z.chargedSumEt()); + } + + if( MET->containsNeutral() ) + Recoil.setSumEt(Recoil.sumEt()-Z.neutralSumEt()); + + reco::METCovMatrix rotatedCovMatrix = rotateToZFrame * Recoil.getSignificanceMatrix(); + Recoil.setSignificanceMatrix( rotatedCovMatrix ); + + return Recoil; +} + +void MVAMET::doCombinations(int offset, int k) +{ + if (k == 0) + { + combinations_.push_back(combination_); + combination_.pop_back(); + return; + } + for (size_t i = offset; i <= allLeptons_.size() - k; ++i) + { + combination_.push_back(allLeptons_[i]); + doCombinations(i+1, k-1); + } + combination_.clear(); +} + +void MVAMET::unpackCompositeCands(edm::Event& evt) +{ + edm::Handle ccands; + evt.getByToken(leptonPermutationsHandle_, ccands); + for (auto ccand = ccands->begin(); ccand != ccands->end(); ++ccand) + { + combination_.clear(); + for (unsigned int ielem = 0; ielem < combineNLeptons_; ++ielem) + { + const reco::Candidate& cand = *(ccand->daughter(ielem)); + combination_.push_back(cand.masterClone().castTo >()); + } + combinations_.push_back(combination_); + } +} + +void MVAMET::handleTaus(edm::Ptr lepton, recoilingBoson& Z, const pat::TauCollection& tauCollection) +{ + recoilComponent rComp(lepton); + for(const auto & tau : tauCollection) + { + if(deltaR2 (*lepton, tau) > 1.e-6) // dR > 0.001 + continue; + for(const auto & candidate : tau.signalCands()) + { + if(abs(candidate->pdgId()) > 11 and abs(candidate->pdgId()) < 16) + continue; + + if(candidate->charge() !=0) + rComp.chargedTauJetCandidates.push_back(candidate); + else + rComp.neutralTauJetCandidates.push_back(candidate); + } + } + Z.addLepton(rComp); +} + +void MVAMET::handleMuons(edm::Ptr lepton, recoilingBoson& Z, const pat::MuonCollection& muCollection) +{ + math::PtEtaPhiELorentzVectorD p4Photon; + recoilComponent rComp(lepton); + for (const auto & muon : muCollection) + { + if(deltaR2 (*lepton, muon) < 1.e-6) // dR < 0.001 + { + p4Photon.SetPt(muon.pfEcalEnergy()/TMath::CosH(muon.p4().eta())); + p4Photon.SetEta(muon.p4().eta()); + p4Photon.SetPhi(muon.p4().phi()); + p4Photon.SetE(muon.pfEcalEnergy()); + } + } + + if(p4Photon.E() > 0 ) + rComp.setP4(lepton->p4() + p4Photon); + else + rComp.setP4(lepton->p4()); + + Z.addLepton(rComp); +} + +void MVAMET::calculateRecoilingObjects(edm::Event &evt, const pat::MuonCollection& muCollection, const pat::TauCollection& tauCollection) +{ + // loop on all the indentified leptons and put them in a common vector allLeptons_ + allLeptons_.clear(); + combinations_.clear(); + for ( std::vector >::const_iterator srcLeptons_i = srcLeptons_.begin(); srcLeptons_i != srcLeptons_.end(); ++srcLeptons_i ) + { + edm::Handle leptons; + evt.getByToken(*srcLeptons_i, leptons); + for (size_t i=0; i < leptons->size(); ++i) + allLeptons_.push_back(leptons->ptrAt(i)); + } + + if (!permuteLeptonsWithinPlugin_) + unpackCompositeCands(evt); + + else if(allLeptons_.size() >= combineNLeptons_) + doCombinations(0, combineNLeptons_); + + if(debug_) + std::cout << allLeptons_.size() << " lead to " << combinations_.size(); + + if(requireOS_) + cleanLeptonsFromSS(); + + if(debug_) + std::cout << " with OS: " << combinations_.size() << std::endl; + + for(auto leptonpair: combinations_) + { + recoilingBoson Z; + for(auto lepton : leptonpair) + { + if(abs(lepton->pdgId()) == 13) + handleMuons( lepton, Z, muCollection); + else if(abs(lepton->pdgId()) == 15) + handleTaus( lepton, Z, tauCollection); + else if(abs(lepton->pdgId())==11) + { + recoilComponent rComp(lepton); + rComp.setP4(lepton->p4()); + Z.addLepton(rComp); + } + else + std::cout << "Warning: unsupported recoiling object found" << std::endl; + } + Bosons_.push_back(Z); + } +} + +void MVAMET::fillEventInformation(edm::Event& evt) +{ + edm::Handle jets; + evt.getByToken(srcJets_, jets); + size_t jetsSize = jets->size(); + for( size_t iJet = 0; iJet <= 1; iJet++) + { + var_["Jet" + std::to_string(iJet)+ "_Pt"] = (jetsSize > iJet) ? jets->at(iJet).p4().pt() : 0; + var_["Jet" + std::to_string(iJet)+ "_Eta"] = (jetsSize > iJet) ? jets->at(iJet).p4().Eta() : 0; + var_["Jet" + std::to_string(iJet)+ "_Phi"] = (jetsSize > iJet) ? jets->at(iJet).p4().Phi() : 0; + var_["Jet" + std::to_string(iJet)+ "_M"] = (jetsSize > iJet) ? jets->at(iJet).p4().M() : 0; + } + + var_["NCleanedJets"] = countJets(*jets, 10); + var_["nCombinations"] = combinations_.size(); + + // treat other collections and save to map + edm::Handle vertices; + evt.getByToken(srcVertices_, vertices); + var_["NVertex"] = vertices->size(); +} + +void MVAMET::TagZ() +{ + if(Bosons_.size() == 0) + return; + + std::vector::iterator taggedBoson = Bosons_.begin(); + for(std::vector::iterator Z= Bosons_.begin(); Z != Bosons_.end(); Z++) + { + if(taggedBoson->p4vec().pt() < Z->p4vec().pt()) + taggedBoson = Z; + } + taggedBoson->setTagged(); +} + +void MVAMET::produce(edm::Event& evt, const edm::EventSetup& es){ + var_.clear(); + Bosons_.clear(); + + // get taus + edm::Handle tauCollectionHandle; + evt.getByToken(srcTaus_, tauCollectionHandle); + const pat::TauCollection tauCollection = *(tauCollectionHandle.product()); + + // get muons + edm::Handle muCollectionHandle; + evt.getByToken(srcMuons_, muCollectionHandle); + const pat::MuonCollection muCollection = *(muCollectionHandle.product()); + + // read in additional taus significance contribution + edm::Handle::type> tausSignificance; + reco::METCovMatrix tausMatrix; + if(useTauSig_) + { + evt.getByToken(srcTausSignificance_, tausSignificance); + tausMatrix(0,0) = (*tausSignificance)(0,0); + tausMatrix(1,0) = (*tausSignificance)(1,0); + tausMatrix(0,1) = (*tausSignificance)(0,1); + tausMatrix(1,1) = (*tausSignificance)(1,1); + } + + //fill allLeptons_ + calculateRecoilingObjects(evt, muCollection, tauCollection); + + //fill event meta information + fillEventInformation(evt); + // create output collections + std::auto_ptr recoilpatMETCollection(new pat::METCollection()); + std::auto_ptr patMETCollection(new pat::METCollection()); + + // stop execution if no recoiling object has been found + if(Bosons_.size() == 0) + { + if(saveMap_) + saveMap(evt); + evt.put(patMETCollection,mvaMETLabel_); + return; + } + + if(saveMap_) + TagZ(); + + std::vector >::const_iterator referenceMET = srcMETs_.begin(); + evt.getByToken(*referenceMET, referenceMETHandle_); + // loop on identified combinations of recoiling objects, here donted as "Z" + for(const auto & Z: Bosons_) + { + metPlus referenceRecoil; + std::vector::const_iterator itMETFlags = srcMETFlags_.begin(); + + // MET flags show what components have to be subtracted to get the recoil from the MET + // MET Flags: 0 -> Neutral + Charged PV + // 1 -> only charged PV + // 2 -> only Neutral PV + // 3 -> not included in MET + // + int i = 0; + for ( std::vector >::const_iterator srcMET = srcMETs_.begin(); srcMET != srcMETs_.end() && itMETFlags!=srcMETFlags_.end(); ++srcMET, ++itMETFlags ) + { + edm::Handle METhandle; + evt.getByToken(*srcMET, METhandle); + assert((*METhandle).size() == 1 ); + + metPlus MET((*METhandle)[0]); + MET.collection_name = srcMETTags_[i].label(); + MET.METFlag = (*itMETFlags); + + // to be removed begin + if(i==0) + { + if(saveMap_) + { + genMET_ = (*METhandle)[0].genMET(); + var_["genMet_Pt"] = genMET_->p4().pt(); + var_["genMet_Phi"] = genMET_->p4().phi(); + } + } + // to be removed end + + auto Recoil = calculateRecoil(&MET, Z, evt); + if(i == 0) + referenceRecoil = Recoil; + + addToMap(Recoil, referenceRecoil); + if(saveMap_) + addToMap(Recoil, Z); + ++i; + } + + // evaluate phi training and apply angular correction + Float_t PhiAngle = GetResponse(mvaReaderPhiCorrection_, variablesForPhiTraining_); + + auto refRecoil = TVector2(referenceRecoil.p4().px(), referenceRecoil.p4().py()); + refRecoil = refRecoil.Rotate(PhiAngle); + reco::Candidate::LorentzVector phiCorrectedRecoil(refRecoil.Px(), refRecoil.Py(), 0, referenceRecoil.sumEt()); + addToMap(phiCorrectedRecoil, referenceRecoil.sumEt(), "PhiCorrectedRecoil"); + + // evaluate second training and apply recoil correction + Float_t RecoilCorrection = GetResponse(mvaReaderRecoilCorrection_, variablesForRecoilTraining_); + refRecoil *= RecoilCorrection; + + // create pat::MET from recoil + pat::MET recoilmvaMET(referenceRecoil); + reco::Candidate::LorentzVector recoilP4(refRecoil.Px(), refRecoil.Py(), 0, referenceRecoil.sumEt()); + recoilmvaMET.setP4(recoilP4); + addToMap(recoilP4, referenceRecoil.sumEt(), "LongZCorrectedRecoil"); + + // evaluate covariance matrix regression + Float_t CovU1Correction = GetResponse(mvaReaderCovU1_, variablesForCovU1_); + Float_t CovU1 = CovU1Correction * refRecoil.Mod(); + Float_t CovU2Correction = GetResponse(mvaReaderCovU2_, variablesForCovU2_); + Float_t CovU2 = CovU2Correction * refRecoil.Mod(); + + //// save results to event + recoilpatMETCollection->push_back(recoilmvaMET); + + // calculate new mvaMET + pat::MET mvaMET(referenceRecoil); + mvaMET.setP4(recoilP4); + mvaMET.setP4(mvaMET.p4() + Z.chargedP4()); + mvaMET.setP4(- mvaMET.p4()); + // copy sumEt from input MET since there's no correction in MVA MET for that + mvaMET.setSumEt((*referenceMETHandle_)[0].sumEt()); + reco::METCovMatrix mvaMETCov; + double cosPhi = std::cos(refRecoil.Phi()); + double sinPhi = std::sin(refRecoil.Phi()); + mvaMETCov(0, 0) = std::pow(CovU1, 2)*cosPhi*cosPhi + std::pow(CovU2, 2)*sinPhi*sinPhi; + mvaMETCov(0, 1) = -std::pow(CovU1, 2) * sinPhi*cosPhi + std::pow(CovU2, 2) * sinPhi*cosPhi; + mvaMETCov(1, 0) = mvaMETCov(0, 1); + mvaMETCov(1, 1) = std::pow(CovU1, 2) * sinPhi*sinPhi + std::pow(CovU2, 2) * cosPhi*cosPhi; + + if(useTauSig_) + mvaMET.setSignificanceMatrix(mvaMETCov + tausMatrix); + else + mvaMET.setSignificanceMatrix(mvaMETCov); + + // add constituent info to pat::MET ; only if combinatorics done here (else potential problems with edm::Ptr cast) + if (permuteLeptonsWithinPlugin_) + { + size_t iCount=0; + for(const auto & lepton: Z.leptons) + mvaMET.addUserCand("lepton" + std::to_string(iCount++), lepton.getSrcLepton()); + } + + // save MVA MET results + mvaMET.addUserFloat("PhiCorrection", PhiAngle); + mvaMET.addUserFloat("LongZCorrection", RecoilCorrection); + mvaMET.addUserFloat("CovU1", CovU1Correction); + mvaMET.addUserFloat("CovU2", CovU2Correction); + patMETCollection->push_back(mvaMET); + + // muon selection for training + if(saveMap_) + { + var_["PhiTrainingResponse"] = PhiAngle; + var_["RecoilTrainingResponse"] = RecoilCorrection; + var_["MVAMET_Pt"] = mvaMET.p4().pt(); + var_["MVAMET_Phi"] = mvaMET.p4().phi() ; + var_["MVAMET_sumEt"] = mvaMET.sumEt(); + + reco::Candidate::LorentzVector dmvamet = genMET_->p4() - mvaMET.p4(); + var_["dmvamet_Pt"] = dmvamet.Pt(); + var_["dmvamet_Phi"] = dmvamet.Phi(); + reco::Candidate::LorentzVector dpfmet = genMET_->p4() - (*referenceMETHandle_)[0].p4(); + var_["dpfmet_Pt"] = dpfmet.Pt(); + var_["dpfmet_Phi"] = dpfmet.Phi(); + if(Z.isTagged()) + { + addToMap(Z); + var_["select"] = Z.select(); + saveMap(evt); + } + } + } + if(debug_) + { + for(const auto & entry : var_) + std::cout << entry.first << " : " << entry.second << std::endl; + } + evt.put(patMETCollection, mvaMETLabel_); +} + +void MVAMET::saveMap(edm::Event& evt) +{ + std::auto_ptr> variableNames(new std::vector); + std::auto_ptr > variables(new std::vector); + + for(const auto & entry : var_){ + variableNames->push_back(entry.first); + variables->push_back(entry.second); + } + + evt.put(variableNames); + evt.put(variables); +} + +void MVAMET::addToMap(const metPlus &recoil, const recoilingBoson &Z) +{ + reco::Candidate::LorentzVector p4(Z.chargedP4() ); + TVector2 diLeptonMomentum(p4.X(), p4.Y()); + TVector2 recoilT2(recoil.p4().X(), recoil.p4().Y()); + TVector2 rotatedRecoil = recoilT2.Rotate(- diLeptonMomentum.Phi()); + auto type = "recoil" + recoil.collection_name; + var_[type + "_LongZ" ] = rotatedRecoil.X(); + var_[type + "_PerpZ" ] = rotatedRecoil.Y(); + +} + +void MVAMET::addToMap(const metPlus &recoil, const metPlus &referenceMET) +{ + auto type = "recoil" + recoil.collection_name; + addToMap(recoil.p4(), recoil.sumEt(), type); + var_[type + "_Cov00" ] = recoil.getSignificanceMatrix()(0,0); + var_[type + "_Cov11" ] = recoil.getSignificanceMatrix()(1,1); + var_[type + "_sumEtFraction" ] = recoil.sumEt()/referenceMET.sumEt(); + var_[type + "_dPhi" ] = TVector2::Phi_mpi_pi(recoil.phi() - referenceMET.phi()); +} + +void MVAMET::addToMap(const recoilingBoson &Z) +{ + reco::Candidate::LorentzVector Zp4 = Z.chargedP4(); + var_["Boson_Pt"] = Zp4.pt(); + var_["Boson_Phi"] = Zp4.Phi(); + var_["Boson_Eta"] = Zp4.Eta(); + var_["Boson_M"] = Zp4.M(); + var_["Boson_sumET"] = Z.chargedSumEt(); + var_["Boson_daughter1"] = float(Z.getPdgId(0)); + var_["Boson_daughter2"] = float(Z.getPdgId(1)); + +} + +void MVAMET::addToMap(const reco::Candidate::LorentzVector p4, const double sumEt, const std::string &type) +{ + var_[type + "_Pt" ] = p4.pt(); + var_[type + "_Phi" ] = p4.phi(); + var_[type + "_sumEt" ] = sumEt; +} + +unsigned int MVAMET::countJets(const pat::JetCollection& jets, const float maxPt){ + int nJets = 0; + for(const auto & jet : jets){ + if(jet.pt() > maxPt) + nJets++; + } + return nJets; +} + +const Float_t MVAMET::GetResponse(const GBRForest * Reader, std::vector &variableNames ){ + + Float_t* mvaInputVector = createFloatVector(variableNames); + double result = Reader->GetResponse(mvaInputVector); + delete mvaInputVector; + return result; +} + + +const GBRForest* MVAMET::loadMVAfromFile(const edm::FileInPath& inputFileName, std::vector& trainingVariableNames, std::string mvaName){ + + if ( inputFileName.location()==edm::FileInPath::Unknown ) + throw cms::Exception("PFMETAlgorithmMVA::loadMVA") << " Failed to find File = " << inputFileName << " !!\n"; + + std::unique_ptr inputFile(new TFile(inputFileName.fullPath().data())); + std::string variableListName = mvaName + "varlist"; + std::unique_ptr> lVec ((std::vector*)inputFile->Get(variableListName.c_str())); + + for(unsigned int i=0; i< lVec->size();++i) + { + trainingVariableNames.push_back(lVec->at(i)); + } + const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data()); + if ( !mva ) + throw cms::Exception("PFMETAlgorithmMVA::loadMVA") << " Failed to load MVA from file = " << inputFileName.fullPath().data() << " !!\n"; + + return mva; +} + +Float_t* MVAMET::createFloatVector(std::vector variableNames){ + Float_t* floatVector = new Float_t[variableNames.size()]; + for(size_t i = 0; i < variableNames.size(); ++i){ + floatVector[i] = var_[variableNames[i]]; + } + return floatVector; +} + + +DEFINE_FWK_MODULE(MVAMET); diff --git a/RecoMET/METPUSubtraction/plugins/MVAMET.h b/RecoMET/METPUSubtraction/plugins/MVAMET.h new file mode 100644 index 0000000000000..885c2f06e09d7 --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/MVAMET.h @@ -0,0 +1,230 @@ +#ifndef MVAMET_h +#define MVAMET_h + +/** \class MVAMET + * */ + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include +#include "DataFormats/METReco/interface/PFMET.h" +#include "DataFormats/METReco/interface/PFMETCollection.h" +#include +#include +#include +#include +#include +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/PatCandidates/interface/CompositeCandidate.h" +#include "CondFormats/EgammaObjects/interface/GBRForest.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" + +#include +#include +#include + +// class storing the MET and recoil information +class metPlus : public pat::MET { +public: + int METFlag; + std::string collection_name; + metPlus() {} + metPlus(pat::MET mother) : pat::MET(mother), + METFlag(-1), + collection_name("unset") {} + //function to return METFlags that tell whether the MET was made from PV charged and neutral components + bool containsCharged() { return ((this->METFlag == 0) || (this->METFlag == 1)); } + bool containsNeutral() { return ((this->METFlag == 0) || (this->METFlag == 2)); } +}; + + +// class to store constituents of the recoil +// because the tau is a composite object this information can be stored and accessed with this class +class recoilComponent { + private: + edm::Ptr srcLepton_; + reco::Candidate::LorentzVector p4_; + public: + std::vector chargedTauJetCandidates; + std::vector neutralTauJetCandidates; + recoilComponent(edm::Ptr srcLepton) : srcLepton_(srcLepton) {} + edm::Ptr getSrcLepton() const { return srcLepton_; } + reco::Candidate::LorentzVector p4() const { return this->p4_; } + void setP4(const reco::Candidate::LorentzVector & p4) { p4_ = p4; } + reco::Candidate::LorentzVector chargedP4() const + { + reco::Candidate::LorentzVector p4; + for(const auto & tauJet: chargedTauJetCandidates) + p4 += tauJet->p4(); + return p4+this->p4(); + } + + reco::Candidate::LorentzVector neutralP4() const + { + reco::Candidate::LorentzVector p4; + for(const auto & tauJet: neutralTauJetCandidates) + p4 += tauJet->p4(); + return p4; + } + + float chargedSumEt() const + { + float sumEt = 0; + for(const auto & tauJet: chargedTauJetCandidates) + sumEt += tauJet->p4().pt(); + return sumEt + this->p4().pt(); + } + + float neutralSumEt() const + { + float sumEt = 0; + for(const auto & tauJet: neutralTauJetCandidates) + sumEt += tauJet->p4().pt(); + return sumEt; + } + int pdgId() const { return this->srcLepton_->pdgId(); } + bool isMuon() const { return (abs(pdgId()) == 13); } +}; + +class recoilingBoson : public reco::Particle { + // tag: taggs the combination with the highest pt + // might be dropped since 3rd lepton veto is applied anyway for the training + bool tag; + + public: + std::vector leptons; + recoilingBoson() : tag(false) {} + void addLepton(recoilComponent rComp) { this->leptons.push_back(rComp); } + bool isDiMuon() const { return (leptons.size() == 2) ? (leptons[0].isMuon() and leptons[1].isMuon()) : false; } + bool select() const { return this->p4vec().M() > 80 && this->p4vec().M() < 100 && this->isDiMuon(); } + void setTagged() { this->tag = true; } + bool isTagged() const { return this->tag; } + + reco::Candidate::LorentzVector p4vec() const { return (this->chargedP4() + this->neutralP4()); } + + reco::Candidate::LorentzVector chargedP4() const + { + reco::Candidate::LorentzVector p4; + for(const auto & lepton: leptons) + p4 += lepton.chargedP4(); + + return p4; + } + + reco::Candidate::LorentzVector neutralP4() const + { + reco::Candidate::LorentzVector p4; + for(const auto & lepton: leptons) + p4 += lepton.neutralP4(); + + return p4; + } + + float chargedSumEt() const + { + float sumEt = 0; + for(const auto & lepton: leptons) + sumEt += lepton.chargedSumEt(); + + return sumEt; + } + + float neutralSumEt() const + { + float sumEt = 0; + for(const auto & lepton: leptons) + sumEt += lepton.neutralSumEt(); + + return sumEt; + } + + double sumEt() const { return (chargedSumEt() + neutralSumEt()); } + int getPdgId(int index) const { return leptons[index].pdgId(); } +}; + +class MVAMET : public edm::stream::EDProducer<> { + + public: + + // basic constructor from parameter set + MVAMET(const edm::ParameterSet&); + ~MVAMET(); + + void produce(edm::Event&, const edm::EventSetup&); + typedef std::vector vInputTag; + float bestMass_; + // create a vector given input variables + Float_t* createFloatVector(std::vector variableNames); + + unsigned int countJets(const pat::JetCollection& jets, const float maxPt); + + // load MVA file produced in the training + const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName, std::vector& trainingVariableNames, std::string mvaName); + // read the response + const Float_t GetResponse(const GBRForest * Reader,std::vector &variableNames ); + + // to correctly create the map of regression input vriables + void addToMap(const reco::Candidate::LorentzVector p4, const double sumEt, const std::string &type); + void addToMap(const metPlus &recoil, const metPlus &referenceMET); + void addToMap(const recoilingBoson &Z); + void addToMap(const metPlus &recoil, const recoilingBoson &Z); + + metPlus calculateRecoil(metPlus* MET, const recoilingBoson &Z, edm::Event& evt); + void TagZ(); +private: + void doCombinations(int offset, int k); + void unpackCompositeCands(edm::Event& evt); + void saveMap(edm::Event& evt); + void calculateRecoilingObjects(edm::Event& evt, const pat::MuonCollection&, const pat::TauCollection& ); + void cleanLeptonsFromSS() + { + combinations_.erase(std::remove_if(combinations_.begin(), combinations_.end(), [](std::vector> pair) { return pair[0]->charge() == pair[1]->charge(); }), combinations_.end()); + } + void handleMuons(edm::Ptr lepton, recoilingBoson& Z, const pat::MuonCollection& ); + void handleTaus(edm::Ptr lepton, recoilingBoson& Z, const pat::TauCollection& ); + void fillEventInformation(edm::Event&); + std::string mvaMETLabel_; + + vInputTag srcMETTags_; + + std::vector > srcMETs_; + edm::EDGetTokenT srcVertices_; + edm::EDGetTokenT srcJets_; + std::vector > srcLeptons_; + edm::EDGetTokenT srcTaus_; + edm::EDGetTokenT srcMuons_; + edm::EDGetTokenT::type> srcTausSignificance_; + bool useTauSig_; + + std::vector srcMETFlags_; + + std::map var_; + + std::vector> allLeptons_; + std::vector>> combinations_; + std::vector> combination_; + + std::vector variablesForPhiTraining_ = {}; + std::vector variablesForRecoilTraining_ = {}; + std::vector variablesForCovU1_ = {}; + std::vector variablesForCovU2_ = {}; + + const GBRForest* mvaReaderPhiCorrection_; + const GBRForest* mvaReaderRecoilCorrection_; + const GBRForest* mvaReaderCovU1_; + const GBRForest* mvaReaderCovU2_; + + bool debug_; + bool saveMap_; + bool produceRecoils_; + std::vector Bosons_; + size_t combineNLeptons_; + bool requireOS_; + edm::Handle referenceMETHandle_; + bool permuteLeptonsWithinPlugin_; + edm::EDGetTokenT leptonPermutationsHandle_; + // to be removed + const reco::GenMET * genMET_; +}; +#endif diff --git a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.cc b/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.cc deleted file mode 100644 index b0b36de55386c..0000000000000 --- a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.cc +++ /dev/null @@ -1,386 +0,0 @@ -#include "RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h" - -using namespace reco; - -const double dR2Min = 0.01*0.01; -const double dR2Max = 0.5*0.5; -const double dPtMatch = 0.1; - -PFMETProducerMVA::PFMETProducerMVA(const edm::ParameterSet& cfg) - : mvaMEtAlgo_(cfg), - mvaMEtAlgo_isInitialized_(false) -{ - srcCorrJets_ = consumes(cfg.getParameter("srcCorrJets")); - srcUncorrJets_ = consumes(cfg.getParameter("srcUncorrJets")); - srcJetIds_ = consumes >(cfg.getParameter("srcMVAPileupJetId")); - srcPFCandidatesView_ = consumes(cfg.getParameter("srcPFCandidates")); - srcVertices_ = consumes(cfg.getParameter("srcVertices")); - vInputTag srcLeptonsTags = cfg.getParameter("srcLeptons"); - for(vInputTag::const_iterator it=srcLeptonsTags.begin();it!=srcLeptonsTags.end();it++) { - srcLeptons_.push_back( consumes( *it ) ); - } - mJetCorrector_ = consumes(cfg.getParameter("corrector")); - minNumLeptons_ = cfg.getParameter("minNumLeptons"); - - globalThreshold_ = cfg.getParameter("globalThreshold"); - - minCorrJetPt_ = cfg.getParameter ("minCorrJetPt"); - useType1_ = cfg.getParameter ("useType1"); - - verbosity_ = ( cfg.exists("verbosity") ) ? - cfg.getParameter("verbosity") : 0; - - produces(); -} - -PFMETProducerMVA::~PFMETProducerMVA(){} - -void PFMETProducerMVA::produce(edm::Event& evt, const edm::EventSetup& es) -{ - // CV: check if the event is to be skipped - if ( minNumLeptons_ > 0 ) { - int numLeptons = 0; - for ( std::vector >::const_iterator srcLeptons_i = srcLeptons_.begin(); - srcLeptons_i != srcLeptons_.end(); ++srcLeptons_i ) { - edm::Handle leptons; - evt.getByToken(*srcLeptons_i, leptons); - numLeptons += leptons->size(); - } - if ( !(numLeptons >= minNumLeptons_) ) { - LogDebug( "produce" ) - << ":" << std::endl - << "Run: " << evt.id().run() << ", LS: " << evt.luminosityBlock() << ", Event: " << evt.id().event() << std::endl - << " numLeptons = " << numLeptons << ", minNumLeptons = " << minNumLeptons_ << " --> skipping !!" << std::endl; - - reco::PFMET pfMEt; - std::auto_ptr pfMEtCollection(new reco::PFMETCollection()); - pfMEtCollection->push_back(pfMEt); - evt.put(pfMEtCollection); - return; - } - } - - //get jet IDs - edm::Handle > jetIds; - evt.getByToken(srcJetIds_, jetIds); - - // get jets (corrected and uncorrected) - edm::Handle corrJets; - evt.getByToken(srcCorrJets_, corrJets); - - edm::Handle uncorrJets; - evt.getByToken(srcUncorrJets_, uncorrJets); - - edm::Handle corrector; - if( useType1_ ) - { - evt.getByToken(mJetCorrector_, corrector); - } - - edm::Handle pfCandidates_view; - evt.getByToken(srcPFCandidatesView_, pfCandidates_view); - - // get vertices - edm::Handle vertices; - evt.getByToken(srcVertices_, vertices); - // take vertex with highest sum(trackPt) as the vertex of the "hard scatter" interaction - // (= first entry in vertex collection) - const reco::Vertex* hardScatterVertex = ( vertices->size() >= 1 ) ? - &(vertices->front()) : nullptr; - - // get leptons - // (excluded from sum over PFCandidates when computing hadronic recoil) - int lId = 0; - bool lHasPhotons = false; - std::vector leptonInfo = computeLeptonInfo(srcLeptons_,*pfCandidates_view,hardScatterVertex, lId, lHasPhotons, - evt); - - // initialize MVA MET algorithm - // (this will load the BDTs, stored as GBRForrest objects; - // either in input ROOT files or in SQL-lite files/the Conditions Database) - if ( !mvaMEtAlgo_isInitialized_ ) { - mvaMEtAlgo_.initialize(es); - mvaMEtAlgo_isInitialized_ = true; - } - - // reconstruct "standard" particle-flow missing Et - CommonMETData pfMEt_data = metAlgo_.run( (*pfCandidates_view), globalThreshold_); - SpecificPFMETData specificPfMET = pfMEtSpecificAlgo_.run( (*pfCandidates_view) ); - const reco::Candidate::LorentzVector p4( pfMEt_data.mex, pfMEt_data.mey, 0.0, pfMEt_data.met); - const reco::Candidate::Point vtx(0.0, 0.0, 0.0 ); - reco::PFMET pfMEt(specificPfMET,pfMEt_data.sumet, p4, vtx); - reco::Candidate::LorentzVector pfMEtP4_original = pfMEt.p4(); - - // compute objects specific to MVA based MET reconstruction - std::vector pfCandidateInfo = computePFCandidateInfo(*pfCandidates_view, hardScatterVertex); - std::vector jetInfo = computeJetInfo(*uncorrJets, corrJets, *jetIds, *vertices, hardScatterVertex, *corrector,evt,es,leptonInfo,pfCandidateInfo); - std::vector vertexInfo = computeVertexInfo(*vertices); - // compute MVA based MET and estimate of its uncertainty - mvaMEtAlgo_.setInput(leptonInfo, jetInfo, pfCandidateInfo, vertexInfo); - mvaMEtAlgo_.setHasPhotons(lHasPhotons); - mvaMEtAlgo_.evaluateMVA(); - pfMEt.setP4(mvaMEtAlgo_.getMEt()); - pfMEt.setSignificanceMatrix(mvaMEtAlgo_.getMEtCov()); - - LogDebug("produce") - << "Run: " << evt.id().run() << ", LS: " << evt.luminosityBlock() << ", Event: " << evt.id().event() << std::endl - << " PFMET: Pt = " << pfMEtP4_original.pt() << ", phi = " << pfMEtP4_original.phi() << " " - << "(Px = " << pfMEtP4_original.px() << ", Py = " << pfMEtP4_original.py() << ")" << std::endl - << " MVA MET: Pt = " << pfMEt.pt() << " phi = " << pfMEt.phi() << " (Px = " << pfMEt.px() << ", Py = " << pfMEt.py() << ")" << std::endl - << " Cov:" << std::endl - <<(mvaMEtAlgo_.getMEtCov())(0,0)<<" "<<(mvaMEtAlgo_.getMEtCov())(0,1)< pfMEtCollection(new reco::PFMETCollection()); - pfMEtCollection->push_back(pfMEt); - evt.put(pfMEtCollection); -} - -std::vector -PFMETProducerMVA::computeLeptonInfo(const std::vector >& srcLeptons_, - const reco::CandidateView& pfCandidates_view, - const reco::Vertex* hardScatterVertex, - int& lId, bool& lHasPhotons, edm::Event& evt ) { - - std::vector leptonInfo; - - for ( std::vector >::const_iterator srcLeptons_i = srcLeptons_.begin(); - srcLeptons_i != srcLeptons_.end(); ++srcLeptons_i ) { - edm::Handle leptons; - evt.getByToken(*srcLeptons_i, leptons); - for ( reco::CandidateView::const_iterator lepton1 = leptons->begin(); - lepton1 != leptons->end(); ++lepton1 ) { - bool pMatch = false; - for ( std::vector >::const_iterator srcLeptons_j = srcLeptons_.begin(); - srcLeptons_j != srcLeptons_.end(); ++srcLeptons_j ) { - edm::Handle leptons2; - evt.getByToken(*srcLeptons_j, leptons2); - for ( reco::CandidateView::const_iterator lepton2 = leptons2->begin(); - lepton2 != leptons2->end(); ++lepton2 ) { - if(&(*lepton1) == &(*lepton2)) { continue; } - if(deltaR2(lepton1->p4(),lepton2->p4()) < dR2Max) { pMatch = true; } - if(pMatch && !istau(&(*lepton1)) && istau(&(*lepton2))) { pMatch = false; } - if(pMatch && ( (istau(&(*lepton1)) && istau(&(*lepton2))) || (!istau(&(*lepton1)) && !istau(&(*lepton2)))) - && lepton1->pt() > lepton2->pt()) { pMatch = false; } - if(pMatch && lepton1->pt() == lepton2->pt()) { - pMatch = false; - for(unsigned int i0 = 0; i0 < leptonInfo.size(); i0++) { - if(std::abs(lepton1->pt() - leptonInfo[i0].p4().pt()) < dPtMatch) { pMatch = true; break; } - } - } - if(pMatch) break; - } - if(pMatch) break; - } - if(pMatch) continue; - reco::PUSubMETCandInfo pLeptonInfo; - pLeptonInfo.setP4( lepton1->p4() ); - pLeptonInfo.setChargedEnFrac( chargedEnFrac(&(*lepton1),pfCandidates_view,hardScatterVertex) ); - leptonInfo.push_back(pLeptonInfo); - if(lepton1->isPhoton()) { lHasPhotons = true; } - } - lId++; - } - - return leptonInfo; -} - - -std::vector -PFMETProducerMVA::computeJetInfo(const reco::PFJetCollection& uncorrJets, - const edm::Handle& corrJets, - const edm::ValueMap& jetIds, - const reco::VertexCollection& vertices, - const reco::Vertex* hardScatterVertex, - const reco::JetCorrector &iCorrector,edm::Event &iEvent,const edm::EventSetup &iSetup, - std::vector &iLeptons,std::vector &iCands) -{ - std::vector retVal; - for ( reco::PFJetCollection::const_iterator uncorrJet = uncorrJets.begin(); - uncorrJet != uncorrJets.end(); ++uncorrJet ) { - // for ( reco::PFJetCollection::const_iterator corrJet = corrJets.begin(); - // corrJet != corrJets.end(); ++corrJet ) { - auto corrJet = corrJets->begin(); - for( size_t cjIdx=0;cjIdxsize();++cjIdx, ++corrJet) { - reco::PFJetRef corrJetRef( corrJets, cjIdx ); - - // match corrected and uncorrected jets - if ( uncorrJet->jetArea() != corrJet->jetArea() ) continue; - if ( deltaR2(corrJet->p4(),uncorrJet->p4()) > dR2Min ) continue; - - // check that jet passes loose PFJet id. - if(!passPFLooseId(&(*uncorrJet))) continue; - - // compute jet energy correction factor - // (= ratio of corrected/uncorrected jet Pt) - //double jetEnCorrFactor = corrJet->pt()/uncorrJet->pt(); - reco::PUSubMETCandInfo jetInfo; - - // PH: apply jet energy corrections for all Jets ignoring recommendations - jetInfo.setP4( corrJet->p4() ); - double lType1Corr = 0; - if(useType1_) { //Compute the type 1 correction ===> This code is crap - double pCorr = iCorrector.correction(*uncorrJet); - lType1Corr = std::abs(corrJet->pt()-pCorr*uncorrJet->pt()); - TLorentzVector pVec; pVec.SetPtEtaPhiM(lType1Corr,0,corrJet->phi(),0); - reco::Candidate::LorentzVector pType1Corr; pType1Corr.SetCoordinates(pVec.Px(),pVec.Py(),pVec.Pz(),pVec.E()); - //Filter to leptons - bool pOnLepton = false; - for(unsigned int i0 = 0; i0 < iLeptons.size(); i0++) { - if(deltaR2(iLeptons[i0].p4(),corrJet->p4()) < dR2Max) {pOnLepton = true; break;} - } - //Add it to PF Collection - if(corrJet->pt() > 10 && !pOnLepton) { - reco::PUSubMETCandInfo pfCandidateInfo; - pfCandidateInfo.setP4( pType1Corr ); - pfCandidateInfo.setDZ( -999 ); - iCands.push_back(pfCandidateInfo); - } - //Scale - lType1Corr = (pCorr*uncorrJet->pt()-uncorrJet->pt()); - lType1Corr /=corrJet->pt(); - } - - // check that jet Pt used to compute MVA based jet id. is above threshold - if ( !(jetInfo.p4().pt() > minCorrJetPt_) ) continue; - - - jetInfo.setMvaVal( jetIds[ corrJetRef ] ); - float chEnF = (uncorrJet->chargedEmEnergy() + uncorrJet->chargedHadronEnergy() + uncorrJet->chargedMuEnergy() )/uncorrJet->energy(); - if(useType1_) chEnF += lType1Corr*(1-jetInfo.chargedEnFrac() ); - jetInfo.setChargedEnFrac( chEnF ); - retVal.push_back(jetInfo); - break; - } - } - - //order jets per pt - std::sort( retVal.begin(), retVal.end() ); - - return retVal; -} - -std::vector PFMETProducerMVA::computePFCandidateInfo(const reco::CandidateView& pfCandidates, - const reco::Vertex* hardScatterVertex) -{ - std::vector retVal; - for ( reco::CandidateView::const_iterator pfCandidate = pfCandidates.begin(); - pfCandidate != pfCandidates.end(); ++pfCandidate ) { - double dZ = -999.; // PH: If no vertex is reconstructed in the event - // or PFCandidate has no track, set dZ to -999 - if ( hardScatterVertex ) { - const reco::PFCandidate* pfc = dynamic_cast( &(*pfCandidate) ); - if( pfc != nullptr ) { //PF candidate for RECO and PAT levels - if ( pfc->trackRef().isNonnull() ) dZ = std::abs(pfc->trackRef()->dz(hardScatterVertex->position())); - else if ( pfc->gsfTrackRef().isNonnull() ) dZ = std::abs(pfc->gsfTrackRef()->dz(hardScatterVertex->position())); - } - else { //if not, then packedCandidate for miniAOD level - const pat::PackedCandidate* pfc = dynamic_cast( &(*pfCandidate) ); - dZ = std::abs( pfc->dz( hardScatterVertex->position() ) ); - //exact dz=zero corresponds to the -999 case for pfcandidate - if(dZ==0) {dZ=-999;} - } - } - reco::PUSubMETCandInfo pfCandidateInfo; - pfCandidateInfo.setP4( pfCandidate->p4() ); - pfCandidateInfo.setDZ( dZ ); - retVal.push_back(pfCandidateInfo); - } - return retVal; -} - -std::vector PFMETProducerMVA::computeVertexInfo(const reco::VertexCollection& vertices) -{ - std::vector retVal; - for ( reco::VertexCollection::const_iterator vertex = vertices.begin(); - vertex != vertices.end(); ++vertex ) { - if(std::abs(vertex->z()) > 24.) continue; - if(vertex->ndof() < 4.) continue; - if(vertex->position().Rho() > 2.) continue; - retVal.push_back(vertex->position()); - } - return retVal; -} -double PFMETProducerMVA::chargedEnFrac(const reco::Candidate *iCand, - const reco::CandidateView& pfCandidates,const reco::Vertex* hardScatterVertex) { - if(iCand->isMuon()) { - return 1; - } - if(iCand->isElectron()) { - return 1.; - } - if(iCand->isPhoton() ) {return chargedFracInCone(iCand, pfCandidates,hardScatterVertex);} - double lPtTot = 0; double lPtCharged = 0; - const reco::PFTau *lPFTau = 0; - lPFTau = dynamic_cast(iCand); - if(lPFTau != nullptr) { - for (UInt_t i0 = 0; i0 < lPFTau->signalPFCands().size(); i0++) { - lPtTot += (lPFTau->signalPFCands())[i0]->pt(); - if((lPFTau->signalPFCands())[i0]->charge() == 0) continue; - lPtCharged += (lPFTau->signalPFCands())[i0]->pt(); - } - } - else { - const pat::Tau *lPatPFTau = nullptr; - lPatPFTau = dynamic_cast(iCand); - if(lPatPFTau != nullptr) { - for (UInt_t i0 = 0; i0 < lPatPFTau->signalCands().size(); i0++) { - lPtTot += (lPatPFTau->signalCands())[i0]->pt(); - if((lPatPFTau->signalCands())[i0]->charge() == 0) continue; - lPtCharged += (lPatPFTau->signalCands())[i0]->pt(); - } - } - } - if(lPtTot == 0) lPtTot = 1.; - return lPtCharged/lPtTot; -} -//Return tau id by process of elimination -bool PFMETProducerMVA::istau(const reco::Candidate *iCand) { - if(iCand->isMuon()) return false; - if(iCand->isElectron()) return false; - if(iCand->isPhoton()) return false; - return true; -} -bool PFMETProducerMVA::passPFLooseId(const PFJet *iJet) { - if(iJet->energy()== 0) return false; - if(iJet->neutralHadronEnergy()/iJet->energy() > 0.99) return false; - if(iJet->neutralEmEnergy()/iJet->energy() > 0.99) return false; - if(iJet->nConstituents() < 2) return false; - if(iJet->chargedHadronEnergy()/iJet->energy() <= 0 && std::abs(iJet->eta()) < 2.4 ) return false; - if(iJet->chargedEmEnergy()/iJet->energy() > 0.99 && std::abs(iJet->eta()) < 2.4 ) return false; - if(iJet->chargedMultiplicity() < 1 && std::abs(iJet->eta()) < 2.4 ) return false; - return true; -} - -double PFMETProducerMVA::chargedFracInCone(const reco::Candidate *iCand, - const reco::CandidateView& pfCandidates, - const reco::Vertex* hardScatterVertex,double iDRMax) -{ - double iDR2Max = iDRMax*iDRMax; - reco::Candidate::LorentzVector lVis(0,0,0,0); - for ( reco::CandidateView::const_iterator pfCandidate = pfCandidates.begin(); - pfCandidate != pfCandidates.end(); ++pfCandidate ) { - if(deltaR2(iCand->p4(),pfCandidate->p4()) > iDR2Max) continue; - double dZ = -999.; // PH: If no vertex is reconstructed in the event - // or PFCandidate has no track, set dZ to -999 - if ( hardScatterVertex ) { - const reco::PFCandidate* pfc = dynamic_cast( (&(*pfCandidate)) ); - if( pfc != nullptr ) { //PF candidate for RECO and PAT levels - if ( pfc->trackRef().isNonnull() ) dZ = std::abs(pfc->trackRef()->dz(hardScatterVertex->position())); - else if ( pfc->gsfTrackRef().isNonnull() ) dZ = std::abs(pfc->gsfTrackRef()->dz(hardScatterVertex->position())); - } - else { //if not, then packedCandidate for miniAOD level - const pat::PackedCandidate* pfc = dynamic_cast( &(*pfCandidate) ); - dZ = std::abs( pfc->dz( hardScatterVertex->position() ) ); - } - } - if(std::abs(dZ) > 0.1) continue; - lVis += pfCandidate->p4(); - } - return lVis.pt()/iCand->pt(); -} - -#include "FWCore/Framework/interface/MakerMacros.h" - -DEFINE_FWK_MODULE(PFMETProducerMVA); diff --git a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h b/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h deleted file mode 100644 index 4e7d9cc39957f..0000000000000 --- a/RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h +++ /dev/null @@ -1,105 +0,0 @@ -#ifndef RecoMET_METPUSubtraction_PFMETProducerMVA_h -#define RecoMET_METPUSubtraction_PFMETProducerMVA_h - -/** \class PFMETProducerMVA - * - * Produce PFMET objects computed by MVA - * - * \authors Phil Harris, CERN - * Christian Veelken, LLR - * - */ - -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/stream/EDProducer.h" -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/EventSetup.h" -#include "FWCore/ParameterSet/interface/ParameterSet.h" - - -#include "DataFormats/Candidate/interface/Candidate.h" -#include "DataFormats/Candidate/interface/CandidateFwd.h" -#include "DataFormats/JetReco/interface/PFJet.h" -#include "DataFormats/JetReco/interface/PFJetCollection.h" -#include "DataFormats/Math/interface/deltaR.h" -#include "DataFormats/METReco/interface/CommonMETData.h" -#include "DataFormats/METReco/interface/PFMET.h" -#include "DataFormats/METReco/interface/PFMETCollection.h" -#include "DataFormats/MuonReco/interface/Muon.h" -#include "DataFormats/PatCandidates/interface/PackedCandidate.h" -#include "DataFormats/PatCandidates/interface/Tau.h" -#include "DataFormats/VertexReco/interface/Vertex.h" - -#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h" - -#include "RecoMET/METAlgorithms/interface/METAlgo.h" -#include "RecoMET/METAlgorithms/interface/PFSpecificAlgo.h" -#include "RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h" -#include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h" - -#include "RecoJets/JetProducers/interface/PileupJetIdAlgo.h" -#include - -namespace reco -{ - class PFMETProducerMVA : public edm::stream::EDProducer<> - { - public: - - PFMETProducerMVA(const edm::ParameterSet&); - ~PFMETProducerMVA(); - - private: - - void produce(edm::Event&, const edm::EventSetup&); - - // auxiliary functions - std::vector computeLeptonInfo(const std::vector >& srcLeptons_, - const reco::CandidateView& pfCandidates, - const reco::Vertex* hardScatterVertex, - int& lId, bool& lHasPhotons, edm::Event & iEvent); - - std::vector computeJetInfo(const reco::PFJetCollection&, const edm::Handle&, - const edm::ValueMap&, const reco::VertexCollection&, - const reco::Vertex*, const reco::JetCorrector &iCorr, - edm::Event & iEvent,const edm::EventSetup &iSetup, - std::vector &iLeptons, - std::vector &iCands); - - std::vector computePFCandidateInfo(const reco::CandidateView&, const reco::Vertex*); - std::vector computeVertexInfo(const reco::VertexCollection&); - double chargedEnFrac(const reco::Candidate *iCand,const reco::CandidateView& pfCandidates,const reco::Vertex* hardScatterVertex); - - bool passPFLooseId(const reco::PFJet *iJet); - bool istau (const reco::Candidate *iCand); - double chargedFracInCone(const reco::Candidate *iCand,const reco::CandidateView& pfCandidates,const reco::Vertex* hardScatterVertex,double iDRMax=0.2); - - // configuration parameter - edm::EDGetTokenT srcCorrJets_; - edm::EDGetTokenT srcUncorrJets_; - edm::EDGetTokenT > srcJetIds_; - //edm::EDGetTokenT srcPFCandidates_; - edm::EDGetTokenT > srcPFCandidatesView_; - edm::EDGetTokenT srcVertices_; - edm::EDGetTokenT mJetCorrector_; - typedef std::vector vInputTag; - std::vector > srcLeptons_; - int minNumLeptons_; // CV: option to skip MVA MET computation in case there are less than specified number of leptons in the event - - bool useType1_; - - double globalThreshold_; - - double minCorrJetPt_; - - METAlgo metAlgo_; - PFSpecificAlgo pfMEtSpecificAlgo_; - PFMETAlgorithmMVA mvaMEtAlgo_; - bool mvaMEtAlgo_isInitialized_; - // PileupJetIdAlgo mvaJetIdAlgo_; - - int verbosity_; - }; -} - -#endif diff --git a/RecoMET/METPUSubtraction/plugins/neutralCandidatePUIDJets.cc b/RecoMET/METPUSubtraction/plugins/neutralCandidatePUIDJets.cc new file mode 100644 index 0000000000000..f1f9ab729781d --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/neutralCandidatePUIDJets.cc @@ -0,0 +1,226 @@ +#ifndef JMEValidator_neutralCandidatePUIDJets_h +#define JMEValidator_neutralCandidatePUIDJets_h + + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "PhysicsTools/SelectorUtils/interface/PFJetIDSelectionFunctor.h" +#include "DataFormats/JetReco/interface/PileupJetIdentifier.h" +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include + + +class neutralCandidatePUIDJets : public edm::stream::EDProducer<> { + + public: + + neutralCandidatePUIDJets(const edm::ParameterSet&); + ~neutralCandidatePUIDJets(){}; + + void produce(edm::Event&, const edm::EventSetup&); + + void beginJob(); + void endJob(){}; + + private: + + edm::InputTag srcJets_; + edm::InputTag srcCandidates_; + std::string neutralParticlesPVJets_ ; + std::string neutralParticlesPUJets_ ; + std::string neutralParticlesUnclustered_ ; + std::string PUJets_; + std::string PVJets_; + std::string jetPUIDWP_; + + std::string jetPUIDMapLabel_ ; + + PileupJetIdentifier::Id jetIdSelection_; + + edm::EDGetTokenT srcJetsToken_; + edm::EDGetTokenT srcCandidatesToken_; + + static std::string jetPUIDNameLabel_ ; + static bool stringInJetCollection_ ; + + float jetPUIDCut_ [4][3]; + +}; +#endif + +bool neutralCandidatePUIDJets::stringInJetCollection_ = false; +std::string neutralCandidatePUIDJets::jetPUIDNameLabel_ = ""; + +neutralCandidatePUIDJets::neutralCandidatePUIDJets(const edm::ParameterSet& iConfig){ + + // new PU Jet ID cuts medium + // https://indico.cern.ch/event/502737/contribution/3/attachments/1234291/1816811/PileupJetID_76X.pdf + jetPUIDCut_[0][0] = -0.61; jetPUIDCut_[0][1] = -0.52; jetPUIDCut_[0][2] = -0.40; jetPUIDCut_[0][3] = -0.36; + jetPUIDCut_[1][0] = -0.61; jetPUIDCut_[1][1] = -0.52; jetPUIDCut_[1][2] = -0.40; jetPUIDCut_[1][3] = -0.36; + jetPUIDCut_[2][0] = -0.61; jetPUIDCut_[2][1] = -0.52; jetPUIDCut_[2][2] = -0.40; jetPUIDCut_[2][3] = -0.36; + jetPUIDCut_[3][0] = -0.20; jetPUIDCut_[3][1] = -0.37; jetPUIDCut_[3][2] = -0.22; jetPUIDCut_[3][3] = -0.17; + jetPUIDCut_[4][0] = -0.20; jetPUIDCut_[4][1] = -0.37; jetPUIDCut_[4][2] = -0.22; jetPUIDCut_[4][3] = -0.17; + + + srcJets_ = iConfig.getParameter("srcJets"); + srcCandidates_ = iConfig.getParameter("srcCandidates"); + srcJetsToken_ = consumes(srcJets_); + srcCandidatesToken_ = consumes(srcCandidates_); + neutralParticlesPVJets_ = iConfig.getParameter("neutralParticlesPVJetsLabel"); + neutralParticlesPUJets_ = iConfig.getParameter("neutralParticlesPUJetsLabel"); + neutralParticlesUnclustered_ = iConfig.getParameter("neutralParticlesUnclusteredLabel"); + PUJets_ = "PUJets"; + PVJets_ = "PVJets"; + + jetPUIDWP_ = iConfig.getParameter("jetPUDIWP"); + if(jetPUIDWP_ != "tight" and jetPUIDWP_ != "medium" and jetPUIDWP_ != "loose" and jetPUIDWP_ != "user") + throw cms::Exception("Configuration") <<" [neutralCandidatePUIDJets] wrong label for jetPUID working point "; + + if(jetPUIDWP_ == "tight") + jetIdSelection_ = PileupJetIdentifier::kTight; + else if(jetPUIDWP_ == "medium") + jetIdSelection_ = PileupJetIdentifier::kMedium; + else if(jetPUIDWP_ == "loose") + jetIdSelection_ = PileupJetIdentifier::kTight; + + + jetPUIDMapLabel_ = iConfig.getParameter("jetPUIDMapLabel"); + + produces >(neutralParticlesPVJets_); + produces >(neutralParticlesPUJets_); + produces >(neutralParticlesUnclustered_); + produces(PUJets_); + produces(PVJets_); + +} + +void neutralCandidatePUIDJets::produce(edm::Event& iEvent, const edm::EventSetup& iSetup){ + + edm::Handle jetCollection; + iEvent.getByToken(srcJetsToken_,jetCollection); + + edm::Handle candCollection; + iEvent.getByToken(srcCandidatesToken_,candCollection); + + std::auto_ptr > neutralParticlesPVJets(new edm::PtrVector); + std::auto_ptr > neutralParticlesPUJets(new edm::PtrVector); + std::auto_ptr > neutralParticlesUnclustered(new edm::PtrVector); + + std::auto_ptr PUJets(new pat::JetCollection); + std::auto_ptr PVJets(new pat::JetCollection); + + // loop on jets + for(auto jet : *jetCollection){ + // look if the value is embedded in pat jets + if(!stringInJetCollection_) + { + if(jetPUIDWP_ != "user") + { + for(size_t iJet = 0; iJet < jet.userIntNames().size(); iJet++) + { + if(jet.userIntNames().at(iJet).find(jetPUIDMapLabel_) != std::string::npos) + { + stringInJetCollection_ = true; + jetPUIDNameLabel_ = jet.userIntNames().at(iJet); + } + } + if(stringInJetCollection_ == false) + throw cms::Exception("neutralCandidatePUIDJets")<<" user int related to jetPUID not found "; + } + else + { + for(size_t iJet = 0; iJet < jet.userFloatNames().size(); iJet++) + { + if(jet.userFloatNames().at(iJet).find(jetPUIDMapLabel_) != std::string::npos) + { + stringInJetCollection_ = true; + jetPUIDNameLabel_ = jet.userFloatNames().at(iJet); + } + } + if(stringInJetCollection_ == false) + throw cms::Exception("neutralCandidatePUIDJets")<<" user float related to jetPUID not found "; + } + } + + // evaluate PU JET ID + bool isPassingPUID = false; + if(jetPUIDWP_ != "user") + { + isPassingPUID = PileupJetIdentifier::passJetId(jet.userInt(jetPUIDNameLabel_), jetIdSelection_); + } + else + { + + int ptBin = 0; + int etaBin = 0; + if ( jet.pt() >= 10. && jet.pt() < 20. ) ptBin = 1; + if ( jet.pt() >= 20. && jet.pt() < 30. ) ptBin = 2; + if ( jet.pt() >= 30. && jet.pt() < 50. ) ptBin = 3; + if ( jet.pt() >= 50. ) ptBin = 4; + if ( std::abs(jet.eta()) >= 2.5 && std::abs(jet.eta()) < 2.75) etaBin = 1; + if ( std::abs(jet.eta()) >= 2.75 && std::abs(jet.eta()) < 3.0 ) etaBin = 2; + if ( std::abs(jet.eta()) >= 3.0 && std::abs(jet.eta()) < 5.0 ) etaBin = 3; + isPassingPUID = jet.userFloat(jetPUIDNameLabel_) > jetPUIDCut_[ptBin][etaBin] ? true : false ; + + } + if(isPassingPUID) + PVJets->push_back(jet); + else + PUJets->push_back(jet); + // loop on constituents + for( auto particle : jet.getJetConstituents()) + { + if(particle->charge() != 0) continue; + if(isPassingPUID) + neutralParticlesPVJets->push_back(particle); + else + neutralParticlesPUJets->push_back(particle); + } + } + + + // loop on pfParticles to determine if unclustered Neutral + size_t indexColl = 0; + for(reco::CandidateView::const_iterator itCand = candCollection->begin(); itCand != candCollection->end(); itCand++) + { + assert(itCand->charge() == 0); + + bool clustered = false; + for(edm::PtrVector::const_iterator iParticle = neutralParticlesPUJets->begin(); iParticle != neutralParticlesPUJets->end(); iParticle++) + { + if( itCand->p4() == iParticle->get()->p4()) + { + clustered = true; + break; + } + } + for(edm::PtrVector::const_iterator iParticle = neutralParticlesPVJets->begin(); iParticle != neutralParticlesPVJets->end(); iParticle++) + { + if( itCand->p4() == iParticle->get()->p4()) + { + clustered = true; + break; + } + } + + if(!clustered) + { + neutralParticlesUnclustered->push_back(edm::Ptr(candCollection, indexColl)); + } + indexColl++; + } + + + iEvent.put(neutralParticlesPVJets,neutralParticlesPVJets_); + iEvent.put(neutralParticlesPUJets,neutralParticlesPUJets_); + iEvent.put(neutralParticlesUnclustered,neutralParticlesUnclustered_); + iEvent.put(PUJets,PUJets_); + iEvent.put(PVJets,PVJets_); + +} + +DEFINE_FWK_MODULE(neutralCandidatePUIDJets); diff --git a/RecoMET/METPUSubtraction/plugins/patElectronIDIsoSelector.cc b/RecoMET/METPUSubtraction/plugins/patElectronIDIsoSelector.cc new file mode 100644 index 0000000000000..44814b6565bac --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/patElectronIDIsoSelector.cc @@ -0,0 +1,330 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDFilter.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/CandidateFwd.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" + +#include +#include + +using namespace std; +using namespace edm; +using namespace pat; + + +//////////////////////////////////////////////////////////////////////////////// +// class definition +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +class patElectronIDIsoSelector : public edm::stream::EDProducer<> { + +public: + // construction/destruction + patElectronIDIsoSelector(const edm::ParameterSet& iConfig); + ~patElectronIDIsoSelector() {;} + + // member functions + void produce(edm::Event& iEvent,const edm::EventSetup& iSetup); + void endJob(); + +private: + // input electron collection to be considered + edm::InputTag src_; + // rho collection in case of effective area + edm::InputTag rho_; + // vertex collection + edm::InputTag vertex_; + + // value map with re-calcuated isolations from iso-deposits .. different wrt to the ones already inside the muon object + edm::InputTag charged_hadron_iso_; + edm::InputTag neutral_hadron_iso_; + edm::InputTag photon_iso_; + + // input tag with the value map for the electron id to be applied + edm::InputTag electron_id_; + + // isolation, pt and eta cut to be applied + double relativeIsolationCut_; + double ptCut_; + double etaCut_; + + // type of relative isolation cut + std::string typeIso_; + + // tokens + edm::EDGetTokenT srcToken_ ; + edm::EDGetTokenT rhoToken_ ; + edm::EDGetTokenT vertexToken_ ; + + edm::EDGetTokenT > electron_idToken_; + + edm::EDGetTokenT > charged_hadron_isoToken_ ; + edm::EDGetTokenT > neutral_hadron_isoToken_ ; + edm::EDGetTokenT > photon_isoToken_ ; + +}; + + + + +//////////////////////////////////////////////////////////////////////////////// +// construction/destruction +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +patElectronIDIsoSelector::patElectronIDIsoSelector(const edm::ParameterSet& iConfig){ + + + if(iConfig.existsAs("src")){ + src_ = iConfig.getParameter("src"); + if(src_ == edm::InputTag("")) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] empty electron tag is given --> please check \n"; + } + else throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] no electron tag is given \n"; + + if(iConfig.existsAs("vertex")){ + vertex_ = iConfig.getParameter("vertex"); + if(vertex_ == edm::InputTag("")) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] empty vertex tag is given --> please check \n"; + } + else throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] no vertex tag is given \n"; + + if(iConfig.existsAs("rho")) + rho_ = iConfig.getParameter("rho"); + + + if(iConfig.existsAs("charged_hadron_iso")) + charged_hadron_iso_ = iConfig.getParameter("charged_hadron_iso"); + + + if(iConfig.existsAs("neutral_hadron_iso")) + neutral_hadron_iso_ = iConfig.getParameter("neutral_hadron_iso"); + + + if(iConfig.existsAs("photon_iso")) + photon_iso_ = iConfig.getParameter("photon_iso"); + + if(iConfig.existsAs("electron_id")){ + electron_id_ = iConfig.getParameter("electron_id"); + if(electron_id_ == edm::InputTag("")) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] no electron id value map provided --> please check \n"; + } + else throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] no electron id value map is given \n"; + + if(iConfig.existsAs("relativeIsolationCut")) + relativeIsolationCut_ = iConfig.getParameter("relativeIsolationCut"); + + if(iConfig.existsAs("ptCut")) + ptCut_ = iConfig.getParameter("ptCut"); + + if(iConfig.existsAs("etaCut")) + etaCut_ = iConfig.getParameter("etaCut"); + + if(iConfig.existsAs("typeIso")) + typeIso_ = iConfig.getParameter("typeIso"); + + // tokens + if(!(src_ == edm::InputTag(""))) + srcToken_ = consumes(src_); + + if(!(rho_ == edm::InputTag(""))) + rhoToken_ = consumes(rho_); + + if(!(vertex_ == edm::InputTag(""))) + vertexToken_ = consumes(vertex_); + + if(!(charged_hadron_iso_ == edm::InputTag(""))) + charged_hadron_isoToken_ = consumes >(charged_hadron_iso_); + + if(!(neutral_hadron_iso_ == edm::InputTag(""))) + neutral_hadron_isoToken_ = consumes >(neutral_hadron_iso_); + + if(!(photon_iso_ == edm::InputTag(""))) + photon_isoToken_ = consumes >(photon_iso_); + + if(!(electron_id_ == edm::InputTag(""))) + electron_idToken_ = consumes >(electron_id_); + + + produces(); +} + + +//////////////////////////////////////////////////////////////////////////////// +// implementation of member functions +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +void patElectronIDIsoSelector::produce(edm::Event& iEvent,const edm::EventSetup& iSetup){ + + edm::Handle ElectronCollectionHandle; + if(!(src_ == edm::InputTag(""))){ + iEvent.getByToken(srcToken_,ElectronCollectionHandle); + if(ElectronCollectionHandle.failedToGet()){ + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] failed to get electron collection \n"; + } + } + + edm::Handle VertexCollectionHandle; + if(!(vertex_ == edm::InputTag(""))) { + iEvent.getByToken(vertexToken_,VertexCollectionHandle); + if(VertexCollectionHandle.failedToGet()) { + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] failed to get vertex collection \n"; + } + } + + + edm::Handle rhoHandle; + double rho = -1; + if(!(rho_ == edm::InputTag(""))) { + iEvent.getByToken(rhoToken_,rhoHandle); + if(!rhoHandle.failedToGet()){ + rho = *rhoHandle; + } + } + + // take the isolation maps + bool isGoodChargeIsoMap = false; + edm::Handle> chargeHadronIsoHandle; + if(!(charged_hadron_iso_ == edm::InputTag(""))){ + iEvent.getByToken(charged_hadron_isoToken_,chargeHadronIsoHandle); + if(!chargeHadronIsoHandle.failedToGet()) + isGoodChargeIsoMap = true; + } + + bool isGoodNeutralIsoMap = false; + edm::Handle> neutralHadronIsoHandle; + if(!(neutral_hadron_iso_ == edm::InputTag(""))){ + iEvent.getByToken(neutral_hadron_isoToken_,neutralHadronIsoHandle); + if(!neutralHadronIsoHandle.failedToGet()) + isGoodNeutralIsoMap = true; + } + + bool isGoodPhotonIsoMap = false; + edm::Handle> photonIsoHandle; + if(!(photon_iso_ == edm::InputTag(""))){ + iEvent.getByToken(photon_isoToken_,photonIsoHandle); + if(!photonIsoHandle.failedToGet()) + isGoodPhotonIsoMap = true; + } + + + // take the value map for ele ID + bool isGoodEleValueMap = false; + edm::Handle> electronIDHandle; + if(!(electron_id_ == edm::InputTag(""))){ + iEvent.getByToken(electron_idToken_,electronIDHandle); + if(!electronIDHandle.failedToGet()) + isGoodEleValueMap = true; + } + + if(isGoodEleValueMap == false) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] failed to get the value map for ele id .. please check \n"; + + if((typeIso_ == "dBetaWeight" and not isGoodNeutralIsoMap) or (typeIso_ == 3 and not isGoodPhotonIsoMap)) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] empty maps for neutrals with PFWeight Dbeta correction .. please check \n"; + + if((typeIso_ == "puppi" and not isGoodNeutralIsoMap) or (typeIso_ == "puppi" and not isGoodPhotonIsoMap) or (typeIso_ == "puppi" and not isGoodChargeIsoMap)) + throw cms::Exception("Configuration")<<"[patElectronIDIsoSelector] empty maps for puppi isolation correction .. please check \n"; + + // Loop on muon collection + size_t electronIndex = 0; + pat::ElectronCollection::const_iterator itElectron = ElectronCollectionHandle->begin(); + + auto_ptr electronColl(new vector); + + for( ; itElectron != ElectronCollectionHandle->end(); itElectron++, electronIndex++){ + + // apply the pt/eta cut on the best track + double electronPt = 0; + reco::GsfTrack electronTrack; + if(itElectron->gsfTrack().isNonnull()){ + electronTrack = *(itElectron->gsfTrack().get()); + electronPt = electronTrack.pt(); + if(electronPt < ptCut_ or fabs(electronTrack.eta()) > etaCut_) continue; + } + else{ + electronPt = itElectron->pt(); + if (itElectron->pt() < ptCut_ or fabs(itElectron->eta()) > etaCut_) continue; + } + + // apply the id starting from the value map + const pat::ElectronRef electronRef(ElectronCollectionHandle,electronIndex); + + if(not (*electronIDHandle)[electronRef]) continue; + + // isolation + + float chargeIso = 0; + if(isGoodChargeIsoMap) + chargeIso = (*chargeHadronIsoHandle)[electronRef]; + else{ + chargeIso = itElectron->userIsolation("PfChargedHadronIso"); + } + + float neutralIso = 0; + if(isGoodNeutralIsoMap) + neutralIso = (*neutralHadronIsoHandle)[electronRef]; + else{ + neutralIso = itElectron->userIsolation("PfNeutralHadronIso"); + } + + float photonIso = 0; + if(isGoodPhotonIsoMap) + photonIso = (*photonIsoHandle)[electronRef]; + else{ + photonIso = itElectron->userIsolation("PfGammaIso"); + } + + if(typeIso_ == "" or typeIso_ == "dBetaWeight" or typeIso_ == "puppi"){ + if( (chargeIso+neutralIso+photonIso)/electronPt >= relativeIsolationCut_) continue; + } + else if(typeIso_ == "dBeta"){ + if( (chargeIso+max(neutralIso+photonIso-0.5*itElectron->userIsolation("PfPUChargedHadronIso"),0.))/electronPt >= relativeIsolationCut_) continue; + } + else if(typeIso_ == ""){ + if( (chargeIso+max(neutralIso+photonIso-rho*TMath::Pi()*0.4*0.4,0.))/electronPt >= relativeIsolationCut_) continue; + } + + pat::Electron newElectron = *(itElectron->clone()); + if(typeIso_ == 3){ + newElectron.addUserFloat("neutralHadronIsoPFWeight03",neutralIso); + newElectron.addUserFloat("photonIsoPFWeight03",photonIso); + } + else if(typeIso_ == 4){ + newElectron.addUserFloat("neutralHadronIsoPuppi03",neutralIso); + newElectron.addUserFloat("photonIsoPuppi03",photonIso); + newElectron.addUserFloat("chargedHadronIsoPuppi03",chargeIso); + } + + electronColl->push_back(newElectron); + } + + iEvent.put(electronColl); + + +} + + +//______________________________________________________________________________ +void patElectronIDIsoSelector::endJob(){} + + +//////////////////////////////////////////////////////////////////////////////// +// plugin definition +//////////////////////////////////////////////////////////////////////////////// + +DEFINE_FWK_MODULE(patElectronIDIsoSelector); diff --git a/RecoMET/METPUSubtraction/plugins/patMuonIDIsoSelector.cc b/RecoMET/METPUSubtraction/plugins/patMuonIDIsoSelector.cc new file mode 100644 index 0000000000000..e13dde9ab8654 --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/patMuonIDIsoSelector.cc @@ -0,0 +1,359 @@ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDFilter.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/CandidateFwd.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" + +#include +#include + +using namespace std; +using namespace edm; +using namespace pat; + + +//////////////////////////////////////////////////////////////////////////////// +// class definition +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +class patMuonIDIsoSelector : public edm::stream::EDProducer<> { + +public: + // construction/destruction + patMuonIDIsoSelector(const edm::ParameterSet& iConfig); + ~patMuonIDIsoSelector() {;} + + // member functions + void produce(edm::Event& iEvent,const edm::EventSetup& iSetup); + void endJob(); + +private: + // input muon collection to be considered + edm::InputTag src_; + // rho collection in case of effective area + edm::InputTag rho_; + // vertex collection + edm::InputTag vertex_; + + // value map with re-calcuated isolations from iso-deposits .. different wrt to the ones already inside the muon object + edm::InputTag charged_hadron_iso_; + edm::InputTag neutral_hadron_iso_; + edm::InputTag photon_iso_; + + // isolation, pt and eta cut to be applied + double relativeIsolationCut_; + double ptCut_; + double etaCut_; + double ptThresholdForHighPt_; + + // tipe of muon id : tightID, mediumID, looseID, softID and highPt + std::string typeID_; + + // type of relative isolation cut + std::string typeIso_; + + // tokens + edm::EDGetTokenT srcToken_ ; + edm::EDGetTokenT rhoToken_ ; + edm::EDGetTokenT vertexToken_ ; + + edm::EDGetTokenT > charged_hadron_isoToken_ ; + edm::EDGetTokenT > neutral_hadron_isoToken_ ; + edm::EDGetTokenT > photon_isoToken_ ; + +}; + + + + +//////////////////////////////////////////////////////////////////////////////// +// construction/destruction +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +patMuonIDIsoSelector::patMuonIDIsoSelector(const edm::ParameterSet& iConfig){ + + + if(iConfig.existsAs("src")){ + src_ = iConfig.getParameter("src"); + if(src_ == edm::InputTag("")) + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] empty muon tag is given --> at least one \n"; + } + else throw cms::Exception("Configuration")<<"[SelectedElectronFEDListProducer] no muon tag is given \n"; + + if(iConfig.existsAs("vertex")){ + vertex_ = iConfig.getParameter("vertex"); + if(vertex_ == edm::InputTag("")) + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] empty vertex tag is given --> at least one \n"; + } + else throw cms::Exception("Configuration")<<"[SelectedElectronFEDListProducer] no vertex tag is given \n"; + + if(iConfig.existsAs("rho")) + rho_ = iConfig.getParameter("rho"); + + + if(iConfig.existsAs("charged_hadron_iso")) + charged_hadron_iso_ = iConfig.getParameter("charged_hadron_iso"); + + + if(iConfig.existsAs("neutral_hadron_iso")) + neutral_hadron_iso_ = iConfig.getParameter("neutral_hadron_iso"); + + + if(iConfig.existsAs("photon_iso")) + photon_iso_ = iConfig.getParameter("photon_iso"); + + + if(iConfig.existsAs("relativeIsolationCut")) + relativeIsolationCut_ = iConfig.getParameter("relativeIsolationCut"); + else + relativeIsolationCut_ = 0.2; + + if(iConfig.existsAs("ptCut")) + ptCut_ = iConfig.getParameter("ptCut"); + else + ptCut_ = 10.; + + if(iConfig.existsAs("etaCut")) + etaCut_ = iConfig.getParameter("etaCut"); + else + etaCut_ = 2.4; + + if(iConfig.existsAs("typeID")) + typeID_ = iConfig.getParameter("typeID"); + else + typeID_ = "TightID"; + + if(iConfig.existsAs("typeIso")) + typeIso_ = iConfig.getParameter("typeIso"); + else + typeIso_ = "dBeta"; + + if(iConfig.existsAs("ptThresholdForHighPt")) + ptThresholdForHighPt_ = iConfig.getParameter("ptThresholdForHighPt"); + else + ptThresholdForHighPt_ = 100; + + // tokens + if(!(src_ == edm::InputTag(""))) + srcToken_ = consumes(src_); + + if(!(rho_ == edm::InputTag(""))) + rhoToken_ = consumes(rho_); + + if(!(vertex_ == edm::InputTag(""))) + vertexToken_ = consumes(vertex_); + + if(!(charged_hadron_iso_ == edm::InputTag(""))) + charged_hadron_isoToken_ = consumes >(charged_hadron_iso_); + + if(!(neutral_hadron_iso_ == edm::InputTag(""))) + neutral_hadron_isoToken_ = consumes >(neutral_hadron_iso_); + + if(!(photon_iso_ == edm::InputTag(""))) + photon_isoToken_ = consumes >(photon_iso_); + + produces(); +} + + +//////////////////////////////////////////////////////////////////////////////// +// implementation of member functions +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +void patMuonIDIsoSelector::produce(edm::Event& iEvent,const edm::EventSetup& iSetup){ + + edm::Handle MuonCollectionHandle; + if(!(src_ == edm::InputTag(""))){ + iEvent.getByToken(srcToken_,MuonCollectionHandle); + if(MuonCollectionHandle.failedToGet()){ + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] failed to get muon collection \n"; + } + } + + edm::Handle VertexCollectionHandle; + if(!(vertex_ == edm::InputTag(""))) { + iEvent.getByToken(vertexToken_,VertexCollectionHandle); + if(VertexCollectionHandle.failedToGet()) { + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] failed to get vertex collection \n"; + } + } + + const reco::Vertex& primaryVertex = VertexCollectionHandle->at(0); + + edm::Handle rhoHandle; + double rho = -1; + if(!(rho_ == edm::InputTag(""))) { + iEvent.getByToken(rhoToken_,rhoHandle); + if(!rhoHandle.failedToGet()){ + rho = *rhoHandle; + } + } + + // take the isolation maps + bool isGoodChargeIsoMap = false; + edm::Handle> chargeHadronIsoHandle; + if(!(charged_hadron_iso_ == edm::InputTag(""))){ + iEvent.getByToken(charged_hadron_isoToken_,chargeHadronIsoHandle); + if(!chargeHadronIsoHandle.failedToGet()) + isGoodChargeIsoMap = true; + } + + bool isGoodNeutralIsoMap = false; + edm::Handle> neutralHadronIsoHandle; + if(!(neutral_hadron_iso_ == edm::InputTag(""))){ + iEvent.getByToken(neutral_hadron_isoToken_,neutralHadronIsoHandle); + if(!neutralHadronIsoHandle.failedToGet()) + isGoodNeutralIsoMap = true; + } + + bool isGoodPhotonIsoMap = false; + edm::Handle> photonIsoHandle; + if(!(photon_iso_ == edm::InputTag(""))){ + iEvent.getByToken(photon_isoToken_,photonIsoHandle); + if(!photonIsoHandle.failedToGet()) + isGoodPhotonIsoMap = true; + } + + if((typeIso_ == "dBetaWeight" and not isGoodNeutralIsoMap) or (typeIso_ == "dBetaWeight" and not isGoodPhotonIsoMap)) + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] empty maps for neutrals with PFWeight Dbeta correction .. please check \n"; + + if((typeIso_ == "puppi" and not isGoodNeutralIsoMap) or (typeIso_ == "puppi" and not isGoodPhotonIsoMap) or (typeIso_ == "puppi" and not isGoodChargeIsoMap)) + throw cms::Exception("Configuration")<<"[patMuonIDIsoSelector] empty maps for puppi isolation correction .. please check \n"; + + // Loop on muon collection + size_t muonIndex = 0; + pat::MuonCollection::const_iterator itMuon = MuonCollectionHandle->begin(); + + auto_ptr muonColl(new vector); + + for( ; itMuon != MuonCollectionHandle->end(); itMuon++, muonIndex++){ + + // apply the pt/eta cut on the best track + double muonPt = 0; + reco::Track muonTrack; + if(itMuon->tunePMuonBestTrack().isNonnull()){ + muonTrack = *(itMuon->tunePMuonBestTrack().get()); + muonPt = muonTrack.pt(); + if(muonPt < ptCut_ or fabs(muonTrack.eta()) > etaCut_) continue; + } + else if(itMuon->muonBestTrack().isNonnull()){ + muonTrack = *(itMuon->muonBestTrack().get()); + muonPt = muonTrack.pt(); + if(muonPt < ptCut_ or fabs(muonTrack.eta()) > etaCut_) continue; + } + else{ + muonPt = itMuon->pt(); + if (itMuon->pt() < ptCut_ or fabs(itMuon->eta()) > etaCut_) continue; + } + + // apply standard muon id + if(typeID_ == "tightID" or typeID_ == "TightID" or typeID_ == "tight" or typeID_ == "Tight"){ + if( muonPt < ptThresholdForHighPt_ ){ + if( not muon::isTightMuon(*itMuon, primaryVertex)) continue; + } + else{ + if( not muon::isTightMuon(*itMuon, primaryVertex) and not muon::isHighPtMuon(*itMuon,primaryVertex)) continue; + } + + } + else if(typeID_ == "mediumID" or typeID_ == "MediumID" or typeID_ == "medium" or typeID_ == "Medium"){ + if( muonPt < ptThresholdForHighPt_ ){ + if( not muon::isMediumMuon(*itMuon)) continue; + } + else{ + if( not muon::isMediumMuon(*itMuon) and not muon::isHighPtMuon(*itMuon,primaryVertex)) continue; + } + } + else if(typeID_ == "looseID" or typeID_ == "LooseID" or typeID_ == "loose" or typeID_ == "Loose"){ + if( muonPt < ptThresholdForHighPt_ ){ + if( not muon::isLooseMuon(*itMuon)) continue; + } + else{ + if( not muon::isLooseMuon(*itMuon) and not muon::isHighPtMuon(*itMuon,primaryVertex)) continue; + } + } + else if(typeID_ == "softID" or typeID_ == "SoftID" or typeID_ == "soft" or typeID_ == "Soft"){ + if( not muon::isSoftMuon(*itMuon,primaryVertex)) continue; + } + else if(typeID_ == "highPt" or typeID_ == "HightPt"){ + if( not muon::isHighPtMuon(*itMuon,primaryVertex)) continue; + } + + // isolation + const pat::MuonRef muonRef(MuonCollectionHandle,muonIndex); + + float chargeIso = 0; + if(isGoodChargeIsoMap) + chargeIso = (*chargeHadronIsoHandle)[muonRef]; + else{ + chargeIso = itMuon->userIsolation("PfChargedHadronIso"); + } + + float neutralIso = 0; + if(isGoodNeutralIsoMap) + neutralIso = (*neutralHadronIsoHandle)[muonRef]; + else{ + neutralIso = itMuon->userIsolation("PfNeutralHadronIso"); + } + + float photonIso = 0; + if(isGoodPhotonIsoMap) + photonIso = (*photonIsoHandle)[muonRef]; + else{ + photonIso = itMuon->userIsolation("PfGammaIso"); + } + + if(typeIso_ == "" or typeIso_ == "dBetaWeight" or typeIso_ == "puppi"){ + if( (chargeIso+neutralIso+photonIso)/muonPt >= relativeIsolationCut_) continue; + } + else if(typeIso_ == "dBeta"){ + if( (chargeIso+max(neutralIso+photonIso-0.5*itMuon->userIsolation("PfPUChargedHadronIso"),0.))/muonPt >= relativeIsolationCut_) continue; + } + else if(typeIso_ == "rhoCorr"){ + if( (chargeIso+max(neutralIso+photonIso-rho*TMath::Pi()*0.4*0.4,0.))/muonPt >= relativeIsolationCut_) continue; + } + + pat::Muon newMuon = *(itMuon->clone()); + if(typeIso_ == "dBetaWeight"){ + newMuon.addUserFloat("neutralHadronIsoPFWeight04",neutralIso); + newMuon.addUserFloat("photonIsoPFWeight04",photonIso); + } + else if(typeIso_ == "puppi"){ + newMuon.addUserFloat("neutralHadronIsoPuppi04",neutralIso); + newMuon.addUserFloat("photonIsoPuppi04",photonIso); + newMuon.addUserFloat("chargedHadronIsoPuppi04",chargeIso); + } + + muonColl->push_back(newMuon); + } + + iEvent.put(muonColl); + + +} + + +//______________________________________________________________________________ +void patMuonIDIsoSelector::endJob(){} + + +//////////////////////////////////////////////////////////////////////////////// +// plugin definition +//////////////////////////////////////////////////////////////////////////////// + +DEFINE_FWK_MODULE(patMuonIDIsoSelector); diff --git a/RecoMET/METPUSubtraction/plugins/tauDecayProducts.cc b/RecoMET/METPUSubtraction/plugins/tauDecayProducts.cc new file mode 100644 index 0000000000000..9ed20132fb8f1 --- /dev/null +++ b/RecoMET/METPUSubtraction/plugins/tauDecayProducts.cc @@ -0,0 +1,111 @@ + + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +//#include "FWCore/ParameterSet/interface/ParameterSet.h" +//#include "FWCore/MessageLogger/interface/MessageLogger.h" +//#include "FWCore/ServiceRegistry/interface/Service.h" + +//#include "DataFormats/Common/interface/Handle.h" +//#include "DataFormats/Candidate/interface/Candidate.h" +//#include "DataFormats/Candidate/interface/CandidateFwd.h" +//#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" +#include "DataFormats/PatCandidates/interface/CompositeCandidate.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "DataFormats/PatCandidates/interface/Tau.h" +#include +#include + +using namespace std; +using namespace edm; +using namespace reco; + + +//////////////////////////////////////////////////////////////////////////////// +// class definition +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +class tauDecayProducts : public edm::EDProducer +{ +public: + // construction/destruction + tauDecayProducts(const edm::ParameterSet& iConfig); + ~tauDecayProducts() {;} + + // member functions + void produce(edm::Event& iEvent,const edm::EventSetup& iSetup); + void endJob(); + +private: + // member data + edm::EDGetTokenT src_; + + std::string moduleName_; + +}; + + + + +//////////////////////////////////////////////////////////////////////////////// +// construction/destruction +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +tauDecayProducts::tauDecayProducts(const edm::ParameterSet& iConfig) + : moduleName_(iConfig.getParameter("@module_label")) +{ + src_ = consumes(iConfig.getParameter("src")); + produces(); +} + + +//////////////////////////////////////////////////////////////////////////////// +// implementation of member functions +//////////////////////////////////////////////////////////////////////////////// + +//______________________________________________________________________________ +void tauDecayProducts::produce(edm::Event& iEvent,const edm::EventSetup& iSetup) +{ + edm::Handle tauCollectionHandle; + iEvent.getByToken(src_, tauCollectionHandle); + const pat::TauCollection tauCollection = *(tauCollectionHandle.product()); + + + + auto_ptr recoCands(new vector); + + for (auto it : tauCollection) + { + for(auto cand : it.signalPFChargedHadrCands() ) + { + recoCands->push_back(cand); + } + for(auto cand : it.signalPFGammaCands() ) + { + recoCands->push_back(cand); + } + } + + //iEvent.put(recoCands,"convertedPackedPFCandidates"); + iEvent.put(recoCands); +} + + +//______________________________________________________________________________ +void tauDecayProducts::endJob() +{ +} + + +//////////////////////////////////////////////////////////////////////////////// +// plugin definition +//////////////////////////////////////////////////////////////////////////////// + +DEFINE_FWK_MODULE(tauDecayProducts); diff --git a/RecoMET/METPUSubtraction/python/LeptonSelectionTools_cff.py b/RecoMET/METPUSubtraction/python/LeptonSelectionTools_cff.py new file mode 100644 index 0000000000000..e5398b63d631e --- /dev/null +++ b/RecoMET/METPUSubtraction/python/LeptonSelectionTools_cff.py @@ -0,0 +1,215 @@ +import FWCore.ParameterSet.Config as cms + +##### to apply muon ID +def applyMuonID(process, + label = "Tight", ##add to the process module + src = 'slimmedMuons', ## input collection + iso_map_charged_hadron = '', ## isolation value map from charged hadrons + iso_map_neutral_hadron = '', ## isolation value map from neutral hadrons + iso_map_photon = '', ## isolation value map from photons + rho = 'fixedGridRhoFastjetAll', + vertex = 'offlineSlimmedPrimaryVertices', + ptVal = 10., + etaVal = 2.4, + typeIso = "dBeta", + relativeIsolationCutVal = 0.12 + ): + + setattr(process,src+label, cms.EDProducer("patMuonIDIsoSelector", + src = cms.InputTag(src), + rho = cms.InputTag(rho), + vertex = cms.InputTag(vertex), + charged_hadron_iso = cms.InputTag(iso_map_charged_hadron), + neutral_hadron_iso = cms.InputTag(iso_map_neutral_hadron), + photon_iso = cms.InputTag(iso_map_photon), + typeID = cms.string(label), + ptCut = cms.double(ptVal), + etaCut = cms.double(etaVal), + typeIso = cms.string(typeIso), + relativeIsolationCut = cms.double(relativeIsolationCutVal) + )) + +##### to apply electron ID +def applyElectronID(process, + label = "Tight", ##add to the process module + src = 'slimmedElectrons', ## input collection + iso_map_charged_hadron = '', ## isolation value map from charged hadrons + iso_map_neutral_hadron = '', ## isolation value map from neutral hadrons + iso_map_photon = '', ## isolation value map from photons + electron_id_map = '', ## value map for electron id + rho = 'fixedGridRhoFastjetAll', + vertex = 'offlineSlimmedPrimaryVertices', + ptVal = 10., ## pt Cut + etaVal = 2.4, ## eta Cut + typeIso = "rhoCorr", + relativeIsolationCutVal = 0.13 + ): + + + setattr(process,src+label, cms.EDProducer("patElectronIDIsoSelector", + src = cms.InputTag(src), + typeID = cms.string(label), ## string to identify the type of the electron ID + rho = cms.InputTag(rho), + vertex = cms.InputTag(vertex), + charged_hadron_iso = cms.InputTag(iso_map_charged_hadron), + neutral_hadron_iso = cms.InputTag(iso_map_neutral_hadron), + photon_iso = cms.InputTag(iso_map_photon), + electron_id = cms.InputTag(electron_id_map), + ptCut = cms.double(ptVal), + etaCut = cms.double(etaVal), + typeIso = cms.string(typeIso), + relativeIsolationCut = cms.double(relativeIsolationCutVal) + )) + + +##### to apply tau ID + lepton cleaning +def applyTauID( process, + label = "Loose", + src = "slimmedTaus", ## input collection + ptCut = 20., ## ptCut + etaCut = 2.4, ## eta Cut + muonCollection = "", ## muon collection for cleaning + electronCollection = "", ## electron collection for cleaning + dRCleaning = 0.3 ## dR for cleaning + ): + + + if label == "Loose" or label == "loose" : + + ## apply tau loose ID + setattr(process,src+label,cms.EDFilter("PATTauRefSelector", + src = cms.InputTag(src), + filter = cms.bool(False), + cut = cms.string('pt > %f & abs(eta) < %f & tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits") < 1.5 & tauID("againstMuonLoose3") > 0.5 & tauID("againstElectronLooseMVA6") > 0.5 '%(ptCut,etaCut)) + )); + + + ## medium tau id + elif label == "Medium" or label == "medium" : + + setattr(process,src+label, cms.EDFilter("PATTauRefSelector", + src = cms.InputTag(src), + filter = cms.bool(False), + cut = cms.string('pt > %f & abs(eta) < %f & tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits") < 1.5 & tauID("againstMuonMedium3") > 0.5 & tauID("againstElectronMediumMVA6") > 0.5 '%(ptCut,etaCut)), + )); + ## tau tight ID + elif label == "Tight" or label == "tight" : + + setattr(process,src+label, cms.EDFilter("PATTauRefSelector", + src = cms.InputTag(src), + filter = cms.bool(False), + cut = cms.string('pt > %f & abs(eta) < %f & tauID("byTightCombinedIsolationDeltaBetaCorr3Hits") < 1.5 & tauID("againstMuonTight3") > 0.5 & tauID("againstElectronTightMVA6") > 0.5 '%(ptCut,etaCut)) + )); + + else : + print "not known string for tau ID .. please check!! " + exit(1); + + + ## cleaning taus from other leptons + tausNotOverlappingWithLeptons = cms.EDProducer("PATTauCleaner", + src = cms.InputTag(src+label), + checkOverlaps = cms.PSet(), + finalCut = cms.string(''), + preselection = cms.string('') + ) + + + if muonCollection != '': + setattr(tausNotOverlappingWithLeptons.checkOverlaps,"muons",cms.PSet( src = cms.InputTag(muonCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + if electronCollection != '': + setattr(tausNotOverlappingWithLeptons.checkOverlaps,"electrons",cms.PSet( src = cms.InputTag(electronCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + + ## copy in the process + setattr(process,src+label+"Cleaned",tausNotOverlappingWithLeptons.clone()); + + +## to clean hets from selected leptons +def cleanJetsFromLeptons (process, + label = "Cleaned", ## label to be added to output collection + jetCollection = '', ## input jet collection + muonCollection = '', ## input muon collection + electronCollection = '', ## input electron collection + tauCollection = '', ## input tau collection + jetPtCut = 30., ## pt threshold on jets + jetEtaCut = 5.0, ## eta cut + dRCleaning = 0.3 ## cleaning cone + ): + + from PhysicsTools.PatAlgos.cleaningLayer1.jetCleaner_cfi import cleanPatJets + + jetsNotOverlappingWithLeptons = cms.EDProducer("PATJetCleaner", + src = cms.InputTag(jetCollection), + preselection = cms.string(''), + checkOverlaps = cms.PSet(), + finalCut = cms.string(('pt > %f && abs(eta) < %f')%(jetPtCut,jetEtaCut)) + ) + + if muonCollection != '': + setattr(jetsNotOverlappingWithLeptons.checkOverlaps,"muons",cms.PSet( src = cms.InputTag(muonCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + if electronCollection != '': + setattr(jetsNotOverlappingWithLeptons.checkOverlaps,"electrons",cms.PSet( src = cms.InputTag(electronCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + + if tauCollection != '': + setattr(jetsNotOverlappingWithLeptons.checkOverlaps,"taus",cms.PSet( src = cms.InputTag(tauCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + + setattr(process,jetCollection+label,jetsNotOverlappingWithLeptons); + + +## clean gen jets from gen leptons +def cleanGenJetsFromGenLeptons (process, + jetCollection = "", + genLeptonCollection = "", + jetPtCut = 10, + jetEtaCut = 5, + dRCleaning = 0.3): + + + jetsNotOverlappingWithLeptons = cms.EDProducer("PATJetCleaner", + src = cms.InputTag(jetCollection), + preselection = cms.string(''), + checkOverlaps = cms.PSet(), + finalCut = cms.string(('pt > %f && abs(eta) < %f')%(jetPtCut,jetEtaCut)) + ) + + + setattr(jetsNotOverlappingWithLeptons.checkOverlaps,"muons",cms.PSet( src = cms.InputTag(genLeptonCollection), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(dRCleaning), + checkRecoComponents = cms.bool(False), + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True))) + + + setattr(process,jetCollection+"Cleaned",jetsNotOverlappingWithLeptons); diff --git a/RecoMET/METPUSubtraction/python/MVAMETConfiguration_cff.py b/RecoMET/METPUSubtraction/python/MVAMETConfiguration_cff.py new file mode 100644 index 0000000000000..ed75c953d65b3 --- /dev/null +++ b/RecoMET/METPUSubtraction/python/MVAMETConfiguration_cff.py @@ -0,0 +1,250 @@ +import FWCore.ParameterSet.Config as cms +import sys + + +import FWCore.ParameterSet.Config as cms +import sys + +from RecoMET.METPUSubtraction.LeptonSelectionTools_cff import applyMuonID +from RecoMET.METPUSubtraction.LeptonSelectionTools_cff import applyElectronID +from RecoMET.METPUSubtraction.LeptonSelectionTools_cff import applyTauID +from RecoMET.METPUSubtraction.LeptonSelectionTools_cff import cleanJetsFromLeptons +from RecoMET.METProducers.PFMET_cfi import pfMet +from JetMETCorrections.Type1MET.correctionTermsPfMetType1Type2_cff import corrPfMetType1 +from JetMETCorrections.Type1MET.correctedMet_cff import pfMetT1 +from PhysicsTools.SelectorUtils.tools.vid_id_tools import switchOnVIDElectronIdProducer, switchOnVIDPhotonIdProducer, DataFormat, setupAllVIDIdsInModule, setupVIDElectronSelection, setupVIDPhotonSelection + +def runMVAMET(process, + srcMuons = "slimmedMuons", ## inputMuonCollection + muonTypeID = "Tight", ## type of muon ID to be applied + typeIsoMuons = "dBeta", ## isolation type to be used for muons + relativeIsoCutMuons = 0.12, + srcElectrons = "slimmedElectrons", + electronTypeID= "Tight", + electronID_map = 'egmGsfElectronIDs:cutBasedElectronID-Spring15-25ns-V1-standalone-tight', + typeIsoElectrons = "rhoCorr", + relativeIsoCutEletrons = 0.12, + srcTaus = "slimmedTaus", + tauTypeID = "Loose", + jetCollectionPF = "slimmedJets", + dRCleaning= 0.3, + jetPtCut = 15, + jetEtaCut =5., + saveMapForTraining = False, + debug = False + ): + + if not hasattr(process,"egmGsfElectronIDs"): + electronIdModules = ['RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Spring15_25ns_V1_cff', + 'RecoEgamma.ElectronIdentification.Identification.heepElectronID_HEEPV51_cff'] + + switchOnVIDElectronIdProducer(process, DataFormat.MiniAOD) + + for idMod in electronIdModules: + setupAllVIDIdsInModule(process, idMod, setupVIDElectronSelection) + + if not hasattr(process,"VersionedPhotonIdProducer"): + switchOnVIDPhotonIdProducer(process, DataFormat.MiniAOD) + + photonIdModules = ['RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_PHYS14_PU20bx25_V2_cff'] + + for idMod in photonIdModules: + setupAllVIDIdsInModule(process, idMod, setupVIDPhotonSelection) + + ## create the Path + process.jmfw_analyzers = cms.Sequence() + if( not hasattr(process, "p")): + process.p = cms.Path() + process.p += process.jmfw_analyzers + # additional contribution from hadronically decaying taus + from RecoMET.METPUSubtraction.tausSignificance import tausSignificance, tauMET, tauPFMET, tauDecayProducts, allDecayProducts + process.tausSignificance = tausSignificance + process.tauMET = tauMET + process.tauMET.srcPFCands = cms.InputTag("packedPFCandidates") + process.tauPFMET = tauPFMET + process.tauDecayProducts = tauDecayProducts + process.allDecayProducts = allDecayProducts + + relativeIsoCutMuonsLoose = relativeIsoCutMuons+0.05; + relativeIsoCutEletronsLoose = relativeIsoCutEletrons+0.05; + + ### run Muon ID + applyMuonID(process, + src = srcMuons, + label = muonTypeID, + iso_map_charged_hadron = '', + iso_map_neutral_hadron = '', + iso_map_photon = '', + typeIso = typeIsoMuons, + relativeIsolationCutVal = relativeIsoCutMuons + ) + + + ### run Electron ID + applyElectronID(process, + label = electronTypeID, + src = srcElectrons, + iso_map_charged_hadron = '', + iso_map_neutral_hadron = '', + iso_map_photon = '', + typeIso = typeIsoElectrons, + electron_id_map = electronID_map, + relativeIsolationCutVal = relativeIsoCutEletrons + ) + + ### run tau ID + applyTauID( process, + label = tauTypeID, + src = srcTaus, + muonCollection = srcMuons+muonTypeID, + electronCollection = srcElectrons+electronTypeID) + + ## jet lepton cleaning + + cleanJetsFromLeptons(process, + label = "Cleaned", + jetCollection = jetCollectionPF, + muonCollection = srcMuons+muonTypeID, + electronCollection = srcElectrons+electronTypeID, + tauCollection = srcTaus+tauTypeID+"Cleaned", + jetPtCut = jetPtCut, + jetEtaCut = jetEtaCut, + dRCleaning = dRCleaning) + + + #### Input definitions like in classic MVA MET + #### tracks from PV + process.pfCHS = cms.EDFilter("CandPtrSelector", + cut = cms.string('fromPV'), + src = cms.InputTag("packedPFCandidates") + ) + process.pfChargedPV = cms.EDFilter("CandPtrSelector", + cut = cms.string('fromPV && charge !=0'), + src = cms.InputTag("packedPFCandidates") + ) + #### tracks not from PV + process.pfChargedPU = cms.EDFilter("CandPtrSelector", + cut = cms.string('!fromPV && charge !=0'), + src = cms.InputTag("packedPFCandidates") + ) + process.pfNeutrals = cms.EDFilter("CandPtrSelector", + cut = cms.string('charge ==0'), + src = cms.InputTag("packedPFCandidates") + ) + #### Neutrals in Jets passing PU Jet ID + #### and Neutrals in Jets not passing PU Jet ID + process.neutralInJets = cms.EDProducer("neutralCandidatePUIDJets", + srcJets = cms.InputTag(jetCollectionPF+"Cleaned"), + srcCandidates = cms.InputTag("pfNeutrals"), + neutralParticlesPVJetsLabel = cms.string("neutralPassingPUIDJets"), + neutralParticlesPUJetsLabel = cms.string("neutralFailingPUIDJets"), + neutralParticlesUnclusteredLabel = cms.string("neutralParticlesUnclustered"), + jetPUDIWP = cms.string("user"), + jetPUIDMapLabel = cms.string("fullDiscriminant")) + + + #### Merge collections to produce corresponding METs + process.pfMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag( + cms.InputTag("pfChargedPV"), + cms.InputTag("pfChargedPU"), + cms.InputTag("neutralInJets", "neutralPassingPUIDJets"), + cms.InputTag("neutralInJets", "neutralFailingPUIDJets"), + cms.InputTag("neutralInJets", "neutralParticlesUnclustered") + )) + #### Track MET + process.pfTrackMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag(cms.InputTag("pfChargedPV"))) + ## No-PU MET + process.pfNoPUMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag( cms.InputTag("pfChargedPV"), + cms.InputTag("neutralInJets", "neutralPassingPUIDJets"))) + ## PU corrected MET + process.pfPUCorrectedMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag( cms.InputTag("pfChargedPV"), + cms.InputTag("neutralInJets", "neutralPassingPUIDJets"), + cms.InputTag("neutralInJets", "neutralParticlesUnclustered"))) + ## PU MET + process.pfPUMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag( cms.InputTag("pfChargedPU"), + cms.InputTag("neutralInJets", "neutralFailingPUIDJets"))) + + ##Experimental + process.pfChargedPUMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag(cms.InputTag("pfChargedPU"))) + process.pfNeutralPUMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag(cms.InputTag("neutralInJets", "neutralFailingPUIDJets"))) + process.pfNeutralPVMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag(cms.InputTag("neutralInJets", "neutralPassingPUIDJets"))) + process.pfNeutralUnclusteredMETCands = cms.EDProducer("CandViewMerger", src = cms.VInputTag(cms.InputTag("neutralInJets", "neutralParticlesUnclustered"))) + + + from PhysicsTools.PatAlgos.producersLayer1.metProducer_cfi import patMETs + patMETsForMVA = patMETs.clone() + patMETsForMVA.computeMETSignificance = cms.bool(False) + patMETsForMVA.addGenMET = cms.bool(False) + patMETsForMVA.srcJets = cms.InputTag(jetCollectionPF) + + setattr(patMETsForMVA,"srcLeptons", cms.VInputTag(srcMuons+muonTypeID,srcElectrons+electronTypeID,srcTaus+tauTypeID+"Cleaned")) + + # get jets for T1 Correction + from RecoJets.JetProducers.ak4PFJets_cfi import ak4PFJets + from JetMETCorrections.Configuration.JetCorrectors_cff import ak4PFCHSL1FastL2L3Corrector,ak4PFCHSL3AbsoluteCorrector,ak4PFCHSL2RelativeCorrector,ak4PFCHSL1FastjetCorrector, ak4PFCHSL1FastL2L3ResidualCorrector, ak4PFCHSResidualCorrector + process.ak4PFCHSL1FastL2L3Corrector = ak4PFCHSL1FastL2L3Corrector + process.ak4PFCHSL1FastL2L3ResidualCorrector = ak4PFCHSL1FastL2L3ResidualCorrector + process.ak4PFCHSResidualCorrector = ak4PFCHSResidualCorrector + process.ak4PFCHSL3AbsoluteCorrector = ak4PFCHSL3AbsoluteCorrector + process.ak4PFCHSL2RelativeCorrector = ak4PFCHSL2RelativeCorrector + process.ak4PFCHSL1FastjetCorrector = ak4PFCHSL1FastjetCorrector + + + for met in ["pfMET", "pfTrackMET", "pfNoPUMET", "pfPUCorrectedMET", "pfPUMET", "pfChargedPUMET", "pfNeutralPUMET", "pfNeutralPVMET", "pfNeutralUnclusteredMET"]: + # create PF METs + setattr(process, met, pfMet.clone(src = cms.InputTag(met+"Cands"), alias = cms.string(met))) + # create Jets + setattr(process, "ak4JetsFor"+met, ak4PFJets.clone(src = cms.InputTag(met+"Cands"))) + setattr(process, "corr"+met, corrPfMetType1.clone(src = cms.InputTag("ak4JetsFor"+met))) + # derive corrections and apply them + setattr(process, met+"T1", pfMetT1.clone( src= cms.InputTag(met), srcCorrections=cms.VInputTag(cms.InputTag("corr"+met, "type1")))) + # convert METs to pat objects + setattr(process, "pat"+met, patMETsForMVA.clone(metSource = cms.InputTag(met))) + setattr(process, "pat"+met+"T1", patMETsForMVA.clone(metSource = cms.InputTag(met+"T1"))) + + process.pfChs = cms.EDProducer("CandViewMerger", src = cms.VInputTag( + cms.InputTag("pfChargedPV"), + cms.InputTag("neutralInJets", "neutralPassingPUIDJets"), + cms.InputTag("neutralInJets", "neutralFailingPUIDJets"), + cms.InputTag("neutralInJets", "neutralParticlesUnclustered") + )) + process.ak4JetsForpfMET.src = cms.InputTag("pfChs") + + ### MVA MET + process.MVAMET = cms.EDProducer("MVAMET", + debug = cms.bool(debug), + requireOS = cms.bool(True), + combineNLeptons = cms.int32(2), + MVAMETLabel = cms.string("MVAMET"), + srcMETs = cms.VInputTag( + cms.InputTag("slimmedMETs"), + cms.InputTag("patpfMET"), + cms.InputTag("patpfMETT1"), + cms.InputTag("patpfTrackMET"), + #cms.InputTag("patpfTrackMETT1"), + cms.InputTag("patpfNoPUMET"), + #cms.InputTag("patpfNoPUMETT1"), + cms.InputTag("patpfPUCorrectedMET"), + #cms.InputTag("patpfPUCorrectedMETT1"), + cms.InputTag("patpfPUMET"), + #cms.InputTag("patpfPUMETT1"), + cms.InputTag("slimmedMETsPuppi"), + #cms.InputTag("patpfChargedPUMET"), + #cms.InputTag("patpfNeutralPUMET"), + #cms.InputTag("patpfNeutralPVMET"), + #cms.InputTag("patpfNeutralUnclusteredMET") + ), + #inputMETFlags = cms.vint32(0,0,0,1,1,0,0,0,0,3,3,0,3,3,2,3), # use this flags if all MET collections are used above + inputMETFlags = cms.vint32(0,0,0,1,0,0,3,0), + srcJets = cms.InputTag(jetCollectionPF+"Cleaned"), + srcVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + srcTaus = cms.InputTag(srcTaus+tauTypeID+"Cleaned"), + srcMuons = cms.InputTag(srcMuons+muonTypeID), + weightFile = cms.FileInPath('RecoMET/METPUSubtraction/data/weightfile.root'), + #srcLeptons = cms.VInputTag("slimmedMuons", "slimmedElectrons", "slimmedTaus"), # to produce all possible combinations + srcLeptons = cms.VInputTag(srcMuons+muonTypeID,srcElectrons+electronTypeID,srcTaus+tauTypeID+"Cleaned"), # to produce a selection specifically designed for trainings + useTauSig = cms.bool(True), + tausSignificance = cms.InputTag('tausSignificance', 'METCovariance'), + saveMap = cms.bool(saveMapForTraining), + permuteLeptonsWithinPlugin = cms.bool(True) + ) diff --git a/RecoMET/METPUSubtraction/python/jet_recorrections.py b/RecoMET/METPUSubtraction/python/jet_recorrections.py new file mode 100644 index 0000000000000..b44079dfc85a2 --- /dev/null +++ b/RecoMET/METPUSubtraction/python/jet_recorrections.py @@ -0,0 +1,47 @@ +import FWCore.ParameterSet.Config as cms +from CondCore.CondDB.CondDB_cfi import * + +def loadLocalSqlite(process, sqliteFilename, tag = 'JetCorrectorParametersCollection_Fall15_25nsV2_MC_AK4PFchs'): + process.load("CondCore.DBCommon.CondDBCommon_cfi") + process.jec = cms.ESSource("PoolDBESSource", + DBParameters = cms.PSet( + messageLevel = cms.untracked.int32(0) + ), + timetype = cms.string('runnumber'), + toGet = cms.VPSet( + cms.PSet( + record = cms.string('JetCorrectionsRecord'), + tag = cms.string(tag), + label = cms.untracked.string('AK4PFchs') + ), + ), + connect = cms.string('sqlite:' + sqliteFilename) + ) + ## add an es_prefer statement to resolve a possible conflict from simultaneous connection to a global tag + process.es_prefer_jec = cms.ESPrefer('PoolDBESSource','jec') + +def recorrectJets(process, isData = False): + ## https://twiki.cern.ch/twiki/bin/view/CMSPublic/WorkBookJetEnergyCorrections#CorrPatJets + from PhysicsTools.PatAlgos.producersLayer1.jetUpdater_cff import updatedPatJetCorrFactors + process.patJetCorrFactorsReapplyJEC = updatedPatJetCorrFactors.clone( + src = cms.InputTag("slimmedJets"), + levels = ['L1FastJet', 'L2Relative', 'L3Absolute'], + payload = 'AK4PFchs' ) # Make sure to choose the appropriate levels and payload here! + from PhysicsTools.PatAlgos.producersLayer1.jetUpdater_cff import updatedPatJets + process.patJetsReapplyJEC = updatedPatJets.clone( + jetSource = cms.InputTag("slimmedJets"), + jetCorrFactorsSource = cms.VInputTag(cms.InputTag("patJetCorrFactorsReapplyJEC")) + ) + if(isData): + process.patJetCorrFactorsReapplyJEC.levels = ['L1FastJet', 'L2Relative', 'L3Absolute', 'L2L3Residual'] +# if( not hasattr(process, "p")): +# process.p = cms.Path() +# process.p += cms.Sequence( process.patJetCorrFactorsReapplyJEC + process. patJetsReapplyJEC ) + +def reapplyPUJetID(process, srcJets = cms.InputTag("slimmedJets")): + from RecoJets.JetProducers.PileupJetID_cfi import pileupJetId + process.pileupJetIdUpdated = pileupJetId.clone( + jets = srcJets, + inputIsCorrected = True, + applyJec = True, + vertexes = cms.InputTag("offlineSlimmedPrimaryVertices") ) diff --git a/RecoMET/METPUSubtraction/python/mitigatedMETSequence_cff.py b/RecoMET/METPUSubtraction/python/mitigatedMETSequence_cff.py deleted file mode 100644 index 50fd22753ee54..0000000000000 --- a/RecoMET/METPUSubtraction/python/mitigatedMETSequence_cff.py +++ /dev/null @@ -1,190 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -from RecoMET.METPUSubtraction.objectSelection_cff import * - - -##================================================ -## MVA MET sequence -##================================================ -from JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff import * -from JetMETCorrections.Configuration.DefaultJEC_cff import * -from RecoJets.JetProducers.PileupJetIDParams_cfi import JetIdParams - -from RecoJets.JetProducers.kt4PFJets_cfi import kt4PFJets -kt6PFJets = kt4PFJets.clone(rParam = cms.double(0.6), doRhoFastjet = True ) - -calibratedAK4PFJetsForPFMVAMEt = cms.EDProducer('PFJetCorrectionProducer', - src = cms.InputTag('ak4PFJets'), - correctors = cms.vstring("ak4PFL1FastL2L3") # NOTE: use "ak5PFL1FastL2L3" for MC / "ak5PFL1FastL2L3Residual" for Data -) - -pfMVAMEt = cms.EDProducer("PFMETProducerMVA", - srcCorrJets = cms.InputTag('calibratedAK4PFJetsForPFMVAMEt'), - srcUncorrJets = cms.InputTag('ak4PFJets'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcVertices = cms.InputTag('offlinePrimaryVertices'), - srcLeptons = cms.VInputTag('selectedElectrons', - 'selectedMuons', - 'selectedTaus', - 'selectedPhotons', - 'selectedJets'), - minNumLeptons = cms.int32(0), - srcRho = cms.InputTag('fixedGridRhoFastjetAll'), - globalThreshold = cms.double(-1.), - minCorrJetPt = cms.double(-1.), - inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_June2013_type1.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_June2013_type1.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root') - ), - loadMVAfromDB = cms.bool(False), - is42 = cms.bool(False), # CV: set this flag to true if you are running mvaPFMET in CMSSW_4_2_x - corrector = cms.string("ak4PFL1Fastjet"), - useType1 = cms.bool(True), - useOld42 = cms.bool(False), - dZcut = cms.double(0.1), - impactParTkThreshold = cms.double(0.), - tmvaWeights = cms.string("RecoJets/JetProducers/data/TMVAClassificationCategory_JetID_MET_53X_Dec2012.weights.xml.gz"), - tmvaMethod = cms.string("JetID"), - version = cms.int32(-1), - cutBased = cms.bool(False), - tmvaVariables = cms.vstring( - "nvtx", - "jetPt", - "jetEta", - "jetPhi", - "dZ", - "beta", - "betaStar", - "nCharged", - "nNeutrals", - "dR2Mean", - "ptD", - "frac01", - "frac02", - "frac03", - "frac04", - "frac05" - ), - tmvaSpectators = cms.vstring(), - JetIdParams = JetIdParams, - verbosity = cms.int32(0) -) - -pfMVAMEtSequence = cms.Sequence( - kt6PFJets* - calibratedAK4PFJetsForPFMVAMEt* - pfMVAMEt -) - - -##================================================ -## Pf No Pileup MET sequence -##================================================ -pfNoPUMEtSequence = cms.Sequence() - -from JetMETCorrections.Configuration.JetCorrectionServices_cff import * -calibratedAK4PFJetsForPFNoPUMEt = cms.EDProducer('PFJetCorrectionProducer', - src = cms.InputTag('ak4PFJets'), - correctors = cms.vstring('ak4PFL1FastL2L3Residual') # NOTE: use "ak4PFL1FastL2L3" for MC / "ak4PFL1FastL2L3Residual" for Data -) -ak4PFJetSequenceForPFNoPUMEt = cms.Sequence(calibratedAK4PFJetsForPFNoPUMEt) -pfNoPUMEtSequence += ak4PFJetSequenceForPFNoPUMEt - -from RecoJets.JetProducers.PileupJetID_cfi import * -puJetIdForPFNoPUMEt = pileupJetId.clone( - algos = cms.VPSet( - full_53x, - cutbased, - PhilV1 - ), -# label = cms.string("fullId"), #MM does not work for weird reasons, cannot be cloned properly - produceJetIds = cms.bool(True), - runMvas = cms.bool(True), - jets = cms.InputTag("calibratedAK4PFJetsForPFNoPUMEt"), - applyJec = cms.bool(False), - inputIsCorrected = cms.bool(True), - ) -pfNoPUMEtSequence += puJetIdForPFNoPUMEt - -from JetMETCorrections.Type1MET.pfMETCorrectionType0_cfi import * -pfNoPUMEtSequence += type0PFMEtCorrection -pfCandidateToVertexAssociationForPFNoPUMEt = pfCandidateToVertexAssociation.clone( - MaxNumberOfAssociations = cms.int32(1), - doReassociation = cms.bool(False), - FinalAssociation = cms.untracked.int32(1), - nTrackWeight = cms.double(0.) -) -pfNoPUMEtSequence += pfCandidateToVertexAssociationForPFNoPUMEt -pfMETcorrType0ForPFNoPUMEt = pfMETcorrType0.clone( - srcPFCandidateToVertexAssociations = cms.InputTag('pfCandidateToVertexAssociationForPFNoPUMEt') -) -pfNoPUMEtSequence += pfMETcorrType0ForPFNoPUMEt - -jvfJetIdForPFNoPUMEt = cms.EDProducer("JVFJetIdProducer", - srcJets = cms.InputTag('calibratedAK4PFJetsForPFNoPUMEt'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcPFCandToVertexAssociations = cms.InputTag('pfCandidateToVertexAssociationForPFNoPUMEt'), - srcHardScatterVertex = cms.InputTag('selectedPrimaryVertexHighestPtTrackSumForPFMEtCorrType0'), - minTrackPt = cms.double(1.), - dZcut = cms.double(0.2), # cm - JVFcut = cms.double(0.75), - neutralJetOption = cms.string("noPU") -) -pfNoPUMEtSequence += jvfJetIdForPFNoPUMEt - -import RecoMET.METProducers.METSigParams_cfi as met_config -pfNoPUMEtData = cms.EDProducer("NoPileUpPFMEtDataProducer", - srcJets = cms.InputTag('calibratedAK4PFJetsForPFNoPUMEt'), - srcJetIds = cms.InputTag('puJetIdForPFNoPUMEt', 'full53xId'), - - minJetPt = cms.double(30.0), - jetIdSelection = cms.string('loose'), - jetEnOffsetCorrLabel = cms.string("ak4PFL1Fastjet"), - srcPFCandidates = cms.InputTag('particleFlow'), - srcPFCandToVertexAssociations = cms.InputTag('pfCandidateToVertexAssociationForPFNoPUMEt'), - srcJetsForMEtCov = cms.InputTag('ak4PFJets'), - minJetPtForMEtCov = cms.double(10.), - srcHardScatterVertex = cms.InputTag('selectedPrimaryVertexHighestPtTrackSumForPFMEtCorrType0'), - dZcut = cms.double(0.2), # cm - resolution = met_config.METSignificance_params, - verbosity = cms.int32(0) -) -pfNoPUMEtSequence += pfNoPUMEtData - -pfNoPUMEt = cms.EDProducer("NoPileUpPFMEtProducer", - srcMEt = cms.InputTag('pfMet'), - srcMEtCov = cms.InputTag(''), # NOTE: leave empty to take MET covariance matrix from reco::PFMET object //MM 08/29/14, bypass hardcoded as this variable has never been used so far - srcMVAMEtData = cms.InputTag('pfNoPUMEtData'), - srcLeptons = cms.VInputTag( 'selectedElectrons', - 'selectedMuons', - 'selectedTaus', - 'selectedPhotons', - 'selectedJets'), -# NOTE: you need to set this to collections of electrons, muons and tau-jets -#passing the lepton reconstruction & identification criteria applied in your analysis - srcMVAMEtDataLeptonMatch = cms.InputTag('pfNoPUMEtData'), - srcType0Correction = cms.InputTag('pfMETcorrType0ForPFNoPUMEt'), - sfNoPUjets = cms.double(1.0), - sfNoPUjetOffsetEnCorr = cms.double(0.0), - sfPUjets = cms.double(1.0), - sfNoPUunclChargedCands = cms.double(1.0), - sfPUunclChargedCands = cms.double(1.0), - sfUnclNeutralCands = cms.double(0.6), - sfType0Correction = cms.double(1.0), - sfLeptonIsoCones = cms.double(0.6), - resolution = met_config.METSignificance_params, - sfMEtCovMin = cms.double(0.6), - sfMEtCovMax = cms.double(1.0), - saveInputs = cms.bool(True), - verbosity = cms.int32(0) -) -pfNoPUMEtSequence += pfNoPUMEt - - -mitigatedMETSequence = cms.Sequence( -selectionSequenceForMVANoPUMET+ -pfMVAMEtSequence+ -pfNoPUMEtSequence -) diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_Data_cff.py b/RecoMET/METPUSubtraction/python/mvaPFMET_Data_cff.py deleted file mode 100644 index a805a479bc18e..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_Data_cff.py +++ /dev/null @@ -1,68 +0,0 @@ - -import FWCore.ParameterSet.Config as cms - -#from RecoMET.METProducers.PFMET_cfi import pfMet -from JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff import * -from JetMETCorrections.Configuration.DefaultJEC_cff import * -from RecoMET.METPUSubtraction.mvaPFMET_leptons_cfi import * -from RecoJets.JetProducers.PileupJetIDParams_cfi import JetIdParams - -calibratedAK5PFJetsForPFMEtMVA = cms.EDProducer('PFJetCorrectionProducer', - src = cms.InputTag('ak5PFJets'), - correctors = cms.vstring("ak5PFL1FastL2L3Residual") # NOTE: use "ak5PFL1FastL2L3" for MC / "ak5PFL1FastL2L3Residual" for Data -) - -pfMEtMVA = cms.EDProducer("PFMETProducerMVA", - srcCorrJets = cms.InputTag('calibratedAK5PFJetsForPFMEtMVA'), - srcUncorrJets = cms.InputTag('ak5PFJets'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcVertices = cms.InputTag('offlinePrimaryVertices'), - srcLeptons = cms.VInputTag(),#"isomuons","isoelectrons","isotaus") # NOTE: you need to set this to collections of electrons, muons and tau-jets - # passing the lepton reconstruction & identification criteria applied in your analysis - srcRho = cms.InputTag('kt6PFJets','rho'), - globalThreshold = cms.double(-1.),#pfMet.globalThreshold, - minCorrJetPt = cms.double(-1.), - inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_Dec2012.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_Dec2012.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root') - ), - corrector = cms.string("ak5PFL1Fastjet"), - useType1 = cms.bool(True), - useOld42 = cms.bool(False), - dZcut = cms.double(0.1), - impactParTkThreshold = cms.double(0.), - tmvaWeights = cms.string("RecoJets/JetProducers/data/TMVAClassificationCategory_JetID_MET_53X_Dec2012.weights.xml.gz"), - tmvaMethod = cms.string("JetID"), - version = cms.int32(-1), - cutBased = cms.bool(False), - tmvaVariables = cms.vstring( - "nvtx", - "jetPt", - "jetEta", - "jetPhi", - "dZ", - "beta", - "betaStar", - "nCharged", - "nNeutrals", - "dR2Mean", - "ptD", - "frac01", - "frac02", - "frac03", - "frac04", - "frac05", - ), - tmvaSpectators = cms.vstring(), - JetIdParams = JetIdParams, - label = cms.string("53XMet"), - verbosity = cms.int32(0) -) - -pfMEtMVAsequence = cms.Sequence( - #(isomuonseq+isotauseq+isoelectronseq)* - calibratedAK4PFJetsForPFMEtMVA* - pfMEtMVA -) diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py b/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py deleted file mode 100644 index 807ab018eadb2..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_cff.py +++ /dev/null @@ -1,100 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -#from RecoMET.METProducers.PFMET_cfi import pfMet -from JetMETCorrections.Configuration.JetCorrectors_cff import * -##from RecoMET.METPUSubtraction.mvaPFMET_leptons_cfi import * -## CV: importing mvaPFMET_leptons_cfi breaks produceAndDiscriminateHPSPFTaus sequence -## (hpsPFTauDiscriminationByDecayModeFinding module gets overwritten by None, -## in case RecoTauTag/Configuration/python/RecoPFTauTag_cff.py is loaded by -## by top-level cfg.py file before RecoMET/METPUSubtraction/python/mvaPFMET_cff.py gets loaded) -##from RecoJets.JetProducers.PileupJetIDCutParams_cfi import full_53x_wp - -calibratedAK4PFJetsForPFMVAMEt = cms.EDProducer('CorrectedPFJetProducer', - src = cms.InputTag('ak4PFJets'), - correctors = cms.VInputTag("ak4PFL1FastL2L3Corrector") # NOTE: use "ak4PFL1FastL2L3Corrector" for MC / "ak4PFL1FastL2L3ResidualCorrector" for Data -) -from JetMETCorrections.Configuration.JetCorrectionServices_cff import ak4PFL1Fastjet -from RecoJets.JetProducers.PileupJetID_cfi import pileupJetIdEvaluator -from RecoJets.JetProducers.PileupJetIDParams_cfi import JetIdParams -puJetIdForPFMVAMEt = pileupJetIdEvaluator.clone( - algos = cms.VPSet( - cms.PSet( - tmvaVariables = cms.vstring( - "nvtx", - "jetPt", - "jetEta", - "jetPhi", - "dZ", - "beta", - "betaStar", - "nCharged", - "nNeutrals", - "dR2Mean", - "ptD", - "frac01", - "frac02", - "frac03", - "frac04", - "frac05" - ), - etaBinnedWeights = cms.bool(False), - tmvaWeights = cms.string("RecoJets/JetProducers/data/TMVAClassificationCategory_JetID_MET_53X_Dec2012.weights.xml.gz"), - tmvaMethod = cms.string("JetID"), - tmvaSpectators = cms.vstring(), - JetIdParams = JetIdParams, - impactParTkThreshold = cms.double(0.), - version = cms.int32(-1), - cutBased = cms.bool(False), - label = cms.string("full") - ) - ), - produceJetIds = cms.bool(True), - runMvas = cms.bool(True), - jets = cms.InputTag("calibratedAK4PFJetsForPFMVAMEt"),#calibratedAK4PFJetsForPFMVAMEt - applyJec = cms.bool(True), - inputIsCorrected = cms.bool(True), - jec = cms.string("AK4PF"), -) - - - - -pfMVAMEt = cms.EDProducer("PFMETProducerMVA", - srcCorrJets = cms.InputTag('calibratedAK4PFJetsForPFMVAMEt'), - srcUncorrJets = cms.InputTag('ak4PFJets'), - srcMVAPileupJetId = cms.InputTag('puJetIdForPFMVAMEt','fullDiscriminant'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcVertices = cms.InputTag('offlinePrimaryVertices'), - srcLeptons = cms.VInputTag(),#"isomuons","isoelectrons","isotaus") # NOTE: you need to set this to collections of electrons, muons and tau-jets - # passing the lepton reconstruction & identification criteria applied in your analysis - minNumLeptons = cms.int32(0), - globalThreshold = cms.double(-1.),#pfMet.globalThreshold, - minCorrJetPt = cms.double(-1.), - inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru_7_4_X_miniAOD_25NS_July2015.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrphi_7_4_X_miniAOD_25NS_July2015.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_7_4_X_miniAOD_25NS_July2015.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_7_4_X_miniAOD_25NS_July2015.root') ), - inputRecords = cms.PSet( - U = cms.string("RecoilCor"), - DPhi = cms.string("PhiCor"), - CovU1 = cms.string("CovU1"), - CovU2 = cms.string("CovU2") - ), - loadMVAfromDB = cms.bool(False), - - corrector = cms.InputTag("ak4PFL1FastjetCorrector"), - useType1 = cms.bool(True), - dZcut = cms.double(0.1), - - verbosity = cms.int32(0) -) - - - -pfMVAMEtSequence = cms.Sequence( - #(isomuonseq+isotauseq+isoelectronseq)* - calibratedAK4PFJetsForPFMVAMEt* - puJetIdForPFMVAMEt* - pfMVAMEt -) diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_db_cfi.py b/RecoMET/METPUSubtraction/python/mvaPFMET_db_cfi.py deleted file mode 100644 index dff77fa4ad420..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_db_cfi.py +++ /dev/null @@ -1,29 +0,0 @@ -from CondCore.DBCommon.CondDBSetup_cfi import * - -## mvaPFMEtGBRForestsFromDB = cms.ESSource("PoolDBESSource", -## CondDBSetup, -## DumpStat = cms.untracked.bool(True), -## toGet = cms.VPSet( -## cms.PSet( -## record = cms.string('GBRWrapperRcd'), -## tag = cms.string('mvaPFMET_53_Dec2012_U'), -## label = cms.untracked.string('mvaPFMET_53_Dec2012_U') -## ), -## cms.PSet( -## record = cms.string('GBRWrapperRcd'), -## tag = cms.string('mvaPFMET_53_Dec2012_DPhi'), -## label = cms.untracked.string('mvaPFMET_53_Dec2012_DPhi') -## ), -## cms.PSet( -## record = cms.string('GBRWrapperRcd'), -## tag = cms.string('mvaPFMET_53_Dec2012_CovU1'), -## label = cms.untracked.string('mvaPFMET_53_Dec2012_CovU1') -## ), -## cms.PSet( -## record = cms.string('GBRWrapperRcd'), -## tag = cms.string('mvaPFMET_53_Dec2012_CovU2'), -## label = cms.untracked.string('mvaPFMET_53_Dec2012_CovU2') -## ) -## ), -## connect = cms.string('sqlite_fip:RecoMET/METPUSubtraction/data/mvaPFMEt_53_Dec2012.db') -## ) diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cff.py b/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cff.py deleted file mode 100644 index 3a2d5d038bf35..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cff.py +++ /dev/null @@ -1,71 +0,0 @@ - -import FWCore.ParameterSet.Config as cms - -#from RecoMET.METProducers.PFMET_cfi import pfMet -from JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff import * -from JetMETCorrections.Configuration.DefaultJEC_cff import * -from RecoMET.METPUSubtraction.mvaPFMET_leptons_cfi import * -from RecoJets.JetProducers.PileupJetIDParams_cfi import JetIdParams - -calibratedAK4PFJetsForPFMEtMVA = cms.EDProducer('PFJetCorrectionProducer', - src = cms.InputTag('ak4PFJets'), - correctors = cms.vstring("ak4PFL1FastL2L3") # NOTE: use "ak4PFL1FastL2L3" for MC / "ak4PFL1FastL2L3Residual" for Data -) - -pfMEtMVA = cms.EDProducer("PFMETProducerMVA", - srcCorrJets = cms.InputTag('calibratedAK4PFJetsForPFMEtMVA'), - srcUncorrJets = cms.InputTag('ak4PFJets'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcVertices = cms.InputTag('offlinePrimaryVertices'), - srcLeptons = cms.VInputTag(),#"isomuons","isoelectrons","isotaus") # NOTE: you need to set this to collections of electrons, muons and tau-jets - # passing the lepton reconstruction & identification criteria applied in your analysis - minNumLeptons = cms.int32(0), - srcRho = cms.InputTag('kt6PFJets','rho'), - globalThreshold = cms.double(-1.),#pfMet.globalThreshold, - minCorrJetPt = cms.double(-1.), - inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_June2013_type1.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_June2013_type1.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root') - ), - loadMVAfromDB = cms.bool(False), - is42 = cms.bool(False), # CV: set this flag to true if you are running mvaPFMET in CMSSW_4_2_x - corrector = cms.string("ak4PFL1Fastjet"), - useType1 = cms.bool(True), - useOld42 = cms.bool(False), - dZcut = cms.double(0.1), - impactParTkThreshold = cms.double(0.), - tmvaWeights = cms.string("RecoJets/JetProducers/data/TMVAClassificationCategory_JetID_MET_53X_Dec2012.weights.xml"), - tmvaMethod = cms.string("JetID"), - version = cms.int32(-1), - cutBased = cms.bool(False), - tmvaVariables = cms.vstring( - "nvtx", - "jetPt", - "jetEta", - "jetPhi", - "dZ", - "beta", - "betaStar", - "nCharged", - "nNeutrals", - "dR2Mean", - "ptD", - "frac01", - "frac02", - "frac03", - "frac04", - "frac05", - ), - tmvaSpectators = cms.vstring(), - JetIdParams = JetIdParams, - verbosity = cms.int32(0) -) - -pfMEtMVAsequence = cms.Sequence( - (isomuonseq+isotauseq+isoelectronseq)* - calibratedAK4PFJetsForPFMEtMVA* - pfMEtMVA - ) - diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cfi.py b/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cfi.py deleted file mode 100644 index 362093dec74cc..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_cfi.py +++ /dev/null @@ -1,119 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -# Single muon for Wjets -isomuons = cms.EDFilter( - "MuonSelector", - src = cms.InputTag('muons'), - cut = cms.string( "(isTrackerMuon) && std::abs(eta) < 2.5 && pt > 9.5"+#17. "+ - "&& isPFMuon"+ - "&& globalTrack.isNonnull"+ - "&& innerTrack.hitPattern.numberOfValidPixelHits > 0"+ - "&& innerTrack.normalizedChi2 < 10"+ - "&& numberOfMatches > 0"+ - "&& innerTrack.hitPattern.numberOfValidTrackerHits>5"+ - "&& globalTrack.hitPattern.numberOfValidHits>0"+ - "&& (pfIsolationR03.sumChargedHadronPt+pfIsolationR03.sumNeutralHadronEt+pfIsolationR03.sumPhotonEt)/pt < 0.3"+ - "&& std::abs(innerTrack().dxy)<2.0" - ), - filter = cms.bool(False) - ) - - -isoelectrons = cms.EDFilter( - "GsfElectronSelector", - src = cms.InputTag('gsfElectrons'), - cut = cms.string( - "std::abs(eta) < 2.5 && pt > 9.5" + - "&& gsfTrack.trackerExpectedHitsInner.numberOfHits == 0" + -# "&& (pfIsolationVariables.chargedHadronIso+pfIsolationVariables.neutralHadronIso)/et < 0.3" + - "&& (isolationVariables03.tkSumPt)/et < 0.2" + - "&& ((std::abs(eta) < 1.4442 " + - "&& std::abs(deltaEtaSuperClusterTrackAtVtx) < 0.007"+ - "&& std::abs(deltaPhiSuperClusterTrackAtVtx) < 0.8" + - "&& sigmaIetaIeta < 0.01" + - "&& hcalOverEcal < 0.15" + - "&& std::abs(1./superCluster.energy - 1./p) < 0.05)"+ - "|| (std::abs(eta) > 1.566 "+ - "&& std::abs(deltaEtaSuperClusterTrackAtVtx) < 0.009"+ - "&& std::abs(deltaPhiSuperClusterTrackAtVtx) < 0.10" + - "&& sigmaIetaIeta < 0.03" + - "&& hcalOverEcal < 0.10" + - "&& std::abs(1./superCluster.energy - 1./p) < 0.05))" - ), - filter = cms.bool(False) - ) - -from RecoJets.Configuration.RecoPFJets_cff import kt6PFJets as dummy -kt6PFJetsForRhoComputationVoronoiMet = dummy.clone( - doRhoFastjet = True, - voronoiRfact = 0.9 - ) - -from RecoTauTag.RecoTau.PFRecoTauDiscriminationByHPSSelection_cfi import hpsSelectionDiscriminator -hpsPFTauDiscriminationByDecayModeFinding = hpsSelectionDiscriminator.clone( - PFTauProducer = cms.InputTag('hpsPFTauProducer') - ) - -from RecoTauTag.RecoTau.TauDiscriminatorTools import requireLeadTrack -# Define decay mode prediscriminant -requireDecayMode = cms.PSet( - BooleanOperator = cms.string("and"), - decayMode = cms.PSet( - Producer = cms.InputTag('hpsPFTauDiscriminationByDecayModeFinding'), - cut = cms.double(0.5) - ) - ) - -from RecoTauTag.Configuration.HPSPFTaus_cff import hpsPFTauDiscriminationByLooseCombinedIsolationDBSumPtCorr3Hits - -hpsPFTauDiscriminationAgainstMuon2 = cms.EDProducer("PFRecoTauDiscriminationAgainstMuon2", - PFTauProducer = cms.InputTag('hpsPFTauProducer'), - Prediscriminants = requireDecayMode.clone(), - discriminatorOption = cms.string('loose'), # available options are: 'loose', 'medium', 'tight' - HoPMin = cms.double(0.2) - ) - - -hpsPFTauDiscriminationByMVAIsolation = cms.EDProducer( - "PFRecoTauDiscriminationByMVAIsolation", - PFTauProducer = cms.InputTag('hpsPFTauProducer'), - rhoProducer = cms.InputTag('kt6PFJetsForRhoComputationVoronoiMet','rho'), - Prediscriminants = requireDecayMode.clone(), - gbrfFilePath = cms.FileInPath('RecoTauTag/RecoTau/data/gbrfTauIso_v2.root'), - returnMVA = cms.bool(False), - mvaMin = cms.double(0.8), - ) - -isotaus = cms.EDFilter( - "PFTauSelector", - src = cms.InputTag('hpsPFTauProducer'), - BooleanOperator = cms.string("and"), - discriminators = cms.VPSet( - cms.PSet( discriminator=cms.InputTag("hpsPFTauDiscriminationByDecayModeFinding"), selectionCut=cms.double(0.5)), - #cms.PSet( discriminator=cms.InputTag("hpsPFTauDiscriminationByMVAIsolation"), selectionCut=cms.double(0.5)), - cms.PSet( discriminator=cms.InputTag("hpsPFTauDiscriminationByLooseCombinedIsolationDBSumPtCorr3Hits"), selectionCut=cms.double(0.5)), - cms.PSet( discriminator=cms.InputTag("hpsPFTauDiscriminationByLooseElectronRejection"), selectionCut=cms.double(0.5)), - cms.PSet( discriminator=cms.InputTag("hpsPFTauDiscriminationAgainstMuon2"), selectionCut=cms.double(0.5)) - ), - cut = cms.string("std::abs(eta) < 2.3 && pt > 19.0 "), - filter = cms.bool(False) - ) - -isomuonseq = cms.Sequence(isomuons) -isoelectronseq = cms.Sequence(isoelectrons) -isotauseq = cms.Sequence( - hpsPFTauDiscriminationByLooseCombinedIsolationDBSumPtCorr3Hits* - #kt6PFJetsForRhoComputationVoronoiMet* - #hpsPFTauDiscriminationByMVAIsolation* - hpsPFTauDiscriminationAgainstMuon2* - isotaus - ) - -leptonSelection = cms.PSet( - SelectEvents = cms.PSet( - SelectEvents = cms.vstring( - 'isomuonseq', - 'isoelectronseq', - 'isotauseq') - ) - ) diff --git a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_data_cff.py b/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_data_cff.py deleted file mode 100644 index a8b995631f394..0000000000000 --- a/RecoMET/METPUSubtraction/python/mvaPFMET_leptons_data_cff.py +++ /dev/null @@ -1,68 +0,0 @@ - -import FWCore.ParameterSet.Config as cms - -#from RecoMET.METProducers.PFMET_cfi import pfMet -from JetMETCorrections.Configuration.JetCorrectionServicesAllAlgos_cff import * -from JetMETCorrections.Configuration.DefaultJEC_cff import * -from RecoMET.METPUSubtraction.mvaPFMET_leptons_cfi import * -from RecoJets.JetProducers.PileupJetIDParams_cfi import JetIdParams - -calibratedAK4PFJetsForPFMEtMVA = cms.EDProducer('PFJetCorrectionProducer', - src = cms.InputTag('ak4PFJets'), - correctors = cms.vstring("ak4PFL1FastL2L3Residual") -) - -pfMEtMVA = cms.EDProducer("PFMETProducerMVA", - srcCorrJets = cms.InputTag('calibratedAK4PFJetsForPFMEtMVA'), - srcUncorrJets = cms.InputTag('ak4PFJets'), - srcPFCandidates = cms.InputTag('particleFlow'), - srcVertices = cms.InputTag('offlinePrimaryVertices'), - srcLeptons = cms.VInputTag("isomuons","isoelectrons","isotaus"),#"muons","hpsPFTauProducer"), # NOTE: you need to set this to collections of electrons, muons and tau-jets - # passing the lepton reconstruction & identification criteria applied in your analysis - srcRho = cms.InputTag('kt6PFJets','rho'), - globalThreshold = cms.double(-1.),#pfMet.globalThreshold, - minCorrJetPt = cms.double(-1.), - inputFileNames = cms.PSet( - U = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_June2013_type1.root'), - DPhi = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_June2013_type1.root'), - CovU1 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - CovU2 = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root') - ), - corrector = cms.string("ak4PFL1Fastjet"), - useType1 = cms.bool(True), - useOld42 = cms.bool(False), - dZcut = cms.double(0.1), - impactParTkThreshold = cms.double(0.), - tmvaWeights = cms.string("RecoJets/JetProducers/data/TMVAClassificationCategory_JetID_MET_53X_Dec2012.weights.xml"), - tmvaMethod = cms.string("JetID"), - version = cms.int32(-1), - cutBased = cms.bool(False), - tmvaVariables = cms.vstring( - "nvtx", - "jetPt", - "jetEta", - "jetPhi", - "dZ", - "beta", - "betaStar", - "nCharged", - "nNeutrals", - "dR2Mean", - "ptD", - "frac01", - "frac02", - "frac03", - "frac04", - "frac05", - ), - tmvaSpectators = cms.vstring(), - JetIdParams = JetIdParams, - verbosity = cms.int32(0) -) - -pfMEtMVAsequence = cms.Sequence( - (isomuonseq+isotauseq+isoelectronseq)* - calibratedAK4PFJetsForPFMEtMVA* - pfMEtMVA - ) - diff --git a/RecoMET/METPUSubtraction/python/objectSelection_cff.py b/RecoMET/METPUSubtraction/python/objectSelection_cff.py deleted file mode 100644 index da3760ece770d..0000000000000 --- a/RecoMET/METPUSubtraction/python/objectSelection_cff.py +++ /dev/null @@ -1,116 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -##====================================== -## Muons -##====================================== - -selectedMuons = cms.EDFilter( - "MuonSelector", - src = cms.InputTag('muons'), - cut = cms.string( "(isTrackerMuon) && std::abs(eta) < 2.5 && pt > 19.5"+#17. "+ - "&& isPFMuon"+ - "&& globalTrack.isNonnull"+ - "&& innerTrack.hitPattern.numberOfValidPixelHits > 0"+ - "&& innerTrack.normalizedChi2 < 10"+ - "&& numberOfMatches > 0"+ - "&& innerTrack.hitPattern.numberOfValidTrackerHits>5"+ - "&& globalTrack.hitPattern.numberOfValidHits>0"+ - "&& (pfIsolationR03.sumChargedHadronPt+std::max(0.,pfIsolationR03.sumNeutralHadronEt+pfIsolationR03.sumPhotonEt - 0.5*pfIsolationR03.sumPUPt))/pt < 0.3"+ - "&& std::abs(innerTrack().dxy)<2.0" - ), - filter = cms.bool(False) - ) - - -##====================================== -## Electrons -##====================================== - - -selectedElectrons = cms.EDFilter( - "GsfElectronSelector", - src = cms.InputTag('gedGsfElectrons'), - cut = cms.string( - "std::abs(eta) < 2.5 && pt > 19.5" + - "&& (gsfTrack.hitPattern().numberOfHits(\'MISSING_INNER_HITS\')<=1 )" + - "&& (pfIsolationVariables.sumChargedHadronPt+std::max(0.,pfIsolationVariables.sumNeutralHadronEt+pfIsolationVariables.sumPhotonEt - 0.5*pfIsolationVariables.sumPUPt))/et < 0.3" + - "&& ((std::abs(eta) < 1.4442 " + - "&& std::abs(deltaEtaSuperClusterTrackAtVtx) < 0.007"+ - "&& std::abs(deltaPhiSuperClusterTrackAtVtx) < 0.8" + - "&& sigmaIetaIeta < 0.01" + - "&& hcalOverEcal < 0.15" + - "&& std::abs(1./superCluster.energy - 1./p) < 0.05)"+ - "|| (std::abs(eta) > 1.566 "+ - "&& std::abs(deltaEtaSuperClusterTrackAtVtx) < 0.009"+ - "&& std::abs(deltaPhiSuperClusterTrackAtVtx) < 0.10" + - "&& sigmaIetaIeta < 0.03" + - "&& hcalOverEcal < 0.10" + - "&& std::abs(1./superCluster.energy - 1./p) < 0.05))" - ), - filter = cms.bool(False) - ) - - -##====================================== -## Taus -##====================================== - -selectedTaus = cms.EDFilter("PFTauSelector", - src = cms.InputTag('hpsPFTauProducer'), - discriminators = cms.VPSet( - cms.PSet( - discriminator = cms.InputTag('hpsPFTauDiscriminationByDecayModeFinding'), - selectionCut = cms.double(0.5) - ), - cms.PSet( - discriminator = cms.InputTag('hpsPFTauDiscriminationByLooseCombinedIsolationDBSumPtCorr3Hits'), - selectionCut = cms.double(0.5) - ) - ), - cut = cms.string("pt > 20. & std::abs(eta) < 2.3") -) - - - -##====================================== -## Photons -##====================================== - -selectedPhotons = cms.EDFilter("PhotonSelector", - src = cms.InputTag("photons"), - cut = cms.string( - "std::abs(eta) < 2.5 && pt > 19.5" + - "&& sigmaIetaIeta < 0.03" + - "&& hadronicOverEm < 0.05" + - "&& hasPixelSeed == 0" + - "&& (chargedHadronIso + neutralHadronIso + photonIso)/pt < 0.2" - ) -) - -##====================================== -## Jets -##====================================== - -#pileup jetId applied per default in the process - -jet_acc = '(pt >= 30 && std::abs(eta)<2.5)' -#jet_id = '' - -selectedJets = cms.EDFilter("PFJetSelector", - src = cms.InputTag("ak4PFJets"), - cut = cms.string( - "(pt >= 30 && std::abs(eta)<2.5)" + - "&& neutralHadronEnergyFraction < 0.99" + - "&& neutralEmEnergyFraction < 0.99" + - "&& getPFConstituents.size > 1" - ) -) - - -selectionSequenceForMVANoPUMET = cms.Sequence( -selectedMuons+ -selectedElectrons+ -selectedTaus+ -selectedPhotons+ -selectedJets -) diff --git a/RecoMET/METPUSubtraction/python/tausSignificance.py b/RecoMET/METPUSubtraction/python/tausSignificance.py new file mode 100644 index 0000000000000..083c680356bef --- /dev/null +++ b/RecoMET/METPUSubtraction/python/tausSignificance.py @@ -0,0 +1,26 @@ +import FWCore.ParameterSet.Config as cms + +from RecoMET.METProducers.PFMET_cfi import pfMet +tauPFMET = pfMet.clone() +tauPFMET.src = cms.InputTag('allDecayProducts') +tauPFMET.alias = cms.string('tauPFMET') + +tauDecayProducts = cms.EDProducer("tauDecayProducts", + src = cms.InputTag("slimmedTaus") ) + +allDecayProducts = cms.EDProducer("CandViewMerger", + src = cms.VInputTag( cms.InputTag("tauDecayProducts"), cms.InputTag("slimmedElectrons"), cms.InputTag("slimmedMuons"))) + +from PhysicsTools.PatAlgos.producersLayer1.metProducer_cfi import patMETs +tauMET = patMETs.clone() +tauMET.computeMETSignificance = cms.bool(True) +tauMET.addGenMET = cms.bool(False) +tauMET.metSource = cms.InputTag("tauPFMET") +tauMET.srcJets = cms.InputTag("slimmedJets") +tauMET.srcLeptons = cms.VInputTag("slimmedElectrons", "slimmedMuons", "slimmedPhotons") + +from RecoMET.METProducers.METSignificance_cfi import METSignificance + +tausSignificance = METSignificance.clone() +tausSignificance.srcMet = cms.InputTag('tauMET') +tausSignificance.srcPFCandidates = cms.InputTag('tauDecayProducts') diff --git a/RecoMET/METPUSubtraction/src/MvaMEtUtilities.cc b/RecoMET/METPUSubtraction/src/MvaMEtUtilities.cc deleted file mode 100644 index c31c3014e9299..0000000000000 --- a/RecoMET/METPUSubtraction/src/MvaMEtUtilities.cc +++ /dev/null @@ -1,284 +0,0 @@ -#include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h" - -#include "DataFormats/Math/interface/deltaR.h" - -#include -#include - -MvaMEtUtilities::MvaMEtUtilities(const edm::ParameterSet& cfg) -{ - // jet id - /* ===> this code parses an xml it uses parameter set - edm::ParameterSet jetConfig = iConfig.getParameter("JetIdParams"); - for(int i0 = 0; i0 < 3; i0++) { - std::string lCutType = "Tight"; - if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium"; - if(i0 == PileupJetIdentifier::kLoose) lCutType = "Loose"; - std::vector pt010 = jetConfig.getParameter >(("Pt010_" +lCutType).c_str()); - std::vector pt1020 = jetConfig.getParameter >(("Pt1020_"+lCutType).c_str()); - std::vector pt2030 = jetConfig.getParameter >(("Pt2030_"+lCutType).c_str()); - std::vector pt3050 = jetConfig.getParameter >(("Pt3050_"+lCutType).c_str()); - for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][0][i1] = pt010 [i1]; - for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][1][i1] = pt1020[i1]; - for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][2][i1] = pt2030[i1]; - for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][3][i1] = pt3050[i1]; - } - */ - //Tight Id => not used - mvaCut_[0][0][0] = 0.5; mvaCut_[0][0][1] = 0.6; mvaCut_[0][0][2] = 0.6; mvaCut_[0][0][3] = 0.9; - mvaCut_[0][1][0] = -0.2; mvaCut_[0][1][1] = 0.2; mvaCut_[0][1][2] = 0.2; mvaCut_[0][1][3] = 0.6; - mvaCut_[0][2][0] = 0.3; mvaCut_[0][2][1] = 0.4; mvaCut_[0][2][2] = 0.7; mvaCut_[0][2][3] = 0.8; - mvaCut_[0][3][0] = 0.5; mvaCut_[0][3][1] = 0.4; mvaCut_[0][3][2] = 0.8; mvaCut_[0][3][3] = 0.9; - //Medium id => not used - mvaCut_[1][0][0] = 0.2; mvaCut_[1][0][1] = 0.4; mvaCut_[1][0][2] = 0.2; mvaCut_[1][0][3] = 0.6; - mvaCut_[1][1][0] = -0.3; mvaCut_[1][1][1] = 0. ; mvaCut_[1][1][2] = 0. ; mvaCut_[1][1][3] = 0.5; - mvaCut_[1][2][0] = 0.2; mvaCut_[1][2][1] = 0.2; mvaCut_[1][2][2] = 0.5; mvaCut_[1][2][3] = 0.7; - mvaCut_[1][3][0] = 0.3; mvaCut_[1][3][1] = 0.2; mvaCut_[1][3][2] = 0.7; mvaCut_[1][3][3] = 0.8; - //Met Id => used - mvaCut_[2][0][0] = -0.2; mvaCut_[2][0][1] = -0.3; mvaCut_[2][0][2] = -0.5; mvaCut_[2][0][3] = -0.5; - mvaCut_[2][1][0] = -0.2; mvaCut_[2][1][1] = -0.2; mvaCut_[2][1][2] = -0.5; mvaCut_[2][1][3] = -0.3; - mvaCut_[2][2][0] = -0.2; mvaCut_[2][2][1] = -0.2; mvaCut_[2][2][2] = -0.2; mvaCut_[2][2][3] = 0.1; - mvaCut_[2][3][0] = -0.2; mvaCut_[2][3][1] = -0.2; mvaCut_[2][3][2] = 0. ; mvaCut_[2][3][3] = 0.2; - - dzCut_ = cfg.getParameter("dZcut"); - ptThreshold_ = ( cfg.exists("ptThreshold") ) ? - cfg.getParameter("ptThreshold") : -1000; -} - -MvaMEtUtilities::~MvaMEtUtilities() -{ -// nothing to be done yet... -} - -bool MvaMEtUtilities::passesMVA(const reco::Candidate::LorentzVector& jetP4, double mvaJetId) -{ - int ptBin = 0; - if ( jetP4.pt() >= 10. && jetP4.pt() < 20. ) ptBin = 1; - if ( jetP4.pt() >= 20. && jetP4.pt() < 30. ) ptBin = 2; - if ( jetP4.pt() >= 30. ) ptBin = 3; - - int etaBin = 0; - if ( std::abs(jetP4.eta()) >= 2.5 && std::abs(jetP4.eta()) < 2.75) etaBin = 1; - if ( std::abs(jetP4.eta()) >= 2.75 && std::abs(jetP4.eta()) < 3.0 ) etaBin = 2; - if ( std::abs(jetP4.eta()) >= 3.0 && std::abs(jetP4.eta()) < 5.0 ) etaBin = 3; - - return ( mvaJetId > mvaCut_[2][ptBin][etaBin] ); -} - -reco::Candidate::LorentzVector MvaMEtUtilities::leadJetP4(const std::vector& jets) -{ - return jetP4(jets, 0); -} - -reco::Candidate::LorentzVector MvaMEtUtilities::subleadJetP4(const std::vector& jets) -{ - return jetP4(jets, 1); -} - -reco::Candidate::LorentzVector MvaMEtUtilities::jetP4(const std::vector& jets, unsigned idx) -{ - reco::Candidate::LorentzVector retVal(0.,0.,0.,0.); - if ( idx < jets.size() ) { - std::vector jets_sorted = jets; - std::sort(jets_sorted.rbegin(), jets_sorted.rend()); - retVal = jets_sorted[idx].p4(); - } - return retVal; -} -unsigned MvaMEtUtilities::numJetsAboveThreshold(const std::vector& jets, double ptThreshold) -{ - unsigned retVal = 0; - for ( std::vector::const_iterator jet = jets.begin(); - jet != jets.end(); ++jet ) { - if ( jet->p4().pt() > ptThreshold ) ++retVal; - } - return retVal; -} -std::vector MvaMEtUtilities::cleanJets(const std::vector& jets, - const std::vector& leptons, - double ptThreshold, double dRmatch) -{ - - double dR2match = dRmatch*dRmatch; - std::vector retVal; - for ( std::vector::const_iterator jet = jets.begin(); - jet != jets.end(); ++jet ) { - bool isOverlap = false; - for ( std::vector::const_iterator lepton = leptons.begin(); - lepton != leptons.end(); ++lepton ) { - if ( deltaR2(jet->p4(), lepton->p4()) < dR2match ) isOverlap = true; - } - if ( jet->p4().pt() > ptThreshold && !isOverlap ) retVal.push_back(*jet); - } - return retVal; -} - -std::vector MvaMEtUtilities::cleanPFCands(const std::vector& pfCandidates, - const std::vector& leptons, - double dRmatch, bool invert) -{ - - double dR2match = dRmatch*dRmatch; - std::vector retVal; - for ( std::vector::const_iterator pfCandidate = pfCandidates.begin(); - pfCandidate != pfCandidates.end(); ++pfCandidate ) { - bool isOverlap = false; - for ( std::vector::const_iterator lepton = leptons.begin(); - lepton != leptons.end(); ++lepton ) { - if ( deltaR2(pfCandidate->p4(), lepton->p4()) < dR2match ) isOverlap = true; - } - if ( (!isOverlap && !invert) || (isOverlap && invert) ) retVal.push_back(*pfCandidate); - } - return retVal; -} -void MvaMEtUtilities::finalize(CommonMETData& metData) -{ - metData.met = sqrt(metData.mex*metData.mex + metData.mey*metData.mey); - metData.mez = 0.; - metData.phi = atan2(metData.mey, metData.mex); -} - -CommonMETData -MvaMEtUtilities::computeCandSum( int compKey, double dZmax, int dZflag, - bool iCharged, bool mvaPassFlag, - const std::vector& objects ) { - - CommonMETData retVal; - retVal.mex = 0.; - retVal.mey = 0.; - retVal.sumet = 0.; - - for ( std::vector::const_iterator object = objects.begin(); - object != objects.end(); ++object ) { - - double pFrac = 1; - - //pf candidates - // dZcut - // maximum distance within which tracks are - //considered to be associated to hard scatter vertex - // dZflag - // 0 : select charged PFCandidates originating from hard scatter vertex - // 1 : select charged PFCandidates originating from pile-up vertices - // 2 : select all PFCandidates - if( compKey==MvaMEtUtilities::kPFCands ) { - if ( object->dZ() < 0. && dZflag != 2 ) continue; - if ( object->dZ() > dZmax && dZflag == 0 ) continue; - if ( object->dZ() < dZmax && dZflag == 1 ) continue; - } - - //leptons - if( compKey==MvaMEtUtilities::kLeptons) { - if(iCharged) pFrac = object->chargedEnFrac(); - } - - //jets - if( compKey==MvaMEtUtilities::kJets) { - bool passesMVAjetId = passesMVA(object->p4(), object->mva() ); - - if ( passesMVAjetId && !mvaPassFlag ) continue; - if ( !passesMVAjetId && mvaPassFlag ) continue; - - pFrac = 1-object->chargedEnFrac();//neutral energy fraction - } - - retVal.mex += object->p4().px()*pFrac; - retVal.mey += object->p4().py()*pFrac; - retVal.sumet += object->p4().pt()*pFrac; - } - - finalize(retVal); - return retVal; -} - - -CommonMETData -MvaMEtUtilities::computeRecoil(int metType) { - - CommonMETData retVal; - - if(metType == MvaMEtUtilities::kPF ) { - //MET = pfMET = - all candidates - // MET (1) in JME-13-003 - retVal.mex = leptonsSum_.mex - pfCandSum_.mex; - retVal.mey = leptonsSum_.mey - pfCandSum_.mey; - retVal.sumet = pfCandSum_.sumet - leptonsSum_.sumet; - } - if(metType == MvaMEtUtilities::kChHS ) { - //MET = - charged HS - // MET (2) in JME-13-003 - retVal.mex = leptonsChSum_.mex - pfCandChHSSum_.mex; - retVal.mey = leptonsChSum_.mey - pfCandChHSSum_.mey; - retVal.sumet = pfCandChHSSum_.sumet - leptonsChSum_.sumet; - } - if(metType == MvaMEtUtilities::kHS ) { - //MET = - charged HS - neutral HS in jets - // MET (3) in JME-13-003 - retVal.mex = leptonsChSum_.mex - (pfCandChHSSum_.mex + neutralJetHSSum_.mex); - retVal.mey = leptonsChSum_.mey - (pfCandChHSSum_.mey + neutralJetHSSum_.mey); - retVal.sumet = pfCandChHSSum_.sumet + neutralJetHSSum_.sumet - leptonsChSum_.sumet; - } - if(metType == MvaMEtUtilities::kPU ) { - //MET = - charged PU - neutral PU in jets - //MET = -recoil in that particular case, - sign not useful for the MVA and then discarded - //motivated as PU IS its own recoil - // MET (4) in JME-13-003 - retVal.mex = -(pfCandChPUSum_.mex + neutralJetPUSum_.mex); - retVal.mey = -(pfCandChPUSum_.mey + neutralJetPUSum_.mey); - retVal.sumet = pfCandChPUSum_.sumet + neutralJetPUSum_.sumet; - } - if(metType == MvaMEtUtilities::kHSMinusNeutralPU ) { - //MET = all candidates - charged PU - neutral PU in jets - // = all charged HS + all neutrals - neutral PU in jets - // MET (5) in JME-13-003 - retVal.mex = leptonsSum_.mex - (pfCandSum_.mex - pfCandChPUSum_.mex - neutralJetPUSum_.mex); - retVal.mey = leptonsSum_.mey - (pfCandSum_.mey - pfCandChPUSum_.mey - neutralJetPUSum_.mey); - retVal.sumet = (pfCandSum_.sumet - pfCandChPUSum_.sumet - neutralJetPUSum_.sumet) -leptonsSum_.sumet; - } - - finalize(retVal); - return retVal; -} - -void -MvaMEtUtilities::computeAllSums(const std::vector& jets, - const std::vector& leptons, - const std::vector& pfCandidates ) { - - cleanedJets_ = cleanJets(jets, leptons, ptThreshold_, 0.5); - - leptonsSum_ = computeCandSum( kLeptons, 0., 0, false , false, leptons ); - leptonsChSum_ = computeCandSum( kLeptons, 0., 0, true , false, leptons); - pfCandSum_ = computeCandSum( kPFCands, dzCut_, 2, false , false, pfCandidates); - pfCandChHSSum_ = computeCandSum( kPFCands, dzCut_, 0, false , false, pfCandidates); - pfCandChPUSum_ = computeCandSum( kPFCands, dzCut_, 1, false , false, pfCandidates); - neutralJetHSSum_ = computeCandSum( kJets, 0., 0, false , true, jets ); - neutralJetPUSum_ = computeCandSum( kJets, 0., 0, false , false, jets ); - -} - -double -MvaMEtUtilities::getLeptonsSumMEX() const { - return leptonsSum_.mex; -} - -double -MvaMEtUtilities::getLeptonsSumMEY() const { - return leptonsSum_.mey; -} - -double -MvaMEtUtilities::getLeptonsChSumMEX() const { - return leptonsChSum_.mex; -} - -double -MvaMEtUtilities::getLeptonsChSumMEY() const { - return leptonsChSum_.mey; -} - - -const std::vector& -MvaMEtUtilities::getCleanedJets() const { - return cleanedJets_; -} diff --git a/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc b/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc deleted file mode 100644 index 4d6c5007e2556..0000000000000 --- a/RecoMET/METPUSubtraction/src/PFMETAlgorithmMVA.cc +++ /dev/null @@ -1,284 +0,0 @@ -#include "RecoMET/METPUSubtraction/interface/PFMETAlgorithmMVA.h" - -#include "FWCore/ParameterSet/interface/FileInPath.h" -#include "FWCore/Utilities/interface/Exception.h" - -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "CondFormats/DataRecord/interface/GBRWrapperRcd.h" - -#include "DataFormats/METReco/interface/CommonMETData.h" - -#include - -#include - -enum MVAType { kBaseline = 0 }; - -const double Pi=std::cos(-1); - -const std::string PFMETAlgorithmMVA::updateVariableNames(std::string input) -{ - if(input=="sumet") return "particleFlow_SumET"; - if(input=="npv") return "nPV"; - if(input=="pfu") return "particleFlow_U"; - if(input=="pfuphi") return "particleFlow_UPhi"; - if(input=="tksumet") return "track_SumET"; - if(input=="tku") return "track_U"; - if(input=="tkuphi") return "track_UPhi"; - if(input=="nopusumet") return "noPileUp_SumET"; - if(input=="nopuu") return "noPileUp_U"; - if(input=="nopuuphi") return "noPileUp_UPhi"; - if(input=="pusumet") return "pileUp_SumET"; - if(input=="pumet") return "pileUp_MET"; - if(input=="pumetphi") return "pileUp_METPhi"; - if(input=="pucsumet") return "pileUpCorrected_SumET"; - if(input=="pucu") return "pileUpCorrected_U"; - if(input=="pucuphi") return "pileUpCorrected_UPhi"; - if(input=="jetpt1") return "jet1_pT"; - if(input=="jeteta1") return "jet1_eta"; - if(input=="jetphi1") return "jet1_Phi"; - if(input=="jetpt2") return "jet2_pT"; - if(input=="jeteta2") return "jet2_eta"; - if(input=="jetphi2") return "jet2_Phi"; - if(input=="nalljet") return "nJets"; - if(input=="njet") return "numJetsPtGt30"; - if(input=="uphi_mva") return "PhiCor_UPhi"; - if(input=="uphix_mva") return "PhiCor_UPhi"; - if(input=="ux_mva") return "RecoilCor_U"; - return input; -} - -const GBRForest* PFMETAlgorithmMVA::loadMVAfromFile(const edm::FileInPath& inputFileName, const std::string& mvaName) -{ - if ( inputFileName.location()==edm::FileInPath::Unknown ) throw cms::Exception("PFMETAlgorithmMVA::loadMVA") - << " Failed to find File = " << inputFileName << " !!\n"; - std::unique_ptr inputFile(new TFile(inputFileName.fullPath().data()) ); - - std::vector *lVec = (std::vector*)inputFile->Get("varlist"); - - if(lVec==nullptr) { - throw cms::Exception("PFMETAlgorithmMVA::loadMVA") - << " Failed to load mva file : " << inputFileName.fullPath().data() << " is not a proper file !!\n"; - } - - std::vector variableNames; - for(unsigned int i=0; i< lVec->size();++i) - { - variableNames.push_back(updateVariableNames(lVec->at(i))); - } - - if(mvaName.find(mvaNameU_) != std::string::npos) varForU_ = variableNames; - else if(mvaName.find(mvaNameDPhi_) != std::string::npos) varForDPhi_ = variableNames; - else if(mvaName.find(mvaNameCovU1_) != std::string::npos) varForCovU1_ = variableNames; - else if(mvaName.find(mvaNameCovU2_) != std::string::npos) varForCovU2_ = variableNames; - else throw cms::Exception("PFMETAlgorithmMVA::loadMVA") << "MVA MET weight file tree names do not match specified inputs" << std::endl; - - - const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data()); - if ( !mva ) - throw cms::Exception("PFMETAlgorithmMVA::loadMVA") - << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n"; - - return mva; -} - -const GBRForest* PFMETAlgorithmMVA::loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName) -{ - edm::ESHandle mva; - es.get().get(mvaName, mva); - return mva.product(); -} - -PFMETAlgorithmMVA::PFMETAlgorithmMVA(const edm::ParameterSet& cfg) - : utils_(cfg), - mvaReaderU_(nullptr), - mvaReaderDPhi_(nullptr), - mvaReaderCovU1_(nullptr), - mvaReaderCovU2_(nullptr), - cfg_(cfg) -{ - mvaType_ = kBaseline; - - loadMVAfromDB_ = cfg.getParameter("loadMVAfromDB"); - -} - -PFMETAlgorithmMVA::~PFMETAlgorithmMVA() -{ - if ( !loadMVAfromDB_ ) { - delete mvaReaderU_; - delete mvaReaderDPhi_; - delete mvaReaderCovU1_; - delete mvaReaderCovU2_; - } -} - -//------------------------------------------------------------------------------- -void PFMETAlgorithmMVA::initialize(const edm::EventSetup& es) -{ - edm::ParameterSet cfgInputRecords = cfg_.getParameter("inputRecords"); - mvaNameU_ = cfgInputRecords.getParameter("U"); - mvaNameDPhi_ = cfgInputRecords.getParameter("DPhi"); - mvaNameCovU1_ = cfgInputRecords.getParameter("CovU1"); - mvaNameCovU2_ = cfgInputRecords.getParameter("CovU2"); - - if ( loadMVAfromDB_ ) { - mvaReaderU_ = loadMVAfromDB(es, mvaNameU_); - mvaReaderDPhi_ = loadMVAfromDB(es, mvaNameDPhi_); - mvaReaderCovU1_ = loadMVAfromDB(es, mvaNameCovU1_); - mvaReaderCovU2_ = loadMVAfromDB(es, mvaNameCovU2_); - } else { - edm::ParameterSet cfgInputFileNames = cfg_.getParameter("inputFileNames"); - - edm::FileInPath inputFileNameU = cfgInputFileNames.getParameter("U"); - mvaReaderU_ = loadMVAfromFile(inputFileNameU, mvaNameU_); - edm::FileInPath inputFileNameDPhi = cfgInputFileNames.getParameter("DPhi"); - mvaReaderDPhi_ = loadMVAfromFile(inputFileNameDPhi, mvaNameDPhi_); - edm::FileInPath inputFileNameCovU1 = cfgInputFileNames.getParameter("CovU1"); - mvaReaderCovU1_ = loadMVAfromFile(inputFileNameCovU1, mvaNameCovU1_); - edm::FileInPath inputFileNameCovU2 = cfgInputFileNames.getParameter("CovU2"); - mvaReaderCovU2_ = loadMVAfromFile(inputFileNameCovU2, mvaNameCovU2_); - } -} - -//------------------------------------------------------------------------------- -void PFMETAlgorithmMVA::setInput(const std::vector& leptons, - const std::vector& jets, - const std::vector& pfCandidates, - const std::vector& vertices) -{ - - - utils_.computeAllSums( jets, leptons, pfCandidates); - - sumLeptonPx_ = utils_.getLeptonsSumMEX(); - sumLeptonPy_ = utils_.getLeptonsSumMEY(); - - chargedSumLeptonPx_ = utils_.getLeptonsChSumMEX(); - chargedSumLeptonPy_ = utils_.getLeptonsChSumMEY(); - - const std::vector jets_cleaned = utils_.getCleanedJets(); - - CommonMETData pfRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kPF ); - CommonMETData chHSRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kChHS ); - CommonMETData hsRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kHS ); - CommonMETData puRecoil_data = utils_.computeRecoil( MvaMEtUtilities::kPU ); - CommonMETData hsMinusNeutralPUMEt_data = utils_.computeRecoil( MvaMEtUtilities::kHSMinusNeutralPU ); - - reco::Candidate::LorentzVector jet1P4 = utils_.leadJetP4(jets_cleaned); - reco::Candidate::LorentzVector jet2P4 = utils_.subleadJetP4(jets_cleaned); - - - var_["particleFlow_U"] = pfRecoil_data.met; - var_["particleFlow_SumET"] = pfRecoil_data.sumet; - var_["particleFlow_UPhi"] = pfRecoil_data.phi; - - var_["track_SumET"] = chHSRecoil_data.sumet/var_["particleFlow_SumET"]; - var_["track_U"] = chHSRecoil_data.met; - var_["track_UPhi"] = chHSRecoil_data.phi; - - var_["noPileUp_SumET"] = hsRecoil_data.sumet/var_["particleFlow_SumET"]; - var_["noPileUp_U"] = hsRecoil_data.met; - var_["noPileUp_UPhi"] = hsRecoil_data.phi; - - var_["pileUp_SumET"] = puRecoil_data.sumet/var_["particleFlow_SumET"]; - var_["pileUp_MET"] = puRecoil_data.met; - var_["pileUp_METPhi"] = puRecoil_data.phi; - - var_["pileUpCorrected_SumET"] = hsMinusNeutralPUMEt_data.sumet/var_["particleFlow_SumET"]; - var_["pileUpCorrected_U"] = hsMinusNeutralPUMEt_data.met; - var_["pileUpCorrected_UPhi"] = hsMinusNeutralPUMEt_data.phi; - - var_["jet1_pT"] = jet1P4.pt(); - var_["jet1_eta"] = jet1P4.eta(); - var_["jet1_Phi"] = jet1P4.phi(); - var_["jet2_pT"] = jet2P4.pt(); - var_["jet2_eta"] = jet2P4.eta(); - var_["jet2_Phi"] = jet2P4.phi(); - - var_["numJetsPtGt30"] = utils_.numJetsAboveThreshold(jets_cleaned, 30.); - var_["nJets"] = jets_cleaned.size(); - var_["nPV"] = vertices.size(); -} - -//------------------------------------------------------------------------------- -std::unique_ptr PFMETAlgorithmMVA::createFloatVector(std::vector variableNames) -{ - std::unique_ptr floatVector(new float[variableNames.size()]); - int i = 0; - for(auto variableName: variableNames) - { - floatVector[i++] = var_[variableName]; - } - return floatVector; -} - -//------------------------------------------------------------------------------- -void PFMETAlgorithmMVA::evaluateMVA() -{ - // CV: MVAs needs to be evaluated in order { DPhi, U1, CovU1, CovU2 } - // as MVA for U1 (CovU1, CovU2) uses output of DPhi (DPhi and U1) MVA - mvaOutputDPhi_ = GetResponse(mvaReaderDPhi_, varForDPhi_); - var_["PhiCor_UPhi"] = var_["particleFlow_UPhi"] + mvaOutputDPhi_; - mvaOutputU_ = GetResponse(mvaReaderU_, varForU_); - var_["RecoilCor_U"] = var_["particleFlow_U"] * mvaOutputU_; - var_["RecoilCor_UPhi"] = var_["PhiCor_UPhi"]; - mvaOutputCovU1_ = std::pow(GetResponse(mvaReaderCovU1_, varForCovU1_)* mvaOutputU_ * var_["particleFlow_U"], 2); - mvaOutputCovU2_ = std::pow(GetResponse(mvaReaderCovU2_, varForCovU2_)* mvaOutputU_ * var_["particleFlow_U"], 2); - - - // compute MET(Photon check) - if(hasPhotons_) { - //Fix events with unphysical properties - double sumLeptonPt = std::max(sqrt(sumLeptonPx_*sumLeptonPx_+sumLeptonPy_*sumLeptonPy_),1.); - if(var_["track_U"]/sumLeptonPt < 0.1 || var_["noPileUp_U"]/sumLeptonPt < 0.1 ) { - mvaOutputU_ = 1.; - mvaOutputDPhi_ = 0.; - } - } - computeMET(); -} -//------------------------------------------------------------------------------- - -void PFMETAlgorithmMVA::computeMET() -{ - double U = var_["RecoilCor_U"]; - double Phi = var_["PhiCor_UPhi"]; - if ( U < 0. ) Phi += Pi; //RF: No sign flip for U necessary in that case? - double cosPhi = std::cos(Phi); - double sinPhi = std::sin(Phi); - double metPx = U*cosPhi - sumLeptonPx_; - double metPy = U*sinPhi - sumLeptonPy_; - double metPt = sqrt(metPx*metPx + metPy*metPy); - mvaMEt_.SetCoordinates(metPx, metPy, 0., metPt); - // compute MET uncertainties in dirrections parallel and perpendicular to hadronic recoil - // (neglecting uncertainties on lepton momenta) - mvaMEtCov_(0, 0) = mvaOutputCovU1_*cosPhi*cosPhi + mvaOutputCovU2_*sinPhi*sinPhi; - mvaMEtCov_(0, 1) = -mvaOutputCovU1_*sinPhi*cosPhi + mvaOutputCovU2_*sinPhi*cosPhi; - mvaMEtCov_(1, 0) = mvaMEtCov_(0, 1); - mvaMEtCov_(1, 1) = mvaOutputCovU1_*sinPhi*sinPhi + mvaOutputCovU2_*cosPhi*cosPhi; -} - -//------------------------------------------------------------------------------- -const float PFMETAlgorithmMVA::GetResponse(const GBRForest * Reader, std::vector &variableNames ) -{ - std::unique_ptr mvaInputVector = createFloatVector(variableNames); - double result = Reader->GetResponse( mvaInputVector.get() ); - return result; -} - -//------------------------------------------------------------------------------- -void PFMETAlgorithmMVA::print(std::ostream& stream) const -{ - stream << ":" << std::endl; - for(auto entry: var_) - stream << entry.first << " = " << entry.second << std::endl; - stream << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl; - stream << " sum(leptons): Pt = " << sqrt(sumLeptonPx_*sumLeptonPx_ + sumLeptonPy_*sumLeptonPy_) << "," - << " phi = " << atan2(sumLeptonPy_, sumLeptonPx_) << " " - << "(Px = " << sumLeptonPx_ << ", Py = " << sumLeptonPy_ << ")" ; - stream << " MVA output: U = " << mvaOutputU_ << ", dPhi = " << mvaOutputDPhi_ << "," << " covU1 = " << mvaOutputCovU1_ << ", covU2 = " << mvaOutputCovU2_ << std::endl; - stream << std::endl; -} - diff --git a/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py b/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py deleted file mode 100644 index c86b0e9a01c2b..0000000000000 --- a/RecoMET/METPUSubtraction/test/mvaMETOnMiniAOD_cfg.py +++ /dev/null @@ -1,83 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -process = cms.Process('MVAMET') - -# import of standard configurations -process.load('Configuration.StandardSequences.Services_cff') -process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi') -process.load('FWCore.MessageService.MessageLogger_cfi') -process.load('Configuration.EventContent.EventContent_cff') -process.load('SimGeneral.MixingModule.mixNoPU_cfi') -process.load('Configuration.StandardSequences.GeometryRecoDB_cff') -process.load('Configuration.StandardSequences.MagneticField_38T_cff') -process.load('Configuration.StandardSequences.EndOfProcess_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') - - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(1000) -) - - - -# Input source -process.source = cms.Source("PoolSource", - secondaryFileNames = cms.untracked.vstring(), - fileNames = cms.untracked.vstring('/store/mc/RunIIFall15MiniAODv2/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-madgraphMLM-pythia8/MINIAODSIM/PU25nsData2015v1_76X_mcRun2_asymptotic_v12-v1/70000/DC327D28-D3B8-E511-9FB7-008CFA0A59C0.root '), - skipEvents = cms.untracked.uint32(0) -) - - -process.options = cms.untracked.PSet( - allowUnscheduled = cms.untracked.bool(True) -) - - -# Output definition - -process.MINIAODSIMoutput = cms.OutputModule("PoolOutputModule", - compressionLevel = cms.untracked.int32(4), - compressionAlgorithm = cms.untracked.string('LZMA'), - eventAutoFlushCompressedSize = cms.untracked.int32(15728640), - outputCommands = cms.untracked.vstring( "keep *_ak4PFJets_*_MVAMET", - "keep *_pfMVAMEt_*_MVAMET", - "keep recoMET_*_*_*", - "keep patMETs_*_*_*", - ), - fileName = cms.untracked.string('miniAODMVAMET.root'), - dataset = cms.untracked.PSet( - filterName = cms.untracked.string(''), - dataTier = cms.untracked.string('') - ), - dropMetaData = cms.untracked.string('ALL'), - fastCloning = cms.untracked.bool(False), - overrideInputFileSplitLevels = cms.untracked.bool(True) -) - - - -process.load("RecoJets.JetProducers.ak4PFJets_cfi") -process.ak4PFJets.src = cms.InputTag("packedPFCandidates") -process.ak4PFJets.doAreaFastjet = cms.bool(True) - -from JetMETCorrections.Configuration.DefaultJEC_cff import ak4PFJetsL1FastL2L3 - -process.load("RecoMET.METPUSubtraction.mvaPFMET_cff") -#process.pfMVAMEt.srcLeptons = cms.VInputTag("slimmedElectrons") -process.pfMVAMEt.srcPFCandidates = cms.InputTag("packedPFCandidates") -process.pfMVAMEt.srcVertices = cms.InputTag("offlineSlimmedPrimaryVertices") - -process.puJetIdForPFMVAMEt.jec = cms.string('AK4PF') -#process.puJetIdForPFMVAMEt.jets = cms.InputTag("ak4PFJets") -process.puJetIdForPFMVAMEt.vertexes = cms.InputTag("offlineSlimmedPrimaryVertices") -process.puJetIdForPFMVAMEt.rho = cms.InputTag("fixedGridRhoFastjetAll") - - - -# Other statements -from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '') - -# Path and EndPath definitions -process.endjob_step = cms.EndPath(process.endOfProcess) -process.MINIAODSIMoutput_step = cms.EndPath(process.MINIAODSIMoutput) diff --git a/RecoMET/METPUSubtraction/test/runMVAMET.py b/RecoMET/METPUSubtraction/test/runMVAMET.py new file mode 100644 index 0000000000000..8efaea7556d94 --- /dev/null +++ b/RecoMET/METPUSubtraction/test/runMVAMET.py @@ -0,0 +1,67 @@ +import sys +import FWCore.ParameterSet.Config as cms +from FWCore.ParameterSet.VarParsing import VarParsing +from PhysicsTools.PatAlgos.tools.tauTools import * +from RecoMET.METPUSubtraction.MVAMETConfiguration_cff import runMVAMET + +options = VarParsing ('python') +options.register ('globalTag',"80X_mcRun2_asymptotic_2016_miniAODv2_v1",VarParsing.multiplicity.singleton,VarParsing.varType.string,'input global tag to be used'); +options.register ('inputFile', 'file:////pnfs/desy.de/cms/tier2/store/mc/RunIISpring16MiniAODv2/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8/MINIAODSIM/PUSpring16RAWAODSIM_reHLT_80X_mcRun2_asymptotic_v14-v1/40000/00E9D1DA-105D-E611-A56E-FA163EE988CA.root', VarParsing.multiplicity.singleton, VarParsing.varType.string, "Path to a testfile") +options.register ("localSqlite", '', VarParsing.multiplicity.singleton, VarParsing.varType.string, "Path to a local sqlite file") +options.register ("reapplyJEC", False, VarParsing.multiplicity.singleton, VarParsing.varType.bool, "Reapply JEC to Jets") +options.register ("reapplyPUJetID", False, VarParsing.multiplicity.singleton, VarParsing.varType.bool, "Reapply PU Jet ID") +options.register ("isData", False, VarParsing.multiplicity.singleton, VarParsing.varType.bool, "Input is data") +options.parseArguments() + +process = cms.Process("MVAMET") + +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_38T_cff') +jetCollection = "slimmedJets" + +if (options.localSqlite == ''): + process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff") + process.GlobalTag.globaltag = options.globalTag +else: + from RecoMET.METPUSubtraction.jet_recorrections import loadLocalSqlite + loadLocalSqlite(process, options.localSqlite) + +if options.reapplyPUJetID: + from RecoMET.METPUSubtraction.jet_recorrections import reapplyPUJetID + reapplyPUJetID(process) + +if options.reapplyJEC: + from RecoMET.METPUSubtraction.jet_recorrections import recorrectJets + recorrectJets(process, options.isData) + jetCollection = "patJetsReapplyJEC" + +if options.reapplyPUJetID: + getattr(process, jetCollection).userData.userFloats.src += ['pileupJetIdUpdated:fullDiscriminant'] + +# configure MVA MET +runMVAMET( process, jetCollectionPF = jetCollection) + +## set input files +process.source = cms.Source("PoolSource") +process.source.fileNames = cms.untracked.vstring(options.inputFile) +## logger +process.load('FWCore.MessageLogger.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 50 + +#! Output and Log +process.options = cms.untracked.PSet(wantSummary = cms.untracked.bool(True)) +process.options.allowUnscheduled = cms.untracked.bool(True) + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(options.maxEvents) +) + +process.output = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('output.root'), + outputCommands = cms.untracked.vstring( + 'keep patMETs_MVAMET_MVAMET_MVAMET', + 'keep *_patJetsReapplyJEC_*_MVAMET' + ), + SelectEvents = cms.untracked.PSet( SelectEvents = cms.vstring('p')) + ) +process.out = cms.EndPath(process.output) diff --git a/RecoMET/METPUSubtraction/test/testMVAMetProducer.py b/RecoMET/METPUSubtraction/test/testMVAMetProducer.py deleted file mode 100644 index 4f4b7f4fff56b..0000000000000 --- a/RecoMET/METPUSubtraction/test/testMVAMetProducer.py +++ /dev/null @@ -1,32 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -process = cms.Process("analysis") -process.load("FWCore.MessageService.MessageLogger_cfi") -process.MessageLogger.cerr.FwkReport.reportEvery = cms.untracked.int32(10) -process.load('Configuration.StandardSequences.Services_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') -process.load('JetMETCorrections.Configuration.JetCorrectionProducers_cff') -process.load('RecoMET.METPUSubtraction.mvaPFMET_cff') - -# Other statements -from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:run2_mc', '') - -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(20) ) - -process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( -"root://eoscms//eos/cms/store/relval/CMSSW_7_5_0_pre4/RelValZMM_13/GEN-SIM-RECO/PU50ns_MCRUN2_75_V0-v1/00000/24B4BCD7-22F8-E411-A740-0025905A48F2.root" - ), - skipEvents = cms.untracked.uint32(0) -) - -process.output = cms.OutputModule("PoolOutputModule", - outputCommands = cms.untracked.vstring('keep *'), - fileName = cms.untracked.string("MVaTest.root") -) - -process.ana = cms.Sequence(process.pfMVAMEtSequence) -process.p = cms.Path(process.ana) -process.outpath = cms.EndPath(process.output) diff --git a/RecoMET/METPUSubtraction/test/writeGBRForests_cfg.py b/RecoMET/METPUSubtraction/test/writeGBRForests_cfg.py deleted file mode 100644 index 4a92d10bd4882..0000000000000 --- a/RecoMET/METPUSubtraction/test/writeGBRForests_cfg.py +++ /dev/null @@ -1,72 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -process = cms.Process("writeGBRForests") - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(1) # CV: needs to be set to 1 so that GBRForestWriter::analyze method gets called exactly once -) - -process.source = cms.Source("EmptySource") - -process.load('Configuration/StandardSequences/Services_cff') - -process.gbrForestWriter = cms.EDAnalyzer("GBRForestWriter", - jobs = cms.VPSet( - cms.PSet( - inputFileName = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmet_53_Dec2012.root'), - inputFileType = cms.string("GBRForest"), - gbrForestName = cms.string("U1Correction"), - outputFileType = cms.string("SQLLite"), - outputRecord = cms.string("mvaPFMET_53_Dec2012_U") - ), - cms.PSet( - inputFileName = cms.FileInPath('RecoMET/METPUSubtraction/data/gbrmetphi_53_Dec2012.root'), - inputFileType = cms.string("GBRForest"), - gbrForestName = cms.string("PhiCorrection"), - outputFileType = cms.string("SQLLite"), - outputRecord = cms.string("mvaPFMET_53_Dec2012_DPhi") - ), - cms.PSet( - inputFileName = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru1cov_53_Dec2012.root'), - inputFileType = cms.string("GBRForest"), - gbrForestName = cms.string("CovU1"), - outputFileType = cms.string("SQLLite"), - outputRecord = cms.string("mvaPFMET_53_Dec2012_CovU1") - ), - cms.PSet( - inputFileName = cms.FileInPath('RecoMET/METPUSubtraction/data/gbru2cov_53_Dec2012.root'), - inputFileType = cms.string("GBRForest"), - gbrForestName = cms.string("CovU2"), - outputFileType = cms.string("SQLLite"), - outputRecord = cms.string("mvaPFMET_53_Dec2012_CovU2") - ) - ) -) - -process.load("CondCore.DBCommon.CondDBCommon_cfi") -process.CondDBCommon.connect = 'sqlite_file:../data/mvaPFMEt_53_Dec2012.db' - -process.PoolDBOutputService = cms.Service("PoolDBOutputService", - process.CondDBCommon, - timetype = cms.untracked.string('runnumber'), - toPut = cms.VPSet( - cms.PSet( - record = cms.string('mvaPFMET_53_Dec2012_U'), - tag = cms.string('mvaPFMET_53_Dec2012_U') - ), - cms.PSet( - record = cms.string('mvaPFMET_53_Dec2012_DPhi'), - tag = cms.string('mvaPFMET_53_Dec2012_DPhi') - ), - cms.PSet( - record = cms.string('mvaPFMET_53_Dec2012_CovU1'), - tag = cms.string('mvaPFMET_53_Dec2012_CovU1') - ), - cms.PSet( - record = cms.string('mvaPFMET_53_Dec2012_CovU2'), - tag = cms.string('mvaPFMET_53_Dec2012_CovU2') - ) - ) -) - -process.p = cms.Path(process.gbrForestWriter)