Skip to content

Commit

Permalink
Merge pull request #13115 from cms-met/METCorDev_Jan16_METSig
Browse files Browse the repository at this point in the history
MET PAT tool update + update of MET significance
  • Loading branch information
davidlange6 committed Mar 17, 2016
2 parents 5f85d4a + 1f1c606 commit 6c4b6cd
Show file tree
Hide file tree
Showing 40 changed files with 1,395 additions and 1,268 deletions.
20 changes: 9 additions & 11 deletions CommonTools/PileupAlgos/plugins/PuppiPhoton.cc
Expand Up @@ -117,19 +117,17 @@ void PuppiPhoton::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
if(matchPFCandidate(&(*itPF),*itPho)) pWeight = weight_;
}
} else {
int iPho = -1;
int iPho = -1;
for(std::vector<uint16_t>::const_iterator itPho = phoIndx.begin(); itPho!=phoIndx.end(); itPho++) {
iPho++;
if(pupCol->refAt(iPF).key() != *itPho) continue;
pWeight = weight_;
if(!useValueMap_) {
double pCorr = phoCands[iPho]->pt()/itPF->pt();
pWeight = pWeight*pCorr;
}
foundPhoIndex.push_back(iPho);
}
iPho++;
if(pupCol->refAt(iPF).key() != *itPho) continue;
pWeight = weight_;
if(!useValueMap_ && itPF->pt() != 0) pWeight = pWeight*(phoCands[iPho]->pt()/itPF->pt());
foundPhoIndex.push_back(iPho); }
}
pVec.SetPxPyPzE(itPF->px()*pWeight,itPF->py()*pWeight,itPF->pz()*pWeight,itPF->energy()*pWeight);

if(itPF->pt() != 0) pVec.SetPxPyPzE(itPF->px()*pWeight,itPF->py()*pWeight,itPF->pz()*pWeight,itPF->energy()*pWeight);

lWeights.push_back(pWeight);
pCand.setP4(pVec);
puppiP4s.push_back( pVec );
Expand Down
1 change: 0 additions & 1 deletion DataFormats/PatCandidates/src/MET.cc
Expand Up @@ -67,7 +67,6 @@ uncertainties_(iOther.uncertainties_),
corrections_(iOther.corrections_),
caloPackedMet_(iOther.caloPackedMet_) {

metSig_ =0.;
initCorMap();
}

Expand Down
File renamed without changes.
222 changes: 222 additions & 0 deletions JetMETCorrections/Type1MET/interface/JetCleanerForType1METT.h
@@ -0,0 +1,222 @@
#ifndef JetMETCorrections_Type1MET_JetCleanerForType1METT_h
#define JetMETCorrections_Type1MET_JetCleanerForType1METT_h

/** \class JetCleanerForType1METT
*
* Clean jets for MET corrections and uncertainties (pt/eta/EM fraction and muons)
* apply also JECs
*
* NOTE: class is templated to that it works with reco::PFJets as well as with pat::Jets of PF-type as input
*
* \authors Matthieu Marionneau ETH
*
*
*
*/

#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 "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/Common/interface/RefToBase.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "DataFormats/METReco/interface/CorrMETData.h"
#include "DataFormats/PatCandidates/interface/Jet.h"

#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
#include "DataFormats/MuonReco/interface/Muon.h"

#include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"

#include <string>
#include <type_traits>


namespace JetCleanerForType1MET_namespace
{
template <typename T, typename Textractor>
class InputTypeCheckerT
{
public:

void operator()(const T&) const {} // no type-checking needed for reco::PFJet input
bool isPatJet(const T&) const {
return std::is_base_of<class pat::Jet, T>::value;
}
};

template <typename T>
class RawJetExtractorT // this template is neccessary to support pat::Jets
// (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
// But it does not handle the muon removal!!!!! MM
{
public:

RawJetExtractorT(){}
reco::Candidate::LorentzVector operator()(const T& jet) const
{
return jet.p4();
}
};
}


template <typename T, typename Textractor>
class JetCleanerForType1METT : public edm::stream::EDProducer<>
{
public:

explicit JetCleanerForType1METT(const edm::ParameterSet& cfg):
moduleLabel_(cfg.getParameter<std::string>("@module_label")),
offsetCorrLabel_(""),
skipMuonSelection_(nullptr)
{
token_ = consumes<std::vector<T> >(cfg.getParameter<edm::InputTag>("src"));

offsetCorrLabel_ = cfg.getParameter<edm::InputTag>("offsetCorrLabel");
offsetCorrToken_ = consumes<reco::JetCorrector>(offsetCorrLabel_);

jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel"); //for MC
jetCorrLabelRes_ = cfg.getParameter<edm::InputTag>("jetCorrLabelRes"); //for data
jetCorrToken_ = mayConsume<reco::JetCorrector>(jetCorrLabel_);
jetCorrResToken_ = mayConsume<reco::JetCorrector>(jetCorrLabelRes_);

jetCorrEtaMax_ = cfg.getParameter<double>("jetCorrEtaMax");

type1JetPtThreshold_ = cfg.getParameter<double>("type1JetPtThreshold");

skipEM_ = cfg.getParameter<bool>("skipEM");
if ( skipEM_ ) {
skipEMfractionThreshold_ = cfg.getParameter<double>("skipEMfractionThreshold");
}

skipMuons_ = cfg.getParameter<bool>("skipMuons");
if ( skipMuons_ ) {
std::string skipMuonSelection_string = cfg.getParameter<std::string>("skipMuonSelection");
skipMuonSelection_.reset( new StringCutObjectSelector<reco::Candidate>(skipMuonSelection_string,true) ) ;
}

produces<std::vector<T> >();
}

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<std::string>("@module_label");
desc.add<edm::InputTag>("src");
desc.add<edm::InputTag>("offsetCorrLabel");
desc.add<edm::InputTag>("jetCorrLabel");
desc.add<edm::InputTag>("jetCorrLabelRes");
desc.add<double>("jetCorrEtaMax",9.9);
desc.add<double>("type1JetPtThreshold");
desc.add<bool>("skipEM");
desc.add<double>("skipEMfractionThreshold");
desc.add<bool>("skipMuons");
desc.add<std::string>("skipMuonSelection");
descriptions.addDefault(desc);
}

private:

void produce(edm::Event& evt, const edm::EventSetup& es)
{

edm::Handle<reco::JetCorrector> jetCorr;
//automatic switch for residual corrections
if(evt.isRealData() ) {
jetCorrLabel_ = jetCorrLabelRes_;
evt.getByToken(jetCorrResToken_, jetCorr);
} else {
evt.getByToken(jetCorrToken_, jetCorr);
}

typedef std::vector<T> JetCollection;
edm::Handle<JetCollection> jets;
evt.getByToken(token_, jets);

//new collection
std::auto_ptr< std::vector<T> > cleanedJets( new std::vector<T>() );

int numJets = jets->size();
for(int jetIndex=0;jetIndex<numJets; ++jetIndex ) {
const T& jet = jets->at(jetIndex);

const static JetCleanerForType1MET_namespace::InputTypeCheckerT<T, Textractor> checkInputType {};
checkInputType(jet);

double emEnergyFraction = jet.chargedEmEnergyFraction() + jet.neutralEmEnergyFraction();
if(skipEM_&&emEnergyFraction>skipEMfractionThreshold_ ) continue;

const static JetCleanerForType1MET_namespace::RawJetExtractorT<T> rawJetExtractor {};
reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);

if ( skipMuons_ ) {
const std::vector<reco::CandidatePtr> & cands = jet.daughterPtrVector();
for ( std::vector<reco::CandidatePtr>::const_iterator cand = cands.begin();
cand != cands.end(); ++cand ) {
const reco::PFCandidate *pfcand = dynamic_cast<const reco::PFCandidate *>(cand->get());
const reco::Candidate *mu = (pfcand != 0 ? ( pfcand->muonRef().isNonnull() ? pfcand->muonRef().get() : 0) : cand->get());
if ( mu != 0 && (*skipMuonSelection_)(*mu) ) {
reco::Candidate::LorentzVector muonP4 = (*cand)->p4();
rawJetP4 -= muonP4;
}
}
}

reco::Candidate::LorentzVector corrJetP4;
if ( checkInputType.isPatJet(jet) ) {
corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4);
}
else {
corrJetP4 = jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
if(corrJetP4.pt()<type1JetPtThreshold_) continue;
}

if(corrJetP4.pt()<type1JetPtThreshold_) continue;

cleanedJets->push_back(jet);
}

evt.put(cleanedJets);
}



std::string moduleLabel_;

edm::EDGetTokenT<std::vector<T> > token_;

edm::InputTag offsetCorrLabel_;
edm::EDGetTokenT<reco::JetCorrector> offsetCorrToken_; // e.g. 'ak5CaloJetL1Fastjet'
edm::InputTag jetCorrLabel_;
edm::InputTag jetCorrLabelRes_;
edm::EDGetTokenT<reco::JetCorrector> jetCorrToken_; // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
edm::EDGetTokenT<reco::JetCorrector> jetCorrResToken_; // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
Textractor jetCorrExtractor_;

double jetCorrEtaMax_; // do not use JEC factors for |eta| above this threshold

double type1JetPtThreshold_; // threshold to distinguish between jets entering Type 1 MET correction
// and jets entering "unclustered energy" sum
// NOTE: threshold is applied on **corrected** jet energy (recommended default = 10 GeV)

bool skipEM_; // flag to exclude jets with large fraction of electromagnetic energy (electrons/photons)
// from Type 1 + 2 MET corrections
double skipEMfractionThreshold_;

bool skipMuons_; // flag to subtract momentum of muons (provided muons pass selection cuts) which are within jets
// from jet energy before compute JECs/propagating JECs to Type 1 + 2 MET corrections
std::unique_ptr<StringCutObjectSelector< reco::Candidate> > skipMuonSelection_;


};

#endif

Expand Up @@ -65,6 +65,7 @@ namespace PFJetMETcorrInputProducer_namespace
return jet.p4();
}
};

}

template <typename T, typename Textractor>
Expand Down
13 changes: 13 additions & 0 deletions JetMETCorrections/Type1MET/plugins/JetCleanerForType1MET.cc
@@ -0,0 +1,13 @@
#include "JetMETCorrections/Type1MET/interface/JetCleanerForType1METT.h"

#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"

#include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"

typedef JetCleanerForType1METT<reco::PFJet, JetCorrExtractorT<reco::PFJet> > PFJetCleanerForType1MET;

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(PFJetCleanerForType1MET);

11 changes: 6 additions & 5 deletions PhysicsTools/PatAlgos/plugins/PATJetUpdater.cc
Expand Up @@ -130,9 +130,10 @@ void PATJetUpdater::produce(edm::Event & iEvent, const edm::EventSetup & iSetup)
Jet ajet( edm::RefToBase<Jet>(jetRef.castTo<JetRef>()) );

if (addJetCorrFactors_) {
unsigned int setindex = ajet.availableJECSets().size();
// Undo previous jet energy corrections
// undo previous jet energy corrections
ajet.setP4(ajet.correctedP4(0));
// clear previous JetCorrFactors
ajet.jec_.clear();
// add additional JetCorrs to the jet
for ( unsigned int i=0; i<jetCorrFactorsTokens_.size(); ++i ) {
const JetCorrFactors& jcf = jetCorrs[i][jetRef];
Expand All @@ -142,13 +143,13 @@ void PATJetUpdater::produce(edm::Event & iEvent, const edm::EventSetup & iSetup)
}
std::vector<std::string> levels = jetCorrs[0][jetRef].correctionLabels();
if(std::find(levels.begin(), levels.end(), "L2L3Residual")!=levels.end()){
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("L2L3Residual"), JetCorrFactors::NONE, setindex);
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("L2L3Residual"));
}
else if(std::find(levels.begin(), levels.end(), "L3Absolute")!=levels.end()){
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("L3Absolute"), JetCorrFactors::NONE, setindex);
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("L3Absolute"));
}
else{
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("Uncorrected"), JetCorrFactors::NONE, setindex);
ajet.initializeJEC(jetCorrs[0][jetRef].jecLevel("Uncorrected"));
if(printWarning_){
edm::LogWarning("L3Absolute not found") << "L2L3Residual and L3Absolute are not part of the jetCorrFactors\n"
<< "of module " << jetCorrs[0][jetRef].jecSet() << ". Jets will remain"
Expand Down
25 changes: 19 additions & 6 deletions PhysicsTools/PatAlgos/plugins/PATMETProducer.cc
Expand Up @@ -47,11 +47,15 @@ PATMETProducer::PATMETProducer(const edm::ParameterSet & iConfig):
if(calculateMETSignificance_)
{
metSigAlgo_ = new metsig::METSignificance(iConfig);
jetToken_ = mayConsume<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("srcJets"));
pfCandToken_ = mayConsume<edm::View<reco::Candidate> >(iConfig.getParameter<edm::InputTag>("srcPFCands"));
rhoToken_ = consumes<double>(iConfig.getParameter<edm::InputTag>("srcRho"));
jetSFType_ = iConfig.getParameter<std::string>("srcJetSF");
jetResPtType_ = iConfig.getParameter<std::string>("srcJetResPt");
jetResPhiType_ = iConfig.getParameter<std::string>("srcJetResPhi");
jetToken_ = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("srcJets"));
pfCandToken_ = consumes<edm::View<reco::Candidate> >(iConfig.getParameter<edm::InputTag>("srcPFCands"));
std::vector<edm::InputTag> srcLeptonsTags = iConfig.getParameter< std::vector<edm::InputTag> >("srcLeptons");
for(std::vector<edm::InputTag>::const_iterator it=srcLeptonsTags.begin();it!=srcLeptonsTags.end();it++) {
lepTokens_.push_back( mayConsume<edm::View<reco::Candidate> >( *it ) );
lepTokens_.push_back( consumes<edm::View<reco::Candidate> >( *it ) );
}
}

Expand Down Expand Up @@ -93,7 +97,7 @@ void PATMETProducer::produce(edm::Event & iEvent, const edm::EventSetup & iSetup

//add the MET significance
if(calculateMETSignificance_) {
const reco::METCovMatrix& sigcov = getMETCovMatrix(iEvent);
const reco::METCovMatrix& sigcov = getMETCovMatrix(iEvent, iSetup);
amet.setSignificanceMatrix(sigcov);
double metSig=metSigAlgo_->getSignificance(sigcov, amet);
amet.setMETSignificance(metSig);
Expand Down Expand Up @@ -160,7 +164,7 @@ void PATMETProducer::fillDescriptions(edm::ConfigurationDescriptions & descripti
}

const reco::METCovMatrix
PATMETProducer::getMETCovMatrix(const edm::Event& event) const {
PATMETProducer::getMETCovMatrix(const edm::Event& event, const edm::EventSetup& iSetup) const {
std::vector< edm::Handle<reco::CandidateView> > leptons;
for ( std::vector<edm::EDGetTokenT<edm::View<reco::Candidate> > >::const_iterator srcLeptons_i = lepTokens_.begin();
srcLeptons_i != lepTokens_.end(); ++srcLeptons_i ) {
Expand All @@ -176,8 +180,17 @@ PATMETProducer::getMETCovMatrix(const edm::Event& event) const {
edm::Handle<edm::View<reco::Candidate> > inputCands;
event.getByToken( pfCandToken_, inputCands );

edm::Handle<double> rho;
event.getByToken(rhoToken_, rho);

JME::JetResolution resPtObj = JME::JetResolution::get(iSetup, jetResPtType_);
JME::JetResolution resPhiObj = JME::JetResolution::get(iSetup, jetResPhiType_);
JME::JetResolutionScaleFactor resSFObj = JME::JetResolutionScaleFactor::get(iSetup, jetSFType_);

//Compute the covariance matrix and fill it
reco::METCovMatrix cov = metSigAlgo_->getCovariance( *inputJets, leptons, *inputCands);
reco::METCovMatrix cov = metSigAlgo_->getCovariance( *inputJets, leptons, *inputCands,
*rho, resPtObj, resPhiObj, resSFObj, event.isRealData());

return cov;
}

Expand Down
7 changes: 6 additions & 1 deletion PhysicsTools/PatAlgos/plugins/PATMETProducer.h
Expand Up @@ -72,8 +72,13 @@ namespace pat {
edm::EDGetTokenT<edm::View<reco::Jet> > jetToken_;
edm::EDGetTokenT<edm::View<reco::Candidate> > pfCandToken_;
std::vector< edm::EDGetTokenT<edm::View<reco::Candidate> > > lepTokens_;
edm::EDGetTokenT<double> rhoToken_;
std::string jetResPtType_;
std::string jetResPhiType_;
std::string jetSFType_;

const reco::METCovMatrix getMETCovMatrix(const edm::Event& event) const;
const reco::METCovMatrix getMETCovMatrix(const edm::Event& event,
const edm::EventSetup& iSetup) const;

};

Expand Down
Expand Up @@ -52,6 +52,10 @@
srcJets = cms.InputTag("selectedPatJets"),
srcPFCands = cms.InputTag("packedPFCandidates"),
srcLeptons = cms.VInputTag("selectedPatElectrons", "selectedPatMuons", "selectedPatPhotons"),
srcJetSF = cms.string('AK4PFchs'),
srcJetResPt = cms.string('AK4PFchs_pt'),
srcJetResPhi = cms.string('AK4PFchs_phi'),
srcRho = cms.InputTag('fixedGridRhoAll'),
parameters = METSignificanceParams
)

Expand Down

0 comments on commit 6c4b6cd

Please sign in to comment.